Modelling the combined effects of land use and climatic changes: coupling bioclimatic modelling with markov-chain cellular automata in a case study in Cyprus

To appear in Ecological Informatics in 2015 (acce pted manuscript ver sion prior to proofs) Modelling the com bined effects of land use and c limatic c hanges: coupling bioclimatic modelling with markov-chain cellular automata in a case study in Cyprus Marianna Louca 1 , Ioannis N. Vogiatzakis 1 , and Aristides Moustakas 2, * 1. School of Pure & Applied Sciences, Open University of C yprus, Nicosia, Cyprus 2. School of Biological a nd Chemical Sciences, Queen Mar y University o f L ondon , Mile End Road, E1 4NS, London, UK * Corresponding author: Aristides (Aris) Moustakas Email: arismoustakas@gmail.com Abstract Environmental change in terms of land use and climatic change s are posing a serious threat to species distributions and extinctions. Thus in order to predict and mi tigate against the effects of environmental ch ange both drivers should be accounted for. Two endemic plant species in the Mediterranean island of Cyprus, Crocus cyprius and Ophrys kotschyi, were used as a c ase study. We have coupled climate change scenarios, and land use change models with species distribution models. Futu re land use sc enarios were modelled by initiall y calculating the rate of current land use changes between two time snapshots (2000 and 2006) on the island, and based on these transition probabilities markov -chain cellular automata were used to generate future land use chan ges for 2050. Climate change scenarios A1B, A2, B1 and B2A were derived from the I PCC reports. Species’ climatic preferences w ere derived from their current distributions using classification trees while habitats preferences were derived from the R ed Data Book of the Flora of C y prus. A bioclimatic model for Crocus cypri us was built usin g mean temperature of wettest quarter, max temperature of w armest month and precipitation seasonality, while for Ophrys kotchyi the bioclimatic model was built using precipitation of wettest month, mean temperature of warmest quarter, isothermalit y , pre cipitation of coldest quarter, and annual precipitation. Sequentially, simulation scenarios were performed regarding future species dist ributions by accounting climate alone and bot h climate and land use changes. The dist ribution of the two species resulting from the biocl imatic models was then filtere d by future land use changes, providing the species’ proje cted potential distribution. T he species’ projected potential distribution varies depending on the t ype and scena rio used, but many of both species’ current sites/locations are projected to be outsi de the ir future potential distribution . Our results demonstrate the importance of includin g both land use and climatic changes in predictive species modeling. Keywords Spatio-temporal computation; Species distribution modelling; land use changes; climatic changes; habitat modification; model coupling; endemic species; Crocus cyprius ; Ophrys kotschyi Introduction Climate and land use changes are two of the main causes for biodiversity loss worldwide, while their combined effects ma y be greater than either of t hese fa ctors acting alone (de Chazal and Rounsevell, 2009). Climate change has already affected species distribution, abundance, phenology and interactions, while even greater impacts are ex pected in the future (Rosenzweig et al., 2007 ) with thre e major options for threatened spe cies: (a) extinction, (b) evolution and subsequent adaptation and (c) shifting its geographic ran ge to more favourable conditions (Moustakas and Evans, 2013). However, t he rate at which climate change is happening today is often faster than the ability of some species to disperse or ad apt, while other factors such as land use change and habitat fragmentation impede thei r ability to move to suitable areas (Thuiller et al., 2005). Climate change models simulate the change in climate due to the accumulation of greenhouse gases, based on the current underst anding of atmospheric phys ics and chemistr y (Hannah, 2010). The Intergovernmental Panel on Climate Change (IPCC) has produced a range of emission scenarios for use in global climate mode ls that predict future climate (IPCC, 2013 ). The l atest I PCC report is based on alternative concentrations of greenhouse gases without being associated with any socio -economic scenario, but instead could result from different combinations of economic, technological, demographic, policy, and institutional futures (IPCC, 2013 ). This change facilitates better integration of so cio- economic factors, such as land use changes into future projections. Species Distribution Models (SDMs) use inform ation on the locations of species and their correspondin g envir onmental covariates, creating statistical functions to be projected in areas or time periods wh ere environmental parameters are known but species dist ribution is unknown, providing inference for potentially suitable sites (Brotons et al., 2004). In addition to climate change, the destruction, fragmentation and degradation of habitats due to changes in land use are among th e strongest pressures on biodiversity (EEA, 2010). In analogy to climate modelling, land use change models use a variety o f approaches to assess and project the future role of land use change on biodiversit y , soil degradation, the ability of biological systems to support human needs and the vuln erability o f places and people to climatic, eco nomic, or sociopolit ical perturbations (Zhang et al., 2014 ). The development of land use change scenarios allows their integration into SDMs along side dynamic climatic v ariables which can significantly improve a model’s explanatory and predictive abilit y at fi ne scales (Martin et al., 2013). Land use can be incorporated int o th e model as static variables that do not change over time (Iverson and Prasad, 2002), or as dynamic variables that change under different scenarios ( Sc hweiger et al., 2012). Despite the fa ct the combined effects of climate and land us e change affect species distributions (Martin et al., 2013), only a small number of SDMs predict species dist ribution based on both of these factors (e.g. Esteve-Selma, et a l., 2012; Sch weiger, et al., 2012 ; Heubes, et al., 2013), while most SDMs that combine climatic variables with land use variables, onl y use dynamic variables for climate, while land use is consider ed stable (M artin et al., 2013 ). D y namic model coupling (Verdin et al., 2014 ) of climate models and land use models can b e emplo yed to account for the int erac tion between both effec ts in projected future conditions (Evans et al., 2013a). Cyprus is a biodiversit y hotspot (M y ers et al., 2000) wh ich is expected to become warmer and dr yer (Hadji nicolaou et al., 2011 ) . At the same time, the incr eased pressure for urban and tourism development in Cy prus is leading to sig nificant cha nges in land uses (Eurostat, 2012). Thus, the island of Cyprus is an ideal stud y area because of the major climate and land us e ch anges expected in the near future, the pres ence of a multitude of threatened endemic spec ies and the absence of simil ar studies in the region to date . We sought to quantify the c ombined effects of c limatic and land use c hanges on two plant endemic species. Materials and Methods Study area Cyprus is the thi rd largest Mediterranean island, with an area of 9251 Km 2 . The climate is Mediterranean, with hot and dry summers from June to September (little or no rainfall, average maximum tempe ratures up to 36°C) , rainy but mild winters from November to March and two short transitional seasons, a utumn and spring. For detailed information regarding the geomorphology and biogeography of the island see Supp. 1A. Target species The target spe cies a re C rocus cyprius Boiss. & Kotschy and Ophrys kotschyi H. Fleischm. & Sofi, both endemic to Cy prus, and categorized as vulnerable under the IUCN classification (Tsintides et al., 2007). The criteria considered to target these two species were : (i) high risk of extinction (ii) endemism (iii) high number of data occurrences/location relative to other available species (iv) si gnificant differences between their distribution s, as Crocus c yprius only occurs in the Troodos Mountains, while Ophrys kotschyi oc curs almost ever y where in Cyprus ex cept the Troodos Mountains. For additional information regarding the target species see Table S1 in Supp. 1A. Data Species distribut ion data were obtained from the Red Data Book of the Fl ora of C y prus (Tsintides et al., 2007), in the form of true presence point s. The data were collected durin g a systematic extensive survey between 2002 and 2006 (Tsintides at al., 2007). There were 102 true presences for C. cyprius and 117 for O. kotschyi .. From these, a set of “theoretical” presences were deriv ed for each species ; these comprised centres of cells of the potential habitat, i.e. the entire area with a suitable al titude, soil and land use, within which at least one tru e presence had been r ecorded. Absence data were created artifici ally, from background data of the entire potential habitat of each species; “ these comprised the centres of cells of the potential habitat where no true presence had b een recorded. This w as done in order to provide a s ample of the set of conditions available to the species in the region and not to pretend that the species is absent in the selected sites (Phi llips , et al., 2009). The data were then wei ghted to simulate prevalence 0.5, i.e. th e total weight of pr esence is equal to the total weight of absences (Barbet-Massin et al., 2012). Bioclimatic da ta were obtained from Worldclim database, Version 1.4 (release 3), which is available on www.worldclim.or g (Hijmans et al., 2005) . The bioclimatic variables used are shown in Table S 2 in Supp. 1B. Future bioclimatic data were also obtained from Worldclim for the year 2050, according to GC M HadCM3 (Hadle y Centre Coupled Model, Version 3) and A1B, A2, B1 and B2A SRES emission scenarios. Each sce nario is based on a different “storyline” and scenario famil y (A1, A2, B1 or B2), representing different demographic, social, ec onomic, technological, and environmental developments (IPCC, 2013). The four storylines combine two sets of di vergent tendencies: one set var y ing between strong economic valu es (A1 and A2 families) an d strong environmental values (B1 and B2 families), the other set between increasing globalization (A1 and B1 families) and increasing regionalization (A2 and B2 families) (Nakicenovic & Swart, 2000 ). . Land use data for 2000 and 2006 were obtained from CORINE database, available on http://www.eea.europa.eu/data-and-maps . The resolution for both the current and the future bioclimatic data was 0.71 Km 2 while for the land-use data 250 m x 250 m. Habitat preferences were determined based on the information provided in The R ed Data Book of the Flora of Cyprus, which was the result of a sy stematic study of all available information on the threatened plants of Cyprus, in combination with field work (Tsint ides et al., 2007). The combinati on of suitable altit ude, soil and land use, as described in Tsintides et al. (2007) was defined as the species current pot ential habitat when usin g current land use and as future potential habitat when using future land use in ArcGIS ( http://www.esri.com ). SDMs The SDMs were created using C lassification Trees (CT), a m achine lear ning method used to c reate predictiv e models (Fig ures 2 a nd 3). The method pre dicts the value of a dependent variable with a finite set of values , from the values of a set of independent variables (J i et al., 2013 ). The m ain ad vantage of this method is that it does not require a specific t ype of data or that they follow a specifi c statistical distribution. The evaluation of the predictive accuracy of the model was measured usin g Cohen's Kappa (Congalton, 1991 ) and the "area under the curve ” (AUC) of Receiver Operating C haracteristic plot (ROC plot) (Fielding and Bell, 1997). The Classification Tree is considered to represent each species’ “bioclimatic envelope” or “bioclimatic space”, which is defined as the climatic component of the fundamental ecological niche, or th e ‘climati c niche’ (Pearson and Dawson, 2003). We used SPSS version 20 for the statistical analy sis ( htt p://www- 01.ibm.com/software/analytics/spss/ ). Land use prediction Future land use w as predicted with an inte gration of Markov chain and Cellular Automata ( Figure 4). Markov chain is a technique that has b een widel y us ed to predict changes in vegetation and which predicts future changes based on the rate of previous changes (A rsanjani et al., 2011 ). The main disadvantage of the Markov chain is its lack of spatial dimension: it gives ac curate information on the transition probabili ties of each land use type to another, bu t provides no information on the spatial distribution of changes (Eastman, 2003). This problem is solved by combining Markov chain with Cellular Automata (Arsanjani et al., 2011). Cellular A utomata are digital entities that have th e ability to change their state based on the previous condition of themselves and their neighbours, based on a specific rule (Moustakas et al., 2006). In order to mod el land use changes we used two different time snapshots 2000 and 2006 of the Cyprus COR I NE land- use maps. For each land use type we calculated the transition probabilities to all the other land use types as well as the probabilit y of no land use changes ( retain current land use t y pe). This was implemented in a spatiall y explicit manner i.e. transition probabilities were calculated for each location, where locations were represented b y 1 km grid cells. Thus, the transition probabilit y depends not only on the cell’s previous state but also on the state and rate of change of the neighbouring cells. Climate and land use change integration We generated future lan d use scenarios for 2050 using CA_MARK OV in Idrisi Selva software ( http://www.clarklabs.org (Eastman, 2003) (F igure 4). We used ArcGIS version 8 (www.esri.com) to map the modelled outputs of the sp ecies’ future pot ential distributions based on bio climatic modelling (bioclimatic envelope and future clima te da ta). We then combined the projected d istribution for 2050 based on bioclimatic space with the current and future potential habitat, p roducing future potential dist ributions based on climate change onl y and on climate and land use changes respectivel y (Figure 5). This was done b y simply taking the cells in which projec ted distribution base d on bioclimatic spac e and potential habita t overlap. The methodology that was used is summarized in Figure 1. Results SDMs The potential distribution of Crocus cyprius can be predicted using only three bioclimatic variables as deduced from the classification tree: Mean Temp erature of Wettest Quarter , Max Temperature of Warmest Month a nd, and Precipitation Seasonalit y; (Fig. 2). Model accurac y was “good” to “substantial” as in dicated b y Cohen’s K appa and AUC values of 0.719 and 0.859 respectivel y (Landis and Koch, 1977 ; Swets, 1988). The potential distribution of Ophrys kotchyi ’s can be predicted using five bioclim atic variables as deduced from the classifi cation tree: Precipitation of Wett est Month, Mean Tempe rature of Warmest Quarter , Isothermality, Precipitation of Coldest Quarter, and Annual Precipitation; (Fig. 3). Model accurac y was also good (Kappa = 0.619, A UC = 0.809). SDM pe rformance measures only relate to current distributions. Land use prediction The map of future land use produced from CORINE 2000 and 2006 maps using CA_MARKOV indicated that the biggest increase in 2050 comp ared to 2006 is predicted for ar tificial surfaces ( Fig . 4 ). Non-irrigated arable land presents the greatest reduction, followed by sclerophyllous vegetation and natural grasslands (Fig. 4). The classes with the lowe st predicted change are broad-leaved forest and mixed forest (Fig. 4). Climate and land use change integration Potential distribution s of Crocus cyprius in present conditions, considering climate change onl y and considering climate and land us e change suggest that in scenario B2A its potential distribution increases, while in all othe r scenarios it decreases, with A1B b eing the most pessimistic (Fig. 5a; Table 1). Th e decreased potential distribution is limited at higher altitudes, while the suit able areas at lower altitudes disappe ar. S hifting is predicted in all scenarios, so that 34 to 68 out of the 102 species’ current sites are excluded from its future potential dist ribution (Fig. 5 a) . The inclusion of land use change does no t cause significant changes (Fig. 5a; Table 1). In the case of Ophrys kotschyi , potential dist ribution is projec ted to de cline in all scenarios, with B2A b eing the most optimistic and A2 the most pessimistic scenario (Fig. 5b , Table 2). The decrease is not li mited to a specifi c ar ea but is sp read throughout the current potential distribution. Furthermore, shifti ng and f ragmentation result in the exclusion of 111 to 116 out of the 117 spe cies’ current sites from its future potential dist ribution (Fig. 5b). The inclusion of land use change causes additional reduction, from 16.2% to 17%, resultin g in the eradication of all of the species’ current locations in scenario B1 and the persistence of one or two locations in all other scenarios (Fig. 5b, Table 2). Discussion According to our model all scenarios predict the disappearance of Crocus cyprius ’ main area of oc currence in the top of the Troodos Mountains, resulting in the ex tinction of many of its curre nt sites. The ex tinction in this region is caused b y t he incr ease of the maximum temperature o f the warmest month, which in all scenarios exceeds the upper limit of 25.68 ⁰ C determined by the classification tree. Th e increased projected potential distribution does not necessarily imply that presen ce sites are safe (Supp. 1C and Table S3). Model outputs for Ophrys kotschyi predict future l ocations of a species adapted to low precipitations and relativel y high temperatures, but not to arid conditions – see Supp. 1C for more details. These conditions are found in the hi lly areas and the central plain of Mesaori a, where this species occurs toda y. Therefore, the combination of decreased precipitation of wettest month, increased mean temperature of warmest quarter and decreased precipitation of coldest quarter in all scenarios is responsible for the reduction of the species’ potential distribution (Table S 4 in Supp. 1C). Scenario A2, which is the most pessimistic, pr edicted the greatest reduction in precipitation of wettest mon th, which is the most im portant predictor of the species distributi on according to the classification tree. The loss, shifting and fragmentation of the O. kotshyi’ s potential distribution are likel y to be d isastrous, since the vast majority of the species’ present sites are not part of its future potential distribution in any scenario. The inclusion of land use change pla y s a different role in each species. For C. cyprius it does not cause an y significant changes compared to considering climate change onl y, while for O. kotschyi it is very im portant, causing an additional reduction in the species’ potential distribution and in the number of p resent sites included in it. This difference is caused by the different requirements of each spe cies regarding land use. For C. cyprius, suitable land uses consist only of coniferous forests, which in Cyprus mainly occur in side Natura protected areas and state land, thus they are not expected to undergo significant cha nge s in the future. In contrast, suitable lan d uses for O. kotschyi consist of 10 different classes, the most important of which are n on-irrigated arable land, scleroph y llous vegetation and natural grasslands. These three classes are expected to suffer the greatest reduction by 2050. To our kno wledge, no similar studies for the t arget species or species belonging to th e same genera were found in the literature. However, the results of this study are consistent with other studi es conducted in the Mediterranea n region. (Esteve-Selma et al., 2012 ) predicted an in crease or decrea se of the pot ential distribution of the endemic Tetraclinis articulata in Southeast Spain, depending on the emissions scenario considered: The potential distribution increases in scenario B2 and severely reduced in scenario A2. This is in agreement with our results for C. cyprius, but not with those of O. kotschyi , whose potential distribution is predicted to decrease in all scenarios. Our results are also consistent with those of (Vennetier and Ripert, 2009), who predict the disappearance of most forest areas with high species richness in Southeast France by 2050, using biocl imatic modelling. The projected limitation of C . cyprius ’ potential distribution at higher altitudes agrees with the predictions of other studies (Bell et al., 2014 ; Ferrarini et al.). In general the coupling of climate change and land us e change resulted in more restricted distributions for O. kotschyi is in agreement with the predictions of other studies (de Chazal and Rounsevell, 2009). Limitations, uncertainty and future directions There are a number of limitations related with the methodology employed in the current study. These include spatial and temporal resolution mismatch between the bioclimatic and the species distribution data, the possibility of sampling bias on the species distribution data, the quality and accuracy of the W orldclim (Hijmans et al., 2005 ) an d COR I NE data (http://www.eea.europa.eu/data-and-maps). In ad dition SDMs are based on a number of unrealistic assumptions ( Evans et al., 2013a). These include considering that the model quantifies the realized ecological niche (Pelliss ier et al., 2013) and that the species are in equilibrium with the environment, ignoring biotic interactions (Matias et al., 2014) and evolutionary and phenotypic changes as an adaptation to climate and land use change (Moustakas and Evans, 2013 ), assuming full dispersion for both species (Rodríguez-Re y et al., 2013) and training the model onl y in the realized environment (Maher et al., 2014). Also, it is considered that the potential distribution is the geographic area that meets one or more components of the fundamental ecological niche, when in fact the real potential distribution includes the entire fundamental ecologica l niche. As a re sult, the actu al future potential distributions are likely to be under- or over- predicted (Jiménez-Valverde et al., 2008). Although widely employed as a measure for SDM performance, AUC has been also criticised (see Lobo et al. 2008). Other important sources of uncertainty were the creation of theoretical presences and absences (Barbet-Massin et al., 2012). In addition predictions are often scale specific while the interaction of species with th eir environment take s pl ace at a variet y o f scales (Bellamy et al., 2013). Finally the inclusion of b oth climatic and l and use var iables introduces a new source of uncertainty, through the different parameters and assumptions it brings (Conlisk et al., 2013) . However simple models do not lead to generality in their predictions and thus increasing complexity may yield more realistic predictions (Evans et al., 2014 ; Evans et al., 2013b). Here model coupling (Verdin et al., 2014) of bioclimatic modelling with the use of markov chain cellular automata (Arsanjani et a l., 2011) was emplo ye d. Alternatives ma y include using individual based models (Gonzalès et al., 2013 ; Zhang et al., 2014 ) to model both land use & climatic changes with a presence only or presence absence model output and compare model outputs (Tonini et al., 2014 ) in an identical grid siz e. Although there is an interaction effect between land use and climate changes, disentangling th ese is conceptually complicated (see Lehsten et al. 2015) and is unaccounted for in this study. The procedure developed herein can be used for any species where data is available and provides a valuable and transfera ble method fo r understanding potentia l shifts in species distributions. In addition the maps produced can guide a ctions on adaptation and mitigation measures to climate change for the species studi ed, providing a useful tool for polic y and decision makers. References Arsanjani, J.J ., Kainz, W., Mousivand, A.J., 2011. Tra cking dynamic land -use change using spatially ex plicit Markov Chain based on cellular automata: the case of Tehran. International Journal of Image and Data Fusion 2, 329-345. Barbet-Massin, M., Jiguet, F., Albert, C .H., Thuiller, W ., 2012. Selecting ps eudo -absences for species distribution models: how, where and how man y ? Methods Ecol. Evol. 3, 327-338. Bell, D.M., Bradford, J.B., Lauenroth, W.K., 2014. Mountain landscapes offer few opportunities for high-el evation tree species migration. Glob. Change. Biol. 20, 1441 - 1451. Bellamy, C., Scott, C., Altringham, J., 2013. Multiscale, presence ‐ only habitat suitability models: fine ‐ resolution maps for eight b at species. J ournal of Applied Ec ology 50, 892 - 901. Brotons, L., Thuiller, W., Ara újo, M.B., Hirzel, A.H., 2004. Presence ‐ abs ence versus presence ‐ only modelli ng methods for predicting bird habitat suitabilit y. Ecograph y 27, 437-448. Congalton, R .G., 1991. A review of assessing the acc urac y o f classifications of remotely sensed data. Remote sensing of e nvironment 37, 35 -46. Conlisk, E., S y phard, A.D., Franklin, J., Flint, L., Flint, A., Regan, H., 20 13. Uncertainty in assessing the impacts of global change with coupled d y namic species distribution and population models. Glob. Change. Biol. 19, 858-869. de Chazal, J., Rounsevell, M.D.A., 2009. Land-use and climate change wit hin assessments of biodiversity change: A review. Global Environmental Change 19, 306-315. Eastman, J.R., 2003. IDRISI Kilimanjaro Tutorial., htt p://gisgee k.pdx.edu/G424- GIS/KilimanjaroTutorial.pdf . EEA, 2010. The European environment – state and outlook 2010: Synthesis. Esteve-Selma, M., Martínez-Fernández, J., Hernández-García, I ., Montávez, J., L ópez- Hernández, J., Calvo, J ., 2012. Potential effects of climatic change on th e distribution of Tetraclinis articulata, an endemic tree from a rid Mediterra nean ecosystems. Climatic Change 113, 663-678. Eurostat, 2012. Agri-environmental indicator - land use change, http://epp.eurostat.ec.europa.eu/statistics_explaine d/index.php/Agri- environmental_indicator_-_land_use_change#Further_information . Evans, M.R. , B enton, T.G., Grimm, V., Lessells, C.M., O’Malley , M.A., Moust akas, A., Weisberg, M., 2014. Dat a availability and model complexit y , generality, and utilit y : a reply to Lonergan. Trends in Ecology & Evolution 29, 302-303. Evans, M.R., Bithell, M., C ornell, S.J ., Dall, S.R.X., Díaz, S., Emmott, S., Ernande, B., Grimm, V., Hodgson, D.J., Lewis, S.L., Mace, G.M., Morecroft, M., M oustakas, A., Murphy, E., Newbold, T., Norris, K.J ., Petchey , O., Smith, M., Travis, J.M.J., Benton, T.G., 2013a. Predictive sy stems ecology . Proceedings of the Royal Society B: Biological Sciences 280. Evans, M.R., G rimm, V., Johst, K., Knuutti la, T., de Langhe, R., Lessells, C.M., Merz, M., O Malley, M.A., Orzack, S.H., Weisberg, M., Wilkinson, D.J., Wolkenhauer, O., Benton, T.G., 2013b. Do sim ple models lead to g enerality in ecolog y ? Trends in ec ology & evolution 28, 578-583. Ferrarini, A., Rossi, G., Mondoni, A., Orsenigo, S., Prediction of climate warming impacts on plant species could be more complex than expected. Evidence from a case study in the Himalaya. Ecological Complexity. Fielding, A.H., Bell, J.F., 1997. A review of methods for the assessm ent of prediction errors in conservation presence/absence models. Environmental conservation 24, 38 -49. Gonzalès, R., Cardille, J.A., P arrott, L., 2013. Agent-based land-use models and farming games on the social web — Fertile ground for a collaborative future? Ecolog ical Informatics 15, 14-21. Hadjinicolaou, P., Giannakopoulos, C., Zerefos, C., Lange, M.A., Pashiardis, S., Lelieve ld, J., 2011. Mid-21st century climate and weather ex tremes in C y prus as p rojected b y six regional climate models. Regional Environmental Chang e 11, 441-457. Hannah, L., 2010. Climate change biology . Academic Press. Hijmans, R .J., Cameron, S.E., Parra, J.L., Jones, P.G., Jarvis, A., 2005. Very high resolution interpolated climate surfaces fo r global land areas. International journal of climatolog y 25, 1965-1978. IPCC, 2013. Scenario Process For AR5., http://sedac.ipcc- data.org/ddc/ar5_scenario_process/index.htm . Iverson, L.R., Prasad, A.M., 2002. Potential redistribution of tree species habitat under five climate change scenarios in the eastern US. Forest Ecology and Management 155, 205 - 222. Ji, Z., Li, N., Xie, W., Wu, J., Zhou, Y., 2013. Comprehensive ass essment of flood risk using the classification and regression tree method. Stoch Environ Res Risk Assess 27, 1815 - 1828. Jiménez-Valverde, A., Lobo, J .M., Hortal, J., 2008. Not as good as t hey seem: the importance of concepts in species distribution modelling. Diversity and Distributions 14, 885 -890. Landis, J.R., Koch, G.G., 1977. The measurement of observer agreement for c ategorical data. biometrics, 159-174. Maher, S.P., Randin, C. F., Guisan, A., Drake, J.M., 2014. Pattern-rec ognition ecological niche models fit to presence-only and presence – absence d ata. Methods Ecol. Evol. 5, 761-770. Martin, Y., Van D yck, H., Dendoncker, N., Titeux, N., 2013. Testing instead of assuming the importance o f land use change scenarios to model species dist ributions under climate change. Global Ecol. & Biogeog. 22, 1204-1216. Matias, M.G., Gravel, D., Guilhaumon, F., Desjar dins-Proulx , P., Loreau, M., Münkemüller, T., Mouquet, N., 2014. Estimates of species extinctions from species – area relationships strongly depend on ecological context. Ecography 37, 431-442. Moustakas, A., Evans, M.R., 2013. I ntegrating Evolution into Eco logica l Modelling: Accommodating Phenotypic Changes in Agent Based Models. PloS One 8, e71125. Moustakas, A., Sil vert, W., Dimitromanolakis, A., 2006. A spatially ex plicit learning model of migratory fish and fishers for evaluating closed are as. Ecol. Model. 192, 245 -258. Myers, N., Mit termeier, R.A., Mit termeier, C.G., Da Fons eca, G.A., Kent, J., 2000. Biodiversity hotspots for conservation priorities. Nature 403, 853-858. Pearson, R.G., Da wson, T.P., 2003. Predicting the impacts of climate change on the distribution of species: are bioclim ate env elope models useful? Global Ecol. & Biogeog. 12, 361-371. Pellissier, L., Bråthen, K.A., Vittoz, P., Yoccoz, N.G., Dubuis, A., Meier, E.S., Zimmermann, N.E., Randin, C.F., Thuiller, W., Garraud, L., 2013. Thermal niches are more conserved at cold than warm limits in arctic ‐ alpine plant species. Global Ecol. & Biogeog. 22, 933-941. Rodríguez-Rey, M., Jiménez-Valverde, A., Ac evedo, P., 2013. Species distribution models predict range ex pansion better than chance but not better than a simple dispersal model. Ecol. Model. 256, 1-5. Rosenzweig, C., Casassa, G., Karol y , D.J., I meson, A., Liu, C., Menzel, A., Rawlins, S., Root, T.L., Se guin, B., Tryjanowski, P ., 2007. Assessment of observed changes and responses in natural and managed s y stems, In: Parr y , M. L., Canziani, O.F., Palutikof, J.P., Linden, P.J .v.d., Hanson, C.E. (Parry, M. L ., Canziani, O.F., Palutikof, J.P., Linden, P.J.v.d., Hanson, C.E.(Parry , M.L., Canziani, O.F., Palutikof, J.P., L inden, P.J.v.d., Hanson, C .E.s), Climate Change 2007: Impacts, Adaptation and Vulnerability Contribution of W orking G roup II to the Fourth Assessment Report of the Intergovernmental P anel on Cli mate C hang e. Cambridge Universit y Press, Cambridge, UK, pp. 79-131. Schweiger, O., Heikkinen, R.K., Harpke, A., Hickler , T., Klotz, S., Kudrna, O., Kühn, I., Pöyry, J ., Settele, J ., 2012. Increasing range mi smatching of interacting species under global ch ange is r elated t o their ecological cha racteristics. Global Ecol. & Biogeog. 21, 88-99. Swets, J.A., 1988. Measuring the accuracy of diagnostic systems. Science 240, 1285-1293. Thuiller, W., Lavorel, S ., Araújo, M.B., S ykes, M.T., P rentice, I.C., 2005. Cli mate change threats to plant diversity in Europe. Proceedings of the National Academy of Scien ces of the United States of America 102, 8245-8250. Tonini, F ., Hochmair, H.H., Scheffrahn, R.H., De Angelis, D.L., 2014. Stochastic spread models: A compa rison between an individual-based and a lattice-based model for assessing the expansion of invasive termites over a landscape. Ecolo gical Informatics 24, 222-230. Tsintides, T., Christodoulou, C., Delipetrou, P., Georghiou, K., 2007. Th e red data book of the flora of Cyprus. Cyprus Forestry Association, L efkosia, 465. Vennetier, M., Ripert, C., 2009. Forest flora turnover with climate change in the Mediterranean region: a case stud y in Sout heastern France. Forest ecology and management 258, S56-S63. Verdin, A., Rajagopalan, B., Kleiber, W ., Ka tz, R., 2014. C oupled stochastic weather generation using spatial and generalized linear models. Stoch Environ Res Risk Assess, 1-10. Zhang, H., Jin, X., Wang, L ., Z hou, Y., Shu, B ., 2014. Multi-agent based modeling of spatiotemporal dynamical urban growth in developing countries: simulating future scenarios of Lianyungang city, China. Stoch Environ Res Risk Assess, 1-16. Figure 1 : A graphical overview of the methodology. Data inputs are shown in dashed boxes, data outputs and intermediates in solid boxes and analy sis methods with arrows. Future potential distribution based on bioclimatic space Current land use Suitable altitude Suitable soil type Current potential habitat “Theoretical” presences and pseudo-bsences Bioclimatic space CA_MARKOV (IDRISI) Future land use Future potential habitat True species’ presence points Current bioclimatic data Future bioclimatic data Current potential habitat Future potential distribution considering climate change and land use change Future potential distribution considering climate change only Multiplication in raster calculator (ArcGIS) ArcGIS Raster calculator (ArcGIS) Classification Tree (SPSS) Multiplication in raster calculator (ArcGIS) Suitable altitude Suitable soil type Figure.2 . Classification tree for Crocus cyprius. As shown, three bioclimatic variables determine the dist ribution of C. cyprius : bio8 (Mean Temperature of Wet test Quarter ), bio5 (Max Temperature of Warmest Mont h) and bio15 (Precipitation Seasonality). If bio8 is less than or equal to 3.92 °C, the species can onl y be present if bio5 is less or equal to 25.68 °C . If bio8 is more than 3.92 °C, the species can be present either when bio15 is less than or equal to 92%. Figure 3 . Classification tree for Ophrys kotschyi. A shown, five bioclimatic variables determine the distribution of O. k otschyi : bio13 (Precipitation of Wettest Month), bio10 (Mean Temperature of Warmest Quarter), bio3 (Isothermality), bio1 9 (Precipitation of Coldest Quarter) and bio12 (Annual P recipitation ). If bio13 is less than or equal to 110.5 mm, the species can be p resent either when bio10 is less than or equal to 26.65 °C and bio3 is less than or equal to 42.5, or when bio10 is more than 26.65 °C and bio19 is between 172.5 mm and 183.5 mm. If bio13 is more than 110.5 mm , the species can onl y be present if bio13 is more than 142.5 mm and bio12 is less than or equal to 622.5 mm. Figure 4. Prediction of land use in 2050, using CORINE land-cover data for 2000 and 2006 and CA-MARKOV. All land use classes and their codes are listed in Table S5 in Supp. 1. CA_MARKOV Figure 5. Potential distribution for Crocus cyprius (a ) and Ophrys kotschyi (b ) in current conditions (top), in 2050 considering only climate change (middle) and in 2050 considering climate and land use change (bottom). All land use classes and their codes are listed in Table S5 in Supp. 1. Current distribution 2050 – Climate change only (a) Crocus cyprius (b) Ophrys kotschyi 2050 – Climate an d land use change Table 1. Number of cells and percentage cover of the potential distribution of Crocus cyprius compared to current potential distribution for each scenario used, considering climate change only and considering climate and land use change. Scenario Number of cells – climate change only Percentage cover – climate chan ge only (%) Number of cells – climate an d lan d use change Percentage cover – climate and land use change (%) Α1Β 22 31 22 31 Α2 53 74 52 72 Β1 31 43 31 43 Β2Α 212 294 204 283 Current conditions 72 100 Table 2. Number of cells and percentage cover of the potential distribution of Ophrys kotschyi compared to current potential distribution for each scena rio used, considering climate change only and considering climate and land use change. Scenario Number of cells – climate change only Percentage cover – climate change only (%) Number of cells – climate and land use change Percentage cover – climate and land use change (%) Α1Β 690 31 573 26 Α2 580 26 479 22 Β1 569 26 477 22 Β2Α 974 44 906 41 Current conditions 2205 100 Supplement 1 Supplement 1A The main geomorphologic z ones of Cy prus are the Troodos Mountains in the southwest, the Pentadaktylos mountain ran ge to the no rth, the central valley of Messaoria and the coastal zone. The main t y pes of land cover includ e high forests on the Troodos and P entadakt ylos ranges, while the lower hills are dominated b y shrubs alternating with built -up area s and cultivations. The plain of Messaoria and the coastal zone are mainl y covered by cultivations and habitations, although some parts of natural or semi -natural ve getation still persist (Tsintides et al., 2007 ). C yprus is predi cted to be severel y af fected by climate change (Lelieveld et al., 2013 ). (Hadjinicolaou et al., 2011), predict the transition of the island to a warmer state for the period from 2026 to 2050, with increase in both the max imum and the minimum temperature by 1 °C to 2 °C and decline in rainfa ll by 8.2%. Table S1: Target species information (IUCN, 2013 ; Tsintides et al., 2007 ). Both are included in Annex II of the Habitats Directive (92/43/EEC), and in Appendix I of the B ern Convention Species Crocus cyprius Ophrys kotschyi Description perennial herb perennial orchid Population over 11,500 plants in the Troodos Mountains at least 30 locations around the island usually fo rming small colonies of 10 to 100 plants Soil type igneous formations limestone formations Altitude 1050 – 1950 m 0 – 900 m Land use type pine forest openings, Juniperus foetidissima forests and edges of peat grasslands phrygana and maquis, grassy slopes, field margins, sparse pine forests and moist places Threats trampling and construction, natural fir es, climate change and military construction land development, tourism infrastructure, overcollection and failure in sex ual reproduction Protection all subpopulations with in the Natura 2000 network All subpopulations in state forests or in Natura 2000 areas ----------------------------------------------------------------- ---------------------------------------------- --------------------- Supplement 1B Table S2. Bioclimatic variables used in the model BIO1 Annual Mean Temperature BIO2 Mean Diurnal R ange (Mean of m onthly (max tem p - m in temp)) BIO3 Isotherm ality (BIO2/BI O7) (* 100) BIO4 Temperature Seasona lity (standa rd deviation *1 00) BIO5 Max Temperature of Warm est Month BIO6 Min Tem perature of Coldest Mon th BIO7 Temperature Annua l Range (BI O5-BIO6) BIO8 Me an Tem perature of Wett est Quarter BIO9 Mean Tem perature of Driest Quar ter BIO10 Mean Tem perature of Warm est Quarter BIO11 Mean Tem perature of Cold est Quarter BIO12 Annual Precip itation BIO13 Precipitation of Wettest Month BIO14 Precipitation of Drie st Month BIO15 Precipitation Sea sonality ( Coefficient of V ariation) BIO16 Precipitation of Wettest Quarter BIO17 Precipitation of D riest Quarter BIO18 Precipitation of Warmest Quarter BIO19 Precipitation of Co ldest Quarter Supplement 1C Crocus cyprius can be found in areas that combine either low mean temp eratures in the wet months (under 3.92 ⁰ C), without extreme high summer temperatures (under 25.68 ⁰ C) or high mean temperature of the wettest quarter (over 3,92 ⁰ C) but relatively low precipitation seasonality (b elow 92%). Max Temperature of Warmest Month and Precipitation S easonality are expected to increase in all scenarios, but A1B and B1 predicted the largest increase (Table S3). Additi onally, Precipitation Seasonality is projected to increase in scenarios A1B and B1, while in A2 and B2A it remains at the same levels as today (Table S3). As a result, A1B has the most pessimistic projections in terms of future suitable area, while B2A the most optimistic (Table S3). Table S3. Projected values of the bioclimatic variables that deter mine the distribution of Crocus cyprius in Cyprus (Worldclim, 2013). Max Temperature of Warmest Month ( ⁰ C) Mean Temperature of Wettest Quarter ( ⁰ C) Precipitation Seasonality (%) Present conditions 24.1-36.5 1.3-13.5 77-103 Α1Β 27.3-39.9 3.1-15.3 88-109 Α2 26.9-39.4 2.9-15.1 77-105 Β1 27.2-39.7 2.7-14.9 82-108 Β2Α 26.3-38.9 2.9-15.9 75-103 Increased projected potential distribution of Crocus cyprius does not necessarily imply that the sp ecies’ present sites are safe. This becomes evident when comparing future potential distributions wit h the species’ curre nt true distribution. Although s cenario A1B predicts the smallest potential distributi on, most of the area where the sp ecies appears today is not affected and therefore does not affect the present sites. Contrastingly in scenario B1, al though the potential di stribution is predicted to be greater than that of A1B, most of the species’ present sites are not included in it. For Ophrys kotschyi the model determines that suitable areas combine: a) Low precipitation of the wettest mont h (under 110.5 mm), low mean temp erature in the summer (under 26.65 ⁰ C) and lower temper ature diurnal range compared to annual range (isothermality under 42.5%), or b) Low pr ecipitation of the wettest mont h (under 110.5 mm), high mean tempe rature in the summ er (over 26.65 ⁰ C ) and low but not extreme low precipitation in the winter (between 172.5 and 183.5 mm), or c) Very high precipitation of the wettest month (over 142 .5 mm) but low precipitation throughout the year (below 622.5 mm). See Table S4 below for more details. Table S4. Projected values of the bioclimatic variables that deter mine the distribution of Ophrys kotschyi in Cyprus (Worldclim, 2013) Isothermality (%) Mean Temperature of Warmest Quarter ( ⁰ C) Precipitation of Wettest Month (mm) Precipitation of Coldest Quarter (mm) Annual Precipitation (mm) Present conditions 29 - 46 18.1 - 27.2 69 - 234 152 - 622 326 - 1018 Α1Β 27 - 43 21.5 - 30.6 61 - 226 136 - 608 272 - 961 Α2 28 - 44 20.8 - 29.9 47 - 213 132 - 581 260 - 947 Β1 29 - 45 20.9 - 29.9 63 - 226 126 - 597 272 - 961 Β2Α 28 - 44 20.6 - 29.8 58 - 197 137 - 528 284 - 893 Table S5. Land use classes and their corresponding c odes. Reclassification code CLC CODE LABEL 1 LABEL 2 LABEL 3 1 Artificial surfaces - - 2 211 Agricultural areas Arable land Non-irrigated arable land 3 212 Arable land Permanently irrigated land 4 221 Permanent crops Vineyards 5 222 Permanent crops Fruit trees and berry plantations 6 223 Permanent crops Olive groves 7 231 Pastures Pastures 8 241 Heterogeneous agricultural areas Annual crops associated with permanent crops 9 242 Heterogeneous agricultural areas Complex cultivation patterns 10 243 Heterogeneous agricultural areas Land principally occupied by agriculture, with significant areas of natural vegetation 11 311 Forest and semi natural areas Forests Broad-leaved forest 12 312 Forests Coniferous forest 13 313 Forests Mixed forest 14 321 Scrub and/or herbaceous vegetation associations Natural grasslands 15 323 Scrub and/or herbaceous vegetation associations Sclerophyllous vegetation 16 324 Scrub and/or herbaceous vegetation associations Transitional woodland-shrub 17 331 Open spaces with little or no vegetation Beaches, dunes, sands 18 332 Open spaces with little or no vegetation Bare rocks 19 333 Open spaces with little or no vegetation Sparsely vegetated areas 20 334 Open spaces with little or no vegetation Burnt areas 21 411 Wetlands Inland wetlands Inland marshes 22 421 Maritime wetlands Salt marshes 23 Water bodies - - References: Hadjinicolaou, P., Giannakopoulos, C., Zerefos, C., Lange, M.A., Pashiardis, S., Lelieve ld, J., 2011. Mid-21st ce ntury climate and we ather ex tremes in C y prus as projected b y six regional climate models. Regional Environmental Chang e 11, 441-457. IUCN, 2013. Crocus cyprius, http://www.iucnredlist.org/details/162216/0 . Lelieveld, J., Hadjinicolaou, P., Kostopoulou, E., Giannakopoulos, C., Poz zer, A., Tanarhte, M., Tyrlis, E., 2013. Model proje cted heat extremes and air poll ution in the eastern Mediterranean and Middle East in the twenty-first century. Reg ional Environmental Change, 1-13. Tsintides, T., Christodoulou, C., Delipetrou, P ., Georghiou, K., 2007. The red data book of the flora of Cyprus. Cyprus Forestry Association, L efkosia, 465. Worldclim, 2013. Global Climate Data. Free climate dat a for ecolo gical modeling and G I S, http://www.worldclim.org .