Assessing climate change risks to the geographical distribution of grass species

ABSTRACT

There is extensive evidence showing the impact of climate change on the biology and biogeography of species. Adopting drought-tolerant plants to conserve water is a potential adaptation to reduce the consequences of climate change. Accordingly, it was hypothesized that climate change would not affect potential distributions of drought-tolerant species. Here, this hypothesis was tested to model the potential distribution of three drought-resistant plant. Here, the potential distribution of Agropyron cristatum, Agropyron desertorum, and Festuca arundinacea was studied in Iran under current and future climate conditions, using 10 species distribution models. Sixty-two climate change scenarios (19 global climate models (GCMs) under four Representative Concentration Pathways (RCPs)) were used to model the potential distribution of the three plants in Iran in the future. The three species have different responses to predicted climate change due to species-environment interactions, species morphological and physiological advancement. The three species showed different responses to predicted climate change due to species-environment interactions. Festuca arundinacea and Agropyron cristatum will, respectively, experience the most and least severe decline in suitable habitats in the next 50 years. This result is because decreased annual precipitation caused an increase in habitat suitability for A. cristatum, while the same variable had the opposite effect for A. desertorum and F. arundinacea. On the other hand, F. arundinacea grows on moist soils that decreased annual precipitation caused a decrease in habitat suitability. Also, our results have clearly shown that plant species drought-stress tolerant are not immune to climate change and their current distributions undergo significant changes as a result of the changing of climate.

KEYWORDS:

• Distribution
• SDM
• climate change
• agropyron cristatum
• agropyron desertorum
• festuca arundinacea

Introduction

There is extensive evidence showing the impact of climate change on the biology and biogeography of species [e.g., 1, 2, 3]. However, predicting how climate change can impact individual species has been challenging.⁴ As a result of climate change, the current habitats of a large number of species will no longer be able to sustain their presence; consequently, species might need to migrate to suitable climates to stand a chance of survival .^2,5–7

There is evidence of species reacting to climate change by migration. In particular, there are extensive data on plant migrations as a result of the past⁸ and present^9–11 changes in climate. Plants seem to be responding to the increased availability of water in lower elevation as rain falls more often.¹² Such observations indicate the complexity of plants response to climate change, and suggest the response might differ for different latitudes, longitudes, elevations, and biomes. Simultaneously, habitat degradation and fragmentation, and invasive species undermine ecosystem resilience and impede the migration of species.¹³ Hardly any region is exempt from the impacts of climate change. The rapid pace of climate change might not allow species to adapt fast enough. Therefore, the aftermath of climate change will be far-reaching and profound .^14,15

Species distribution has always been one of the most important ecological aspects of biodiversity protection which has been affected by climate change through time. There is considerable support for the importance of environmental variables in limiting species distributions.¹⁶ Managers and conservationists widely use statistical species distribution models (SDMs) to predict the potential range shifts in species distributions under climate change scenarios.^17,18 SDMs can appropriately assess the impact of climate change on species distribution on both terrestrial [e.g., 19, 20] and marine species [e.g., 21]. A set of environmental conditions suitable for a given species has been defined as a bioclimatic envelope.²² It has been obtained by associating the current climate variables with species occurrences using statistically or theoretically derived methods. Assessing shifts of the bioclimatic envelope under climate change scenarios can predict any shift in the species current range.²¹

There is an important question that is how plant species will respond to climate change in the future, especially drought-resistant plant species. The drought-tolerant plant species grow naturally in different climates. Adopting drought-tolerant plants to conserve water is a potential adaptation to reduce the consequences of climate change .^23,24

Accordingly, it was hypothesized that climate change will not affect potential distributions of drought-tolerant species. Here, this hypothesis was tested using ten different algorithms to model the potential distribution of three drought-resistant plant species, namely Agropyron cristatum, Agropyron desertorum, and Festuca arundinacea in Iran under current and future climate conditions. Tall fescue (F. arundinacea) is widely distributed in temperate regions. It is a perennial bunchgrass with the ability to grow in a range of climates and has shown high tolerance to diseases and weeds.^25,26 These properties have made F. arundinacea popular in agriculture, landscaping, and phytoremediation. However, hot summer temperature adversely affects this species and reduces its growth and productivity.²⁷ Genus Agropyron includes the wild relatives of wheat. Species in the genus exhibit a host of useful traits such as tolerance to biotic and abiotic stressors such as low temperatures, salt, pathogens, and pests.^28,29 Agropyron cristatum and A. desertorum represent the typical traits of genus Agropyron and are found in the majority of temperate regions. All three species are used as forage for livestock and wildlife and have found applications in weed control, soil conservation, and watershed management. All three species are widely distributed in Iran. Agropyron desertorum, and F. arundinacea are distributed in the north, northeast, northwest and center and A. cristatum is mainly distributed in the northern regions of Iran.^30,31 Our goals include: (1) to determine the current and future distribution of the plant species in Iran, (2) to determine the major environmental parameters that dictate habitat suitability for the plant species, and (3) to evaluate the effect of climate change on habitat suitability for the plant species in the future.

