Modeling the 20th-century distribution changes of Microgyne trifurcata, a rare plant of the southern South American grasslands

Microgyne trifurcata is a rare native plant species from one of the areas with the highest human impact on the environment in southern South America. Its habitat, mostly grasslands suitable for agriculture, has been increasingly covered by crops since the late 1800s. Microgyne trifurcata provides an excellent case study to understand how different environmental variables have affected the distribution area of a rare species. This study aims to estimate the impact of topoclimatic and land-use changes in the distribution of Microgyne trifurcata throughout the twentieth century. We carried out recent past and present distribution modeling using the Ensembles of Small Models (ESM) methodology. In this spatio-temporal study, we included climatic, topographic, and land-use variables. We classified the occurrences into two periods of the twentieth century. The first dates from 1901 to 1940, and the second, from 1960 to 2000, when the main cropping changes of the area occurred. The projected area between 1960 and 2000 provides for this species new suitable habitats toward the northeast of the area of study. Our results highlight the importance of assessing the combined impacts of climate and land-use changes on species distributions over time. This study shows that the potential area of Microgyne trifurcata decreased and underwent fragmentation throughout the twentieth century when these variables combined are used to model its distribution. Our outcomes prompt future studies on the vulnerability of Microgyne trifurcata to outline conservation strategies.


