Habitat quality modelling and effect of climate change on the distribution of Centaurea pabotii in Iran

Climate change resulting from increased greenhouse gases affects the distribution of weeds by commonly expanding and shifting their future distribution. In this study, habitat distribution of Behbahanian Knapweed (Centaurea pabotii) was modelled as an endemic weed of wheat fields in four provinces in the southwest of Iran. Then, the current and the predicted future distributions were compared under two scenarios based on the lowest and highest carbon dioxide emissions. Field survey was carried out during March-May of 2015-2017 to collect presence points of C. pabotii. Habitat modelling was done using MaxEnt software using eight environmental variables and 25 presence points. To predict the future distribution, modelling projection of CCSM4 was performed for the year 2070 under scenarios of representative concentration pathways (RCP) 2.6 and RCP 8.5 using the current and the projected future bioclimatic variables in MaxEnt. Our results revealed that the suitable area of distribution will be approximately doubled in the future for both scenarios and will be shifted to lower latitudes and higher altitudes. Also, in the most western province of the study area, a new isolated and large suitable area will occur in the future. Therefore, it was suspected that this plant will be expanded to the wheat fields of this province. Expanding and shifting in the distribution of C. pabotii should be taken into consideration by agricultural managers in Iran. Additional keywords: greenhouse gases; MaxEnt; weeds; wheat fields. Abbreviations used: AUC (Area Under the Curve); BP (Bushehr Province); DEM (Digital Elevation Model); FP (Fars Province); KBP (Kohgiluyeh and Boyer-Ahmad Province); KP (Khuzestan Province); ROC (Receiver Operating Characteristic); SDM (Species Distribution Model). Authors’ contributions: KA: habitat and corridor modelling. AZ: drafting of the manuscript. KN: field survey and plant identification. KDH: field survey. All authors read and approved the final manuscript. Citation: Almasieh, K.; Zoratipour, A.; Negaresh, K.; Delfan-Hasanzadeh, K. (2018). Habitat quality modelling and effect of climate change on the distribution of Centaurea pabotii in Iran. Spanish Journal of Agricultural Research, Volume 16, Issue 3, e0304. https://doi.org/10.5424/sjar/2018163-13098 Received: 25 Feb 2018. Accepted: 03 Oct 2018. Copyright © 2018 INIA. This is an open access article distributed under the terms of the Creative Commons Attribution 4.0 International (CC-by 4.0) License. Funding: Research and Technology of Agricultural Sciences and Natural Resources University of Khuzestan (951/29). Competing interests: The authors have declared that no competing interests exist. Correspondence should be addressed to Kamran Almasieh: almasieh@ramin.ac.ir; Kamran.almasieh@gmail.com


Introduction
Over the last two decades, climate change has had a tremendous impact on natural and human-made ecosystems (IPCC, 1995).Climate change caused by human activities will increase the concentration of greenhouses gases like carbon dioxide (CO 2 ) and methane (CH 4 ) in the future (Patterson, 1995;IPCC, 2014).Increasing CO 2 , as the most important greenhouse gas, affects plants by changing the species distribution, alteration in reproduction timing and length of growing seasons (Peters et al., 2014;Remya et al., 2015).Commonly, the distribution of plants will be expanded to higher latitudes and elevations (Hijmans & Graham, 2006).
Climate change directly influences arable weeds through raising the temperature and changing the precipitation pattern (Peters et al., 2014).Weed plants generally expand their geographic distribution as a result of climate change; however, in some cases, there will be a contraction in their distribution (Wilson et al., 2009;Ramesh et al., 2017).Distribution shifting has been recorded for European arable weed species in recent decades by many researchers (e.g., Andreasen & Streibig, 2011;Kolarova et al., 2013;Salonen et al., 2013;Rasmussen et al., 2017).Moreover, in the southern hemisphere, for instance in South Australia, most species are expected to shift to southern parts of their current distribution in the future (Scott et al., 2014).
Species distribution models (SDMs) are widely used in habitat modelling of species (Fitzpatrick et al., 2013;Hu et al., 2015, Almasieh et al., 2016).There are several methods used for habitat modelling (e.g., GLM, GAM, BIOCLIM, GARP, MaxEnt and ENFA); some, like GLM, require both presence and absence points and others, like MaxEnt, require only presence points (Phillips et al., 2006).Recently, SDMs were used to predict the effects of climate change on plant species (Thuiller et al., 2005).These models assess the relationship between species presence and associated climatic variables to predict the current and the future distribution of species (Elith et al., 2006;Remya et al., 2015).The MaxEnt model has proved effective in prediction of future distribution of weeds (Wilson et al., 2009).
Centaurea sp.L., with ca.104 species, is among the largest genera of the Asteraceae family in Iran (Wagenitz, 1980;Negaresh & Rahiminejad, 2015a, 2016).Behbahanian knapweed (Centaurea pabotii Wagenitz) as a Zagrosian species is an endemic species to the south and southwest of Iran, which is found near roads and farms, especially wheat fields (Negaresh & Rahiminejad, 2015b).Due to the high potential dispersal of this arable weed (Negaresh & Rahiminejad, 2015b), we assumed that the distribution of C. pabotii might increase in the future as a result of climate change.Therefore, in this study, our aim was to model the habitat distribution of this plant in Iran and compare the current distribution with the future distribution under two scenarios based on the lowest and highest carbon dioxide emissions.Understanding the effects of climate change on the distribution of this weed plant of wheat fields will be useful for agricultural managers in Iran.