Methods

Study area

Iran is located in the Middle East between latitudes 24° and 40°N, and longitudes 44° and 64°E at 1,648,195 square kilometers (Figure 1). Iran has a broad spectrum of climatic conditions across regions with significant annual precipitation variability (120 to 2000 mm yr⁻¹) and temperature variability (−20 to 50°C). Previous studies revealed that wet regions of the country would receive more precipitation while dry regions will receive less. Analysis of daily precipitation intensities indicated more frequent and larger‐intensity floods in the wet regions and more prolonged droughts in the dry regions [e.g., 32, 33].

Figure 1. Study area location and studied species occurrences

Presence records

Occurrence records and distribution of the species were obtained from herbariums of Ferdowsi University of Mashhad, Tehran University, and University of Birjand. Relevant literature^34–37 was also consulted, including Flora Iranica³¹ and Flora of Iran.³⁰ Herbaria data were obtained from field samplings between 2009 and 2019. The coordinates of all the occurrence points were recorded using a hand-held multichannel Global Positioning System (GPS) receiver with a positional accuracy of ±5 m. The spatially correlated presence records were removed using spatial autocorrelation and Moran’s I test. The number of presence records was 156 for Agropyron cristatum, 171 for Agropyron desertorum, and 154 for Festuca arundinacea (Figure 1).

Environmental variables

Topography, geology, soil and climate variables were supplied as input for SDMs. The National Cartographic Center of Iran (NCC) and Geological Survey and mineral exploration of Iran (GSI) provided topography, geology, and soil variables at 1:25000 scale. Mean elevation and slope were the only topography variables used in modeling. The digital elevation model was used to generate mean elevation and slope variables using ARCGIS software version 10.2 at 1 km² spatial resolution. The soil variable includes soil orders that consist of 20 classes. The classes include Alfisol, Aridisol, Entisol/Alfisol, Entisol/Inceptisol, Inceptisol, Inceptisol/Vertisol, Mollisol, Badland, Coastal Sand, Dune Land, Kalut, Marsh, Playa, Rock Outcrop/Entisol, Rock Outcrop/Inceptisol, Rocky Land, Salt Flat, Salt Pug, Urban, Water Body. The geology variable of the study area has 557 classes. The geologic map shows the types and ages of rocks and younger deposits found at or near the Earth surface. For climate variables, the set of 19 climate variables was provided by the Worldclim database³⁸ at 1 km² spatial resolution. Monthly solar radiation for a 12-month period was also obtained from Worldclim database at 1 km² spatial resolution. All variables were used as raster with 1 km² spatial resolution prepared by a resample analysis in ARCGIS software version 10.2. The climate data contains monthly air temperature and precipitation, which has been obtained from a global network using long-term time series from 1950 to 2000.

The climate projections for the year 2070 under 62 climate change scenarios, including four representative concentration pathways (RCPs) and 19 general circulation models (GCMs), was the base for the future climate data. Different scenarios were used to select the best-performing projection. These data were used to study habitat suitability for the selected species in the present and future. Representative concentration pathways were proposed by the IPCC in its fifth report in 2014 as predictions of greenhouse gas concentrations.³⁹ The four concentrations (490, 650, 850, and >1370 ppm CO2-equivalent) are titled RCP 26, 45, 60, and 85, respectively.⁴⁰ GCMs simulate the general circulation patterns of a planet atmosphere or oceans using mathematical models (see Table S4). Topography variables were excluded from the forecasting model of the future because future changes are expected to happen in the climate variables linked to topography and not topography by itself.²⁰

To omit redundant variables, the correlation between variables was tested using Pearson’s correlation analysis [‘Raster Correlations and Summary Statistics’ implemented in SDM Toolbox: cf. 41]. If −0.7 > r > +0.7, one variable was omitted from modeling. Doing so produced a set of final 12 variables used as input for habitat modeling (Table 1).

Table 1. Environmental variables to model habitat suitability of the studied species