Introduction
The extensive cropping on the native grasslands of southern South America, especially in Argentina, began about 1875, and since then, wheat, maize, sunflower, and, lately, soybean crops have expanded continuously, covering most of the areas suitable for agriculture (Bilenca and Miñarro 2004;Magrín et al. 2005). This activity has had a dramatic impact on the native grasslands by fragmenting, reducing, and modifying the original habitats (Andrade et al. 2015;Baldi and Paruelo 2008;Baldi et al. 2006;Bilenca et al. 2012;Guerschman and Paruelo 2005;Staude et al. 2018).
Like cropping, climate change is among the most significant factors threatening biodiversity (Newbold 2018). Southern South America, in particular, faced changes in climatic patterns during the twentieth century. For instance, the grasslands of this region experienced changes in rainfall patterns during the second half of the twentieth century (Berbery et al. 2006). An increase in late spring and summer rainfall (Minetti et al. 2003) caused a westward displacement of the transition area to semi-arid regions (i.e., Monte biogeographic province; Cabrera and Willink 1973), these regions being the boundary of rainfed agriculture (Berbery et al. 2006;Magrín et al. 2005). Concurrent with these climate changes, radical variations in the distribution of rangelands and grazing lands, together with increasing crop production, affected the grasslands of southern South America in the last century (Magrín et al. 2005). Thus, pastures and annual forage crops were progressively replaced by wheat-soybean relay cropping and maize and sunflower crops. Simultaneously, the animal stocking rate increased in those lands with less agricultural aptitude (Abba et al. 2015;Paruelo et al. 2006;Viglizzo et al. 2011).
Some studies approached the consequences of both land-use and climate changes on the southern South American grasslands (e.g., Bilenca et al. 2009;Viglizzo and Frank 2006). Others (Baldi and Paruelo 2008;Baldi et al. 2006;Bilenca et al. 2012) focused on native grassland fragmentation due to an increase in cropping. However, up to the present, there have been no studies integrating land-use and climate changes in the last century applying methodologies such as Species Distribution Modeling (SDM).
Species with particular habitat requirements and distributions that also show an accurate collecting record over time could help decode the effects of increased land-use and climate changes on range dynamics (Yang et al., 2018;Zhang et al., 2020). Microgyne trifurcata Less. (1832) is a species of Asteraceae endemic to central-eastern Argentina, southern Brazil, and Uruguay (Heiden and Sancho 2016;Sancho 2014;Sancho and Ariza Espinar 2003;Sancho et al. 2006). Most of the range of M. trifurcata involves the so-called Pampas region, which is dominated by natural grasslands and is one of the areas with the highest human impact on the environment of southern South America (e.g., Magrín et al. 2005;Miñarro and Bilenca 2008). According to our field observations carried out since 2009, M. trifurcata is apparently becoming more restricted as land use in the area increases. In this context of climate and land-use changes, M. trifurcata is an excellent case study to understand, by means of applying recent past and present distribution modeling, how these changes have impacted its distribution.
Microgyne trifurcata is a perennial herb or small sub-shrub with glandular leaves, distinctly trifid apically, solitary capitula, and marginal white corollas that occurs in dry, clayey, calcareous and rocky soils, and roadsides up to 1300 m a.s.l. (Sancho et al. 2006) (Fig. 1a-f). The populations of M. trifurcata are rather small and usually consist of a few sparsely arranged individuals where specific local environmental characteristics are met. Due to these peculiar population features, we consider M. trifurcata a rare species. Rare species, either naturally or as a result of the impact of human-led activities, are characterized by restricted geographic ranges, habitat specialization (i.e., species that occur only in one or a few specific habitat types), and small population sizes (Aitken et al. 2007;Medrano and Herrera 2008;Rabinowitz 1981). Datasets on rare species distribution are usually characterized by a small number of observations often gathered over relatively long periods of time (e.g., throughout the twentieth century), limited spatial accuracy, and a lack of valid absences (Wisz and Guisan 2009). However, a small number of available observations can frequently provide a reasonably comprehensive view of the species' distribution (Lomba et al. 2010). For specific cases such as those of Microgyne trifurcata, the use of Ensembles of Small Models (ESMs) provides a novel strategy by predicting the distribution of rare species (Lomba et al. 2010;Scherrer et al. 2019). The smaller the sample size of a species is, the more significant is the gain in ESMs performance compared to standard models (Breiner et al. 2018). On the other hand, discrepancies between different species distribution models (SDMs) can be significant, making the selection of the most appropriate model difficult (Cuyckens et al. 2016;Renner and Warton 2013). In this context, ensemble forecasting approaches may be a helpful tool (Cianfrani et al. 2011;Hao et al. 2019).
Historical ecological data can help ecologists integrate future perspectives within historical contexts in order to provide insights into the long-term  (Beaugrand et al. 2015;Scherrer et al. 2017). However, data limitations, including incomplete and spatially biased datasets must be taken into consideration when long-term ecological data (LTE: spanning [ 50 years), mainly historical data, are included (Boakes et al. 2010;McClenachan et al. 2012). Besides, multiple resources and species distribution models have proven to be effective methods for reconstructing historical distributions from LTE data (Zhang et al. 2017). Within this framework, Microgyne trifurcata, a species with historical records dating back to 1878, provides an accurate time frame for applying SDM methods.
This study aims to estimate the impact of climate and land-use changes in the distribution of Microgyne trifurcata throughout the twentieth century. Specifically, we carry out recent past and present SDM applying the ESMs framework, first including only topoclimatic variables and second, topoclimatic and land-use data.
We expect our results to show that changes in climate and land-use were reflected in changes in the potential distribution area of Microgyne trifurcata throughout the last century. Finally, we presume that our outcomes will provide information on the vulnerability of Microgyne trifurcata useful to outline conservation strategies in the future.

Study area
The study area corresponds to the currently known distribution of Microgyne trifurcata and includes central and eastern Argentina, southeastern Brazil, and Uruguay (-67°to -50°Longitude West and -27°to -41°Latitude South; Fig. 2a).
Microgyne trifurcata occurs in areas usually sparsely covered by natural grasslands or, less commonly, in clearings of xerophytic native forests. Most of these areas are known as Río de la Plata grasslands Laterra and Rivas 2005;Miñarro and Bilenca 2008;Soriano et al. 1992) or Pampas region (e.g., Magrín et al. 2005). The study area is mostly a vast plain interrupted only by scattered hill ranges, rarely reaching more than 900 m a.s.l. On the western range of the species distribution, the area includes higher elevations reaching up to 1300 m a.s.l.
(the highest elevation known for Microgyne trifurcata). The climate in the area is temperate, with mean temperatures varying between 14°C in the south and 18°C in the north. Annual rainfall decreases from 1300 mm in the northeast to 500 mm in the southwest, although interannual variability of rainfall is quite frequent in the Pampas with extensive rainfall or drought (Bilenca and Miñarro 2004).
From a biogeographic point of view, the study area mostly agrees with the Pampeana province and partly with the Espinal province and the Chaqueña province (Cabrera and Willink 1973). From a conservation perspective, these grassland ecosystems are under considerable environmental pressure, and to make matters worse they are very poorly represented in the Systems of Protected Areas of the Southern Cone countries, a situation that dramatically increases its vulnerability (Hoekstra et al. 2005). Only ca. 1.5% of the grasslands in the Pampas region in Argentina (Miñarro and Bilenca 2008), 0.21% of those in Uruguay, and 0.36% of those in Rio Grande do Sul (Brazil) are included in protected areas under sustainable use programs (Baldi et al. 2017;Bilenca and Miñarro 2004).