Study area
The study area included Khuzestan Province (KP), Bushehr Province (BP), Fars Province (FP), and Kohgiluyeh and Boyer-Ahmad Province (KBP) with a total area of about 223000 km 2 in the southwest of Iran (Fig. 1).In this area, there were two relatively distinct topographic situations: (i) vast arid plains in the south of KP and west of BP with high temperatures in summer averaging ~37°C and mild winters averaging ~17°C (IRIMO, 2010), and (ii) mountainous areas in the north of KP, east of BP as well as the other two provinces (FP and KBP) with moderate temperatures in summer averaging about 25°C and cold winters averaging about 6°C (IRIMO, 2010).

Data collection
Field survey was done in the study area during March-May of 2015-2017 and presence points of C. pabotii were obtained.This species was classified by the Flora Iranica (Wagenitz, 1980) as well as examination of relevant taxonomic literature (Negaresh  & Rahiminejad, 2015b& Rahiminejad, , 2016)).Presence points were recorded by a GPS with an error of less than 10 m.The spring season was chosen to conduct field surveys because growing and flowering of C. pabotii in this season facilitated the identification of the species.Altogether, we obtained 25 presence points (Fig. 1).

Habitat modelling
Habitat quality map was obtained to show areas with high quality for C. pabotii using MaxEnt v3.4.1 (Phillips et al., 2017).MaxEnt only needs a few presence points to create pseudo-absence points (here we selected 10,000 points), and compares environmental variables in presence cells with values of the same variables in pseudo-absence points to find the discrimination between them (Phillips et al., 2006;Almasieh et al., 2016).Finally, the software creates a continuous quality map.The area under the receiver operating characteristic (ROC) curve (AUC) in MaxEnt (0.5-1, in which ˈ0.5ˈ indicates random and ˈ1ˈ indicates complete discrimination of presence points from pseudo-absence points) was used to evaluate model performance for both habitat quality modelling and climate change predictions.Jackknife test within MaxEnt was used to evaluate the relative importance of each variable and response curves within MaxEnt were used to show the probability of species presence for each variable considering the correlation between each variable and other variables.We randomly selected 75% of the presence points as the training data set and the other 25% as test data (Pearson et al., 2007).
Environmental variables including topographic, climatic and land-cover were used for habitat modelling (Table 1).Digital Elevation Model (DEM) at about 1-km spatial resolution was used to create the slope variable using Spatial Analyst tool.Solar radiation variable (kJ m -2 day -1 ), which was classified for each month of the year (averaged period of 1970-2000) (Fick & Hijmans, 2017), was used for modelling.We averaged solar radiations of April-May using Raster Calculator tool, because this quarter of the year was necessary for growth of C. pabotii.Regarding C. pabotii's dependence on wheat fields, we extracted agricultural lands, the majority of which were wheat fields, from the land-cover variable (FRWMO, 2010) and distance from agricultural lands was created using Euclidean Distance tool.All the tools mentioned are available in ArcGIS 10.2.
Among 19 climate variables (Hijmans et al., 2005), we chose annual mean temperature (BIO1), temperature annual range (BIO7), annual precipitation (BIO12), precipitation of the wettest month (BIO13) and precipitation of the wettest quarter (BIO16).We included those climate variables which we assumed had the greatest physiological impact on C. pabotii (Bradley, 2014); Annual temperature and precipitation (BIO1 and BIO12), maximum temperature of the warmest month minus minimum temperature of the coldest month (BIO7), as well as precipitation of the wettest periods (BIO13 and BIO16) which were important for growth of C. pabotii were chosen.We removed BIO1 due to its high correlation with DEM.All variables had a correlation of less than 0.85 (Yang et al., 2013).
We categorized continuous habitat quality map based on maximum training sensitivity and specificity logistic threshold in MaxEnt (Jimenez-Valverde & Lobo, 2007).Therefore, 0.248 was considered as a threshold of suitable and non-suitable categories.We also categorized the suitable category into three subcategories: 0.246-0.4as low quality, 0.4-0.6 as moderate quality, and 0.6-1 as high quality (Hu et al., 2015).