Environmental variables	Unit
Monthly solar radiation	kJ m-2 day-1
Precipitation of Driest Quarter	mm
Mean Diurnal Range	°C
Temperature Seasonality	-
Precipitation of Wettest Quarter	mm
Mean Temperature of Coldest Quarter	°C
Mean Temperature of Warmest Quarter	°C
Annual Precipitation	mm
Elevation	m
Slope	degree
Geology	577 classes
Soil order	20 classes

SDM analysis

Ten different algorithms in four broad categories (machine learning, classification, regression-based, and enveloping) were used to determine habitat suitability for the species (Table 2). The general framework for this study is shown in Figure S1.

Table 2. The SDM algorithms used in this research and their required type of dependent variables

Category	Algorithm	Type	Reference
Machine Learning	Artificial neural networks [ANN)	P/A	42
	Boosted regression trees [BRT]	P/A	43
	Maximum entropy [MaxEnt]	P/B	44
	Multivariate adaptive regression splines [MARS]	P/A	45
	Random forest [RF]	P/A	46
Classification	Classification and regression trees [CART]	P/A	47
	Flexible discriminant analysis [FDA]	P/A	48
Regression-based	Generalized additive models [GAM]	P/A	49
	Generalized linear models [GLM]	P/A	49
Enveloping	Surface range envelope [SRE]	P/B	50

P: Presence; A: Absence; B: Background

To develop the algorithms for the studied species, the Biomod2 package, version 3.3–7 in R version 3.1.25¹⁴ was used.⁵¹ Since both presence and absence of species are required to fit SDMs, and because only the presence records were available, pseudo-absences were generated randomly with an equal number of presences.⁵² In our study, bootstrapping with 10 cycles was performed. The following two indices of model performance were used to evaluate the models by randomly assigning 70% and 30% of presence records as training and test data, respectively. The first index is the area under the ROC–plot curve [AUC; 53] in which 0 and 1 mean the worst and best predictions, respectively. The second index is the True Skill Statistics [TSS; 54], in which −1 and 1 mean the worst and best predictions, respectively. The best model was selected by the highest TSS and AUC values among models and climate change scenarios. Here, model performance is either categorized as poor (< 0.4), moderate (0.4–0.7), and good (>0.7). The importance of each environmental variable for each species in the models was assessed in Biomod2 by randomizing each variable individually and then projecting the model with the randomized variable while keeping the other variables unchanged. The model predictions containing the randomized variable were further correlated with those of the original models. Finally, the importance of the variable was calculated as one minus the correlation; higher values indicate predictors that are more important for the model.⁵¹ Finally, the best-performing model with current data was used to project it into future scenarios. Then, the obtained habitat types (poor, moderate, and good) were compared between current and future conditions to assess the possible change in habitats. New future climate conditions, that is, those not observed for an area currently, were detected using mobility-oriented parity (MOP) and the multivariate environmental similarity surface (MESS) .⁵⁵ MOP identifies future environmental conditions not available in current climate conditions, while MESS assesses the level of similarity between climate scenarios. Results showed that to establish those areas of strict model extrapolation from those areas with environmental conditions present in Iran. MOP and MESS were performed in R with a sample of 10% of the pixels in the environmental variables.

Results

The area under the ROC–plot curve and true skill statistics values were between 0.50 and 0.94 for the models in current and future times .The BCC-CSM1-1, RCP45 scenario reached the highest TSS and AUC values for the three species in future times . MaxEnt had the highest TSS and AUC values among the models for A. cristatum and A. desertorum, and random forest for F. arundinacea in the current and future times.

The importance of environmental variables varied for different species (Figure 2) and demonstrated different patterns (Figure 2). Environmental variables that contributed the most to modeling the potential distribution of all species were annual precipitation, the mean temperature of the warmest quarter, mean temperature of coldest quarter, and geology. Meanwhile, some variables have a high contribution to one species only. For example, solar radiation is important for modeling habitat for F. arundinacea, while precipitation is important for A. cristatum.

Figure 2. Importance of environmental variables for A. desertorum, A. cristatum, and F. arundinacea. (error bars indicate differences between current and future condition for each species, the importance of environmental variables generated by MaxEnt model for A. cristatum and A. desertorum and RF model for F. arundinacea under current and future condition (BCC-CSM1-1, RCP45 scenario))