Occurrences data
We obtained historical occurrence data from specimens of five herbaria (BAB, CORD, LP, MO, MVFA, SI; Thiers, continuously updated), field observations, and specific literature (Sancho et al. 2006;Sancho 2014). We completed a database of 116 records of Microgyne trifurcata. These records span from 1878 to 2018. To carry out a long-term comparison, we classified all of the records into two periods according to the available variables' databases, one at the beginning of the twentieth century spanning from 1901 to 1940 and a second at the end of the twentieth century, from 1960-2000.
We defined these two periods based on 1. the high impact of the cropping boom on the study area in the last century, especially from the 1970s (Magrín et al. 2005) and 2. the climate changes undergone in the study area, especially since the early 1970s (Barros et al. 2015;Hoffmann 1989). In order to highlight the differences among variables within each period, we left out from the analyses the records between 1940 and 1960. Duplicated occurrences were removed within each time period dataset, and those finally selected were filtered to avoid spatial autocorrelation using the Species Occurrences Disaggregation function of 'ecospat' package (Di Cola et al. 2017) in R 3.4.0 (R Core Team 2017). This procedure removes species occurrences which are closer to each other than a specified distance threshold, in this case 10 km (Slodowicz et al. 2018). The final database included a total of 74 occurrences: 35 for the first period  and 39 for the second period   (Fig. 2b).