Climate change prediction
We used MaxEnt to predict the effect of climate change on the distribution of C. pabotii in the future.MaxEnt can predict future distribution of species more reliable than other prediction models (Hijmans & Graham, 2006).We performed modelling projection in MaxEnt using current bioclimatic variables (i.e., BIO1, BIO7, BIO12, BIO13 and BIO16) as well as future

Habitat modelling
The AUC of ROC curves were 0.984 and 0.978 for training and test data respectively, which indicated the high performance of the model.The high quality category was mostly located in KP and also several small ranges of the high quality category were located in KBP and BP.FP was the province with the smallest area of high quality category (Fig. 2).
BIO13 variable was the most important variable that solely provided the most useful information (Fig. 3A).C. pabotii preferred precipitation of more than 100 mm in the wettest month (Fig. 3B).Also, DEM, BIO16 and BIO12 variables were the second, third and fourth most important variables, respectively (Fig. 3A).C. pabotii inhabited areas with an elevation of less than 1000 m and preferred precipitation of more than 200 mm in the wettest quarter and 300 mm as annual precipitation (Fig. 3B).Distance from agricultural lands was the variable that, if omitted, the highest reduction in gain of the model would appear.Therefore, distance from agricultural lands provided the greatest volumes of information, which were not present in other variables (Fig. 3A).Due to the dependence of C. pabotii to agricultural lands, the probability of the presence of this species decreased as distance from agricultural lands increased (Fig. 3B).bioclimatic variables.We downloaded the projection model of CCSM4 for the year 2070 (average 2060-2080), which was from the fifth phase of coupled model inter-comparison project (CMIP5).Four representative concentration pathways (RCPs) were available based on prediction of greenhouse gases concentration (i.e., RCPs 2.6, 4.5, 6.0 and 8.5).Each pathway explained the radiative forcing value as the amount of energy entering the atmosphere and the quantity of reflex.By the end of the year 2100, the radiative forcing levels are estimated to reach 2.6, 4.5, 6.0 and 8.5 W/m 2 of surface (Van Vuuren et al., 2011;Yousefi et al., 2015).We chose RCPs 2.6 and 8.5 as the lowest and highest radiative forcing levels for this study.
We assessed changes in C. pabotii distribution due to climate change for the future under both scenarios of RCPs 2.5 and 8.5.For that, we categorized habitat quality map of the current and future distribution into two categories (i.e., non-suitable and suitable) based on maximum training sensitivity and specificity logistic threshold in MaxEnt.Then, we calculated the suitable area for the present and the future, as well as the changes in the suitable area in the future under both scenarios.Also, to assess the trend of shifts in the suitable area, we calculated the centroids of merged suitable polygons for present and future maps using Hawths tools (Beyer, 2004) and finally, we compared the spatial change of centroids of suitable polygons (Hu et al., 2015).

Climate change
The AUC of ROC curves were 0.963 and 0.884 for training and test data of scenario of RCP 2.6, and 0.934 and 0.928 for training and test data of scenario of RCP 8.5 respectively, which showed the high performance of model.Our results revealed that the suitable area of the C. pabotii distribution will increase in the future in all provinces.Nevertheless, KP showed the highest increase in suitable area in the future under both scenarios (Table 2 & Fig. 4).We observed a shift for C. pabotii toward lower latitudes in the future based on centroids of polygons, as well as a shift toward higher elevations (averaging 770 m to 1406 m in  both scenarios).The suitable area in the future will approximately double under the two scenarios and will expand to the surrounding areas, especially toward the center of KP and north of BP.Also, a new relatively large and isolated suitable polygon was observed in FP under both scenarios (Fig. 4).We observed more expansion of the suitable area in RCP 8.5 (Fig. 4B) with respect to RCP 2.6 (Fig. 4A), particularly more expansion toward the center of KP.