Response curves revealed the direction of effect for the most important variables (Figure 3). The increased mean temperature of the coldest quarter and mean temperature of the warmest quarter caused a decrease in habitat suitability for the three species. Also, increased annual precipitation caused an increase in habitat suitability for A. desertorum and F. arundinacea, while the same variable had the opposite effect on A. cristatum. Comparing predictive maps with elevation classes indicated that all species would disappear in sites located below 500 m, and expand its range toward sites located between 500 and 3500 m during future projections (Figure 4).

Figure 3. Relationships between environmental variables and the probability of presence for A. desertorum, A. cristatum, using MaxEnt model and F. arundinacea using RF model under current and future condition (BCC-CSM1-1, RCP45 scenario) (red curves show the mean response and blue margins are ± SD)

Figure 4. Habitat suitability shifts of A. desertorum, A. cristatum, and F. arundinacea in relation to elevation under current and future condition. (the data generated by MaxEnt model for A. cristatum and A. desertorum and RF model for F. arundinacea under current and future condition (BCC-CSM1- 1, RCP45 scenario))

Our findings showed that suitable habitats (moderate and good habitat suitability) for all three species would decline in the future (Figure 5–7). Results also showed that 26% of F. arundinacea’s moderate habitats and 23% of its good habitats would remain stable by 2070, while 66% of the moderate habitats and 46% of the good habitats for the species will be lost . Also, the moderate and good habitats for A. cristatum, respectively, 46% and 35% will maintain their current suitability by 2070, while 41% of the moderate habitats and 5% of the good habitats will no longer be able to support the species (Table 3). Of the moderate habitats for A. desertorum, 26% will retain their current level of suitability by 2070, and 30% of the good habitats for the species will still be able to support the species by 2070. However, 64% of the moderate habitats and 16% of the good habitats will be lost in the next 50 years (Table 3). The results of the multivariate environmental similarity surface analyses showed broad areas with similar environmental conditions between current and future climate conditions in all 62 climate change scenarios. The mobility-oriented parity analyses also revealed that the future projection at BCC-CSM1-1, RCP45 has the minimum extrapolated areas among the 62 climate change scenarios.

Figure 5. Habitat suitability map for F. arundinacea in the current and future time using RF model under current and future condition (BCC-CSM1-1, RCP45 scenario)

Figure 6. Habitat suitability map for A. cristatum in the current and future time using MaxEnt model under current and future condition (BCC-CSM1-1, RCP45 scenario)

Figure 7. Habitat suitability map for A. desertorum in the current and future time using MaxEnt model under current and future condition (BCC-CSM1-1, RCP45 scenario)

Table 3. Changes in area and proportion of suitable habitats for A. desertorum, A. cristatum, using MaxEnt model and F. arundinacea using RF model under current and future condition (BCC-CSM1-1, RCP45 scenario] (km² and %)

			Current time
	Species	Habitat	Poor	Moderate	Good
Future time		Poor	794559.8 (96.2%)	282072.4 (66.6%)	163634.2 (46.5%)
	F. arundinacea	Moderate	27375.4 (3.3%)	111683.7 (26.4%)	106671.2 (30.3%)
		Good	4083.2 (0.3%)	29141.9 (6.8%)	81335.2 (23.1%)
		Poor	1173947.1 (97.2%)	99939.7 (41.7%)	9194.7 (5.9%)
	A. cristatum	Moderate	29287.0 (2.4%)	112382.4 (46.9%)	89569.7 (58.1%)
		Good	3981.5 (0.5%)	26947.8 (11.2%)	55306.8 (35.8%)
		Poor	713654.2 (98.1%)	312546.1 (64.8%)	65893.1 (16.8%)
	A. desertorum	Moderate	12734.1 (1.7%)	126485.2 (26.2%)	204239.5 (52.2%)
		Good	1036.8 (0.1%)	43056.5 (8.9%)	120911.4 (30.9%)

Discussion

Based on our projections, all species will retreat to the central parts of their current suitable habitats (moderate and good habitat suitability) under climate change. The extent and spatial pattern of distribution shifts varied substantially, according to species (Figure 5–7). The species also respond similarly to variables in the current and future times. However, the optimal intervals might be different. For A. cristatum, under climate change scenarios, suitable habitats shift from the mean temperature of the warmest quarter of 10–18°C in current conditions to 20–25°C in the future conditions. Plant species could circumvent the impacts of climate change through migration. However, migration is limited by the displacement capability of species. Some species show the geographical pattern of adaptation by migrating to a higher elevation.^56,57 In this study, all species showed a range shift toward higher elevation (Figure 4). Such phenomenon was reported for mountain plant species in northern Africa (e.g., Rosa arabica), ^58,59 Middle East² and in the Mediterranean mountains.⁶⁰ Our study indicated that habitats in areas lower than 500 m are unsuitable for these plant species. In the near future, it would be expecting a fluctuation in the species occurrence between 500 and 3500 m.