Environmental data
The set of variables for modeling the distribution of Microgyne trifurcata included topoclimatic and landuse variables.
The final set of topoclimatic variables covered one topographic variable and four climatic variables. We obtained the data of climatic variables from the monthly climatic variables of the CRU-TS v3.10.01 Historic Climate Database for GIS (Climatic Research Unit Time-Series 2012, University of East Anglia) over the last century. CRU-TS provides high-resolution gridded datasets for a considerable time period  global land data for multiple variables) which is appropriate and fits well our data in order to answer the questions proposed in this study.
We selected the topographic variable Openness index (generated from DEM GMTED2010, Danielson and Gesch 2011; with SAGA software, Conrad et al. 2015) for the analysis because it fits the environmental characteristics of Microgyne trifurcata (Sancho 2014;Sancho et al. 2006).
The land-use variables were obtained from the land-use harmonization project (Hurtt et al. 2011). We used the LUH2 v2_high historical land-use projection dataset, which consists of estimates of the proportion of each terrestrial grid cell in each of six major landuse classes: primary vegetation, secondary vegetation, plantation forest, cropland, pasture, and urban land.
The grid cell resolution for climatic, topographic, and land-use variables was initially 0.5°latitude 9 0.5°longitude.
Regarding the climatic variables, data of two time periods (1901 to 1940 and 1960 to 2000) were extracted for the study area from monthly climatic variables of CRU-TS v3.10.01 dataset. Preliminarily, we obtained 19 bioclimatic variables for each period with 'dismo' package ) in R (R Development Core Team, Vienna, Austria), using the  ; red dots indicate occurrences in the second period  monthly temperature and precipitation variables. Mean values of the bioclimatic variables for each period (1901-1940 and 1960-2000) were used to account for interannual climate variation and potential lags in species' responses to climate (De Frenne et al. 2013).
Due to the static nature of topography throughout the time frame covered in our study, we used the same raster layer, Openness index, to model the potential distribution in both periods.
As done for climatic variables, we selected three land-use variables significant to the study area: pasture area (lands primarily used for domesticated forage plants for livestock), rangeland (areas with primarily native vegetation suitable for grazing use), and C3 annuals crop area (e.g., Miñarro and Bilenca 2008;Viglizzo et al. 2011Viglizzo et al. ) from a dataset covering 1901Viglizzo et al. to 1940Viglizzo et al. and 1960Viglizzo et al. to 2000 We first assessed the correlation among the topoclimatic variables by applying a pairwise correlation analysis (Pearson's product moment correlation coefficient). The outcomes allowed discrimination of those variables with high correlation ([ 0.8;Sokal and Rohlf 1979). Based on these results, the final set of five variables used to model the potential distribution of Microgyne trifurcata were four climate variables (temperature of the wettest quarter, Bio8; precipitation of wettest quarter, Bio16; precipitation of driest quarter, Bio17; potential evapotranspiration) and one topographic variable (Openness index). The selected variables are reliable descriptors of the climate change observed during the twentieth century in the study area. The variables related to precipitation (Bio16 and Bio17) as well as evapotranspiration are those that have shown greater variation (seasonal and annual increase) from the beginning of the twentieth century to the present time (Barros et al. 2015). In addition, they have been used to account for the past (Minoli et al. 2018;Rodrigues Silva et al. 2018) and future climate changes (Ferretti et al. 2018;Minoli et al. 2018) in the study area.
With the same method used for topoclimatic variables, the correlation among the final topoclimatic and land-use variables was, in addition, assessed to avoid multicollinearity between these sets of variables. The results showed these variables to be uncorrelated.
We used the same set of variables (topoclimatic ? land use) for both models 1901-1940 and 1960-2000. The data were resampled to a 1 km pixel resolution using bilinear interpolation with 'raster' R package (Hijmans 2017). To assess the changes over time that could impact on the distribution of Microgyne trifurcata, we calculated the differences between the mean values of the variables from the first  and second  time periods.

Modeling framework
To obtain the potential distributions of Microgyne trifurcata, we developed two sets of models. One set was based on topoclimatic variables, while a second set was based on topoclimatic variables ? land-use variables. Both sets of models were calculated for each period of the twentieth century. The first ranged from 1901 to 1940, and second, from 1960 to 2000.
As Microgyne trifurcata is rare, and its populations small and scattered throughout the distribution range, the occurrence data for this species are relatively low. The combination of limited occurrence data and a large number of predictors in SDM can lead to model over fitting (Breiner et al. 2015). To overcome this bias, we followed the ensemble of small models (ESMs) framework (Breiner et al. 2015), which is suitable for modeling rare species (Breiner et al. 2015(Breiner et al. , 2018, available in the 'ecospat 2.2.0' R package (Di Cola et al. 2017). The ESMs framework involves building many models with small subsets of predictors, two as implemented here (i.e., bivariate models). Each of these small models is then evaluated according to a cross-validated model performance measure. In the last step, all small models are averaged into a weighted ensemble, where the weights are based on the small model performances, to produce the final ESMs (Breiner et al. 2015(Breiner et al. , 2018. We used four different modeling algorithms for our study: Artificial Neural Networks (ANN), Generalized Boosting Model (GBM), Generalized Linear Model (GLM), and Maximum Entropy (Maxent, Max. P as implemented by Phillips et al. 2006) in an ensemble modeling scheme. These algorithms were among the best performing techniques in ESMs (Breiner et al. 2018) and represented the two main families of models (regression-based and machine-learning methods). The ensemble forecasting integrates the results of multiple SDM algorithms into a single geographical projection for each time period of the twentieth century considered in the analysis of Microgyne trifurcata. This approach reduces the uncertainties associated with a single-model algorithm (Araújo and New 2007;Araújo et al. 2019).
The model's predictive performance was assessed using the area under the curve (AUC) of the receiver operating characteristic (ROC) (Phillips et al. 2006) and TSS (true skill statistics). Weights of the bivariate models with an AUC lower than 0.5, hence worse than random (counter-predictions), were set to 0. We calibrated ESMs by fitting numerous bivariate models, which were then averaged into an ensemble model using weights based on model performances. For each ESMs, we ran 100 split-samples (70% calibration/ 30% evaluation) using one set of 1000 pseudoabsences created at random (Elith et al. 2006) in the entire study area. Finally, an ensemble prediction was created by averaging the 10 runs weighted by AUC (Fielding and Bell 1997).
To estimate the surface area of suitable habitats for Microgyne trifurcata for each period, the continuous maps of species distribution generated by our ESMs were converted into binary maps (suitable/unsuitable) by applying a threshold and discarding the lowest 10% of training locations, hence encompassing 90% of training presences (minimal predicted area, MPA90; Engler et al. 2004). The size (in km 2 ) of the suitable areas obtained for both periods and both groups of variables was calculated and compared.

Model transferability and niche overlap
Two tests were carried out to assess the suitability of our data for predicting the temporal projections for Microgyne trifurcata: temporal transferability and niche overlap.
Temporal transferability: The temporal transferability of the model is the ability to predict species distribution in a different time frame (Regos et al. 2019). In our case, the model calibrated on the first period  should be able to predict the distribution of the second period  simply by projecting it with the new climate data and vice versa. In order to carry out the transferability test, we used the same settings described above for the ESMs. We used the AUC scores to assess the performance of models during the cross-validation procedure as these values have been used as accuracy metrics to evaluate model transferability (Dobrowski et al. 2011;Fielding and Bell 1997).
Niche overlap: The incongruence of environmental niches of Microgyne trifurcata modeled in two periods could reflect either environmental niche shifts during the last 100 years (highly unlikely) or biased datasets. To avoid these issues, we carried out the Niche Overlap test implemented in the new R package, 'humboldt' (Brown and Carnaval 2019).

Temporal transferability
All the models performed well during the crossvalidation procedure. Similar values of AUC for models calibrated with data of the first period (0.77) were obtained for models projected with data from the second period (0.81).

Niche overlap
We obtained a high value of analogous climate space percentage (83.6%) and a low value of Niche Overlap index (D = 0.231), which indicates that the two sets of occurrences do not overlap with each other in the available climate space (see Online Figs. 1, 2).

Comparison of environmental variables between the two periods
The differences between the first  and second  periods showed changes over time for the variables included in the analyses (see Online Table 1 and Online Figs. 3, 4). The temperature values of the wettest quarter increased in the central, northeastern, and southwestern parts of the study area during the second period. For the same variable, the temperature values decreased toward the northwestern part of the study area during the same period (see Online Fig. 3a). An overall increase was observed in the precipitation values of the wettest quarter, with the highest values in the north-central part of the study area (see Online Fig. 3b). Likewise, the precipitation values of the driest quarter increased during the second period toward the eastern part of the study area, less so toward the southwestern part and, instead, the values decreased in the central part and toward the northwestern part of the study area (see Online Fig. 3c). Also, the potential evapotranspiration values increased eastward and decreased toward the southwestern part of the study area (see Online Fig. 3d).
Regarding the differences between the mean values of land-use variables for the first and second periods, an increase of the area of annual crops of C3 plants was verified mainly in the central and southern parts of the study area (see Online Fig. 4a). Also, the results showed that the pasture areas increased mainly in central Uruguay and south of Buenos Aires Province and decreased in the central part of the study area and the surroundings of the Uruguay and Río de La Plata rivers (which define the boundaries between Argentina and Uruguay) throughout the two periods (see Online Fig. 4b). Finally, the rangelands expanded during the second period mainly in Uruguay and south of Buenos Aires Province (see Online Fig. 4c).

Species distribution modeling analyses
The AUC values for ESMs using topoclimatic variables only varied from 0.78 to 0.83 for all algorithms and for the ensemble prediction, thus showing a good performance. For ESMs using topoclimatic ? landuse variables, the AUC values were 0.72-0.86. The TSS values for all algorithms and for the ensemble prediction were 0.5-0.57 for ESMs using topoclimatic variables only and 0.50-0.65 for ESMs using topoclimatic ? land-use variables.
The importance of the variables used to obtain the potential distribution models was different according to both the period analyzed and the variables used in the analyses (i.e., topoclimatic variables or topoclimatic ? land-use variables) (Fig. 3). The temperature of the wettest quarter (Bio8) was the most relevant variable for the ESMs of the second period , obtained either with topoclimatic or topoclimatic ? land-use variables. Regarding the ESMs of the first period , in the models obtained with topoclimatic variables only, precipitation of wettest quarter (Bio16) was the most significant variable. On the other hand, C3 annuals crop area was the most relevant variable in models obtained with topoclimatic ? land-use variables.
The suitable area for Microgyne trifurcata provided by ESMs based on the topoclimatic variables in the first period  resulted in a relatively discontinuous set of areas (Fig. 4a) including Northeastern San Luis, Central Cordoba, Southeastern Santa Fe, northeastern parts of the study area, and areas associated to Ventania mountain ranges in Buenos Aires, Southwestern and Eastern Entre Ríos, and Southwestern Uruguay.
For the second period , the analyses provided a suitable area overall agreeing with that of the first period (Fig. 4a), although the general limits greatly expanded eastward to Central Entre Ríos in Argentina, and Central and Western Uruguay. The suitable area in central Córdoba and San Luis obtained for the second period decreased more than half the size with respect to the area obtained in the first period, somewhat restricting Microgyne trifurcata to the Sierras Pampeanas environments (as described by Crisci et al. 2001).
The estimated suitable area obtained for the second period by including only the topoclimatic variables was larger (about a 22%) than that of the first period (154,185 km 2 and 188,225 km 2 for the first and second period, respectively) ( Table 1).
In the analyses using topoclimatic ? land-use variables, for the first period, the suitable area obtained by ESMs yielded similar results to those obtained by including only the topoclimatic variables, although with some modifications (Fig. 4b). For instance, the center of the potential area expanded more widely toward Southern Santa Fe and Northeastern Buenos Aires.
Regarding the second period, the projected models by ESMs show similar overall results to those using only the topoclimatic variables in the second period (Fig. 4b). Thus, the projected area expands eastward to Central Entre Ríos in Argentina, and Western and Central Uruguay. Also, the projected area for the second period no longer includes Ventania mountain ranges. Besides, the projected area provided by ESMs drastically decreases in Central Córdoba, Southern Santa Fe, Western Entre Ríos, and Northeastern Buenos Aires.
In the analyses including the topoclimatic ? landuse variables, for the second period , the area was smaller (about 25%) than that obtained in the first period (115,812 km 2 and 86,559 km 2 for the first and second period, respectively).
Overall, similarly to the results obtained applying topoclimatic variables, the models obtained by adding the land-use variables provided and increased area in central Uruguay and surroundings of Uruguay River, indicating suitable environments toward the northeastern part of the study area. However, the projected area obtained for the second period with topoclimatic ? land-use variables drastically decreases (54%) if compared with that obtained for the same period using only topoclimatic variables (Fig. 5a, b).

Discussion
As we expected, the results of this study show that changes in climate and land-use were reflected in changes of potential distribution area of Microgyne trifurcata throughout the last century. Our results also highlight the importance of assessing the combined impacts over time of climate and land-use changes on species distributions. Methodologically, the ESM strategy provided models with good performance to assess the potential distribution changes of Microgyne trifurcata. Regarding land-use variables, modifications over time of pasture and rangelands areas did not impact the projected area of this species as much as the increment in cropping areas.

Implications of climate changes in the projected area
Climate change over the last century in the study area could partly explain the modifications in the suitable projected area for Microgyne trifurcata over the analyzed time. In Argentina, climate studies show that during the twentieth century, Buenos Aires province faced a continuous increase in precipitation and temperature (Barros et al. 2015;Berbery et al. 2006;Menéndez 2006) that affected the Pampean region (Magrín et al. 2005). Precipitation increased from 1916 to 1991 in the study area. From the 1950s, the increment also was registered throughout most of the country (Grimm et al. 2000). Some sites of the study area showed a pattern of increasing precipitation and minimum temperatures and decreasing maximum temperatures (Messina 1999). Moreover, since the early 1970s, westward shifts of the isohyets in the study area resulted in considerably increments in the mean annual precipitation in the sub-humid and semiarid zones of the Pampean region (Hoffmann 1989). Our results show similar climate patterns in, for Fig. 3 Importance of topoclimatic and land-use variables used for modeling the potential distribution of Microgyne trifurcata. a Topoclimatic variables for the first period . b Topoclimatic variables of the second period . c Topoclimatic ? land-use variables of the first period . d Topoclimatic ? land-use variables of second period . Bio8 mean temperature of wettest quarter, Bio16 precipitation of wettest quarter, Bio17 precipitation of the driest quarter, ET potential evapotranspiration, RL rangelands, PT pastures, C3 C3 annuals crops instance, the increase of precipitation in of wettest quarter, mean temperatures of the wettest quarter, and potential evapotranspiration (see Online Fig. 3a, b, and d). The suitable area projected for Microgyne trifurcata agrees with these climate changes over time by providing new available environments toward the northeastern part of the study area in Central Uruguay and along the Uruguay River (Fig. 4). Apparently, the mean temperature of the wettest quarter had an outstanding effect on Microgyne trifurcata by somewhat establishing new eastside boundaries for the projected area (see Online Fig. 3a).

Implications of land-use changes in the projected area
Our results for Microgyne trifurcata show that the increase of C3 annual plant crops such as those that occurred in central Argentina and Santa Fe and Buenos Aires provinces from the 1970s (see Online Fig. 4a) did not promote suitable environments for this species. On the contrary, Microgyne trifurcata has been usually found in rangelands (Sancho et al. 2016) whose expansion (see Online Fig. 4c) has not drastically affected the distribution of this species. Despite the apparent negative effect of cattle grazing, ecosystems under this environmental pressure represent suitable environments for Microgyne trifurcata. Our field observations through time on Microgyne trifurcata populations confirmed that this species is missing in annual crop areas but found in grazed areas. Our outcomes agree with these observations by showing an increase in the projected suitable area in rangelands from the 1960s up to the present (see Online Fig. 4c).
Simultaneously with the climate changes, drastic modifications in crop production and distribution and land-use were also registered for the study area (Bilenca and Miñarro 2004;Miñarro and Bilenca 2008). Nowadays, Argentina's Pampean region accounts for more than 90% of national grain production, with soybean, wheat, maize, and sunflower as the main crops (Magrín et al. 2005). Both, crop production and the relative importance of crops have changed with time (Magrín et al. 2005), particularly during the last 30 years (some of these changes can be observed in our Online Fig. 4a). From the 1970s, the soybean expansion increased continuously in both planted areas and yields (Hall et al. 1992;Miñarro and Bilenca 2008) covering most of the agricultural lands and even roadsides.
In Argentina, the most significant changes in the annual crop areas occurred in Buenos Aires, Cordoba, Entre Ríos, and Santa Fe provinces (Magrín et al. 2005). In the case of Uruguay, the total crop area increased in most departments, even though the area devoted to annual crops did not change or slightly decreased over the last half of the twentieth century (Miñarro and Bilenca 2008). The relatively low expansion of annual croplands in Uruguay contrasts with the high increase of forest plantation, especially eucalyptus and pines (Achkar et al. 2012;Tiscornia et al. 2014). These areas are mostly located in western and central Uruguay. All the above indicates that the expansion of soybean planted areas in Argentina and forest plantations in Uruguay represents the most important transformation affecting natural pastures of these countries in the recent decades (Achkar et al. 2012;Brazeiro et al. 2008;Paruelo et al. 2006).
Studies on changes in the distribution of species due to climate and land-use variations in the Pampas region are scarce and have been focused, for the most part, on animal species distributions (e.g., Abba et al. 2015;Bilenca et al. 2017;Ferretti et al. 2018;Guerrero 2014;Medone et al. 2015;Minoli et al. 2018;Nori et al. 2013). Most of the species treated by these authors are distributed in subtropical forests of Northeastern Argentina, Eastern Paraguay, and Southern Brazil and show southward expansion of their distributional boundaries throughout the twentieth century. The available studies involving species from the Pampas and northern Patagonia in Argentina, b Fig. 4 ESMs obtained with the two groups of variables applied to Microgyne trifurcata. a Models obtained with topoclimatic variables only. b Models obtained with topoclimatic ? land-use variables. References are indicated in each map Table 1 Suitable area (km 2 ) obtained with ESMs for each period (1901-1940; 1960-2000)  however, show restriction in species distribution, probably as a response to human disturbance and/or recent climate change (Agnolin and Lucero 2014;Teta et al. 2014). Unlike the above studies, our work on this rare species of the Pampean region projects new suitable habitats to the northeastern parts of the study area along the Uruguay River and in Central Uruguay ( Fig. 4a and b). Our study is a starting point to elucidate if species that occur in similar environments to those of Microgyne trifurcata that underwent the same climate and land-use changes have experienced similar distributional patterns over time. As the driving forces of distributional shifts of the species of the Pampean region are highly complex, the consequences of climate and land-use changes over habitat availability should be analyzed specifically. Areas like those of Northern Buenos Aires and Southern Santa Fe provinces, mainly devoted to agriculture, were not available in the second half of the twentieth century in the projected models that included the land-use variables (Fig. 4b).
The results for Microgyne trifurcata indicate that ESMs modeling framework underlines the effects of land-use changes through time, supporting the importance of assessing the combined impacts of climate and land-use changes on species distribution. The extent to which the factors that shaped the distribution of Microgyne trifurcata have impacted other taxa occurring in the same area depends on multiple variables. One of the factors influencing species spatial distribution is their dispersal capabilities. In the case of Microgyne trifurcata, the presence of pappus, the typical dispersal structure of Asteraceae, potentially assures the wind dispersal of the fruits (Funk et al. 2009). However, the capability of fruits to get established in the ground is a different matter, which in the case of Microgyne trifurcata, needs to be further investigated.
Our results show a decrease in the potential area of Microgyne trifurcata throughout the last century if topoclimatic ? land-use variables are included in the analysis (the projected area is 25% smaller in the second period). This together with our field observations on the same populations since 2009, suggest recent habitat restriction for the same area. Populations once more widely distributed (e.g., those from Espinal forest in Western Uruguay) are now restricted to suitable microhabitats like rocky grasslands, that are often surrounded by increasing soy fields. Indeed, a decrease in the number of individuals per population was verified during fieldwork, but these impressions remain subjective and need to be quantified.

Future perspectives
As explained above, some of the limitations of SMD rely on a significant number of occurrences to obtain reliable models. The ESM tool deals with those constraints by predicting the distribution of rare and most threatened species (Breiner et al. 2015). Our study using this methodology has shown that ESM could be applied to other rare species of Rio de La Plata grasslands allowing those most vulnerable to climate and land-use changes to be identified.
Our results, together with ongoing phylogeographic studies (Viera Barreto et al. in prep.) devoted to understand the evolution of Microgyne trifurcata distribution, on one hand, and the recent variation in population sizes on the other, will surely assist on the planning of conservation strategies for this species.
The Rio de La Plata grasslands suffered significantly from climate change during the last century. In addition, they were subjected to considerable pressure from human-led activities such as agriculture and forest plantations. The studies involving these matters in this area have approached the land-use issue from different perspectives (e.g., Bilenca et al. 2009;Viglizzo and Frank 2006). Our study shows the importance of carrying out integral analyses of species distribution changes by including both climate changes and human pressure in a time frame, from the SDM perspective. This approach could be the starting point for similar integrative studies involving species distribution changes in this area over the last century as a way to understand how future changes in climatic variables or use of land will impact their current distribution. comments on the manuscript. This research was funded by