Discussion
Our habitat quality modelling revealed that KP had the largest high quality category, which was in accordance with the previous findings about the distribution of the species (Negaresh & Rahiminejad, 2015b).There were also several small patches of modelled high quality categories in other provinces.
BIO13 variable had the highest contribution to habitat modelling in MaxEnt.It seems that precipitation in the wettest month was necessary for the growth season of C. pabotii.DEM was the second most important variable in habitat modelling.Negaresh & Rahiminejad (2015b) reported a range of 30 to 450 m for the distribution of this species; whereas, our response curves revealed that there was no C. pabotii presence observed in elevations of more than 1000 m.Therefore, Zagros Mountains with elevations of more than 1000 m are not appropriate habitats for this species.This plant is dependent on wheat fields and in the absence of the variable of distance from agricultural lands, our MaxEnt model had the highest reduction in gain.
Based on the current and the future model, KP still will have the highest distribution of C. pabotii in the future.Therefore, KP wheat fields seem to be mostly affected by this arable weed in the future.In this study, we considered four provinces in the southwest of Iran as the study area.We collected presence points in the most western border of FP.However, we considered the whole area of FP because we believe that as a result of climate change, C. pabotii distribution will expand into this province.Our results showed that the suitable area in FP will be approximately seven times larger in the future based on both scenarios (Table 2).Centaurea pabotii's distribution will expand to FP through a continuous expansion into four provinces, as well as an isolated large patch in FP in the future (Fig. 4).
Based on the current and the future model, C. pabotii will shift to higher elevations (with a mean of about 1400 m) that revealed that more mountainous areas of Zagros Mountains could be inhabited by C. pabotii in the future.Unlike other studies that reported a shift of species distribution toward higher latitudes in both plants and animals (e.g., Hu et al., 2015;Yousefi et al., 2015), we observed a distribution shift toward lower latitudes, because KP has vast plains adjacent to the Persian Gulf at lower elevations and higher latitudes compared to the other two provinces of the study area.Zagros Mountains span these provinces (FP and KBP) in the southwest and extend toward lower latitudes in Iran (Fig. 1).Therefore, it is anticipated that C. pabotii will be distributed to lower latitudes.
In summary, results of this study proved our hypothesis about the future expansion of C. pabotii distribution as an arable weed of wheat fields in Iran.This plant will tend to inhabit higher elevations and lower latitudes due to the temperature rise in the future.On the other hand, it seems that crops, such as wheat, will reproduce even more in the future resulting in higher CO 2 (Ziska et al., 2005;Tokatlidis, 2013;Peters et al., 2014) and the competition of wheat and C. pabotii in the southwest of Iran will be unknown in the future.In this respect, strong consideration should be taken into account by agricultural managers in this area of Iran.

Figure 1 .
Figure 1.Study area based on the distribution of Centaurea pabotii and location of presence points in four provinces in the southwest of Iran.

Figure 2 .
Figure 2. Categorized modelled habitat quality of C. pabotii in the study area based on 25 presence points and eight environmental variables in MaxEnt model.

Figure 3 .
Figure 3. (A): Jackknife test of the importance of environmental variables in MaxEnt model; when only the BIO13 variable was considered, the model produced the highest gain compared to when other variables were used separately.Excluding distance from agricultural lands, the highest reduction in gain of the model would appear.(B): Response curves of the C. pabotii presence to environmental variables within MaxEnt model, for each variable with respect to its correlation with other variables.

Figure 4 .
Figure 4. Effect of climate change on the potential distribution of C. pabotii in the study area based on CCSM4 projection model under two scenarios: RCPs 2.6 (A) and 8.5 (B).

Table 1 .
Environmental variables used for habitat modeling of Centaurea pabotii in Iran.

Table 2 .
Potential suitable area (in parentheses % of the total study area) and changes in potential suitable area of Centaurea pabotii in year 2070 due to climate change based on projection model of CCSM4 under two scenarios of RCPs 2.6 and 8.5 in MaxEnt model.