Adaptation and acclimation can be interpreted indirectly by SDMs results under climate change.^61–63 Moreover, as a general principle, high genetic variation in natural populations enhances the potential for disease and pest resistance and adaptation to environmental change, including climate change.⁶⁴ Previous studies have shown that these three species have high levels of genetic variation in Iran.^65–68 The rate of genetic variation determines the ability of species to sense changes in the habitat and make a plastic response.⁶⁹ Therefore, SDMs should also consider genetic variation in plant species, which could improve accuracy .⁷⁰

The different responses depend not only on the climate variables that SDMs used but also on many non-climate factors as well as the individual plant morphological and physiological advancement and its ability to cope with an adverse climate.⁷¹ The three species showed different responses to predicted climate change due to species-environment interactions. Festuca arundinacea and Agropyron cristatum will, respectively, experience the most and least severe decline in suitable habitats in the next 50 years. This result is because decreased annual precipitation caused an increase in habitat suitability for A. cristatum, while the same variable had the opposite effect for A. desertorum and F. arundinacea. On the other hand, F. arundinacea grows on moist soils that decreased annual precipitation caused a decrease in habitat suitability. Agropyron desertorum, and F. arundinacea are predicted to lose more than 50% of their current distribution ranges (Table 3). The lost suitable habitats are related to 15, 18, 16% of the actual distribution of A. cristatum, A. desertorum, and F. arundinacea, respectively. Therefore, future climate change can significantly impact the shifting distributions of these species. Moderate and good habitats for these species are mostly located in rangelands that are particularly at risk of overgrazing. Overgrazing can substantially worsen the negative impact of climate change, especially in arid habitats, where fragile assemblages of plant and threatened animal species are present, representing a serious threat to maintaining the food chain’s function.⁷² Iran has the maximum density of domestic sheep among arid countries in the world, ahead of Sudan and Australia. The sheep and goat population in Iran has been doubled in the last 40 years.⁷³ Iran has twice as many domesticated animals as the capacity of its rangelands allows. In some areas, the livestock population has reached five times the current carrying capacity.⁷²

SDMs have limitations and assumptions. For example, SDMs can be influenced by sampling bias, ⁷⁴ sample size, ⁷⁵ multicollinearity, and spatial autocorrelation.⁷⁶ Also, SDMs can be used to forecast species distribution under future conditions but the level of uncertainty of those models can be very high.^77,78 Using MaxEnt in habitat evaluation studies has become a common practice.⁷⁹ Nevertheless, MaxEnt has drawbacks such as presence-only models underperform in species distribution modeling compared to presence-absence models, since absence data can help indicate absence points for species⁸⁰ and MaxEnt can take advantage of museum and atlas data, which provide presence-only data. However, there are biases in data provided by museums, herbariums, and atlases toward accessible areas or areas where species are expected to occur. Our data face the same challenges. Rather than being collected with prior sampling strategies such as random or systematic sampling, data in museums and atlases are gathered without particular sampling schemes, and without any particular purpose in mind.⁸¹ Biases in sampling intensity can arise because data has not been collected according to a prior plan (such as random or stratified sampling).⁸² As opposed to presence-absence models, presence-only models such as MaxEnt are strongly affected by such sampling biases. On the other hands, MaxEnt has also been shown to be superior to other species distribution models, specifically for evaluating expanded uncertainty in modeling.

Our results have demonstrated that even plant species drought-stress tolerant are not immune to climate change and will see their current distributions undergo significant changes due to climatic changes. It should be noted that this research used the BCC-CSM1-1, RCP45 scenario because the future projection in BCC-CSM1-1, RCP45 has the minimum extrapolation areas among the 62 climate change scenarios. However, future economic development, especially in China and the United States, will determine the direction of development in the world. As a result, an increase in greenhouse gas emissions and an increase in the Earth’s temperature are expected. These future phenomena will possibly cause further impacts on plant distributions. These results are crucial for notifying long-term management activations aimed to reduce the effect of wildlife-human conflict.

Data Availability

Some or all data, models, or code that support the findings of this study are available from the corresponding author upon reasonable request.

Acknowledgments

This work was supported by Iran National Science Foundation [grant numbers 96015179].

Supplementary material

Supplemental data for this article can be accessed on the publisher’s website.