The Potential Distribution of Tree Species in Three Periods of Time under a Climate Change Scenario

Species distribution models have become some of the most important tools for the assessment of the impact of climatic change, and human activity, and for the detection of failure in silvicultural or conservation management plans. In this study, we modeled the potential distribution of 13 tree species of temperate forests distributed in the Mexican state Durango in the Sierra Madre Occidental, for three periods of time. Models were constructed for each period of time using 19 climate variables from the MaxEnt (Maximum Entropy algorithm) modelling algorithm. Those constructed for the future used a severe climate change scenario. When comparing the potential areas of the periods, some species such as Pinus durangensis (Martínez), Pinus teocote (Schiede ex Schltdl. & Cham.) and Quercus crassifolia (Bonpl.) showed no drastic changes. Rather, the models projected a slight reduction, displacement or fragmentation in the potential area of Pinus arizonica (Engelm.), P. cembroides (Zucc), P. engelmanni (Carr), P. leiophylla (Schl), Quercus arizonica (Sarg), Q. magnolifolia (Née) and Q. sideroxila (Humb. & Bonpl.) in the future period. Thus, establishing conservation and reforestation strategies in the medium and long term could guarantee a wide distribution of these species in the future.


Introduction
In conservation biology, estimating areas of potential distribution of species through the modeling of an ecological niche, has had a variety of applications for learning the current and future state of species [1].Yet, human activities have produced significant changes in the distribution of species in the different ecosystems of the world [2], with the temperate forests of Mexico being one of the most affected, due to the excessive extraction of some species of commercial interest, mainly of the Pinus genus like Pinus durangensis Martínez, P. arizonica Engelm., P. cooperi C. E. Blanco and P. engelmannii Carr in the Mexican northwest [3,4].On the other hand, the changes in land-use and activities related to the extraction of timber species, are among the main causes of the disturbance of habitat and the loss of diversity in the temperate forests of Mexico [5][6][7][8].Also, climate change is another phenomenon to be considered when studying the distribution of plants since many habitats have undergone severe changes.Even so, the survival of some plant species is at risk because of this [9].Projecting the potential distribution of species using climatic variables becomes very important to evaluate and foresee any alteration of either natural or anthropogenic origin [1,10,11].
One way of estimating the potential distribution of a species for a future period is by modeling their ecological niche considering that the geographical distribution of a given species is not random, but obeys environmental factors such as altitude, topographic position, temperature, humidity and precipitation, among others.Likewise, the distribution of plants is associated with physiographic (aspect and slope) and edaphic factors [12][13][14].
One of the most used tools to model the distribution of species in geographic space and their environmental tolerances, is the Maximum Entropy algorithm (MaxEnt), whose predictions start with the principle of maximum probability of occurrence of the species from the presence data and the bioclimatic variables associated with each sampling site [15].
The MaxEnt method has been used successfully in studies related to the analysis of biodiversity, the conservation of the species niche, the identification of priority conservation areas and the implementation of actions focused on preventing the establishment of invasive species in a locality [16][17][18].The maps of potential distribution projected by MaxEnt and the spatial analysis by means of a GIS have proved their effectiveness by providing a reasonable approximation of the niche of the species, even with a small sample size [19][20][21][22].
An effective way to measure the effect of environmental variations on the geographic distribution of high-interest plants is to evaluate any change manifested when going from one period of time, with specific climatic characteristics, to another, such as the period of the most recent glaciation and the contemporary period.In this way, it is possible to predict a possible displacement, fragmentation or reduction of the potential area in the face of future climate change scenarios.
This study had two objectives: (i) to identify if there is a significant change in the potential distribution of 13 tree species (highly valued both economically and ecologically) native to the temperate forests of Durango, northwestern Mexico, as a function of 19 climatic variables, considering three periods of time; and (ii) to identify species with risks of decreasing their potential distribution area in the future (2080 to 2100), under a scenario of severe climate change.In this study the spatial resolution of the data was 1 km 2 .The hypothesis was that the potential distribution of each species does not change significantly between one period of time and another.

Study Region
The modeled region corresponds to a portion of the mountain system of the Sierra Madre Occidental (SMO) in the State of Durango (northwest of Mexico), between geographical coordinates (WGS 84) 26 • 50 and 22 • 17 N and 107 • 09 and 102 • 30 W, covering an area of about 6.33 million ha (Figure 1).In this region, forests of oak-pine, pine-oak, oak or pine, together with other species typical of the temperate climate, are mainly distributed.Portions of temperate mesophytic as well as tropical deciduous and semi-deciduous forests are also found on the western slope of the region [23].The annual rainfall fluctuates from 250 to 1444 mm and the mean annual temperature from 8.3 to 26.2The presence records of the 13 species were collected using two data sources.The first data-set was obtained from on-line herbarium specimens provided by The National Herbarium (MEXU http: //www.ib.unam.mx/botanica/herbario/).The second data-set was obtained directly in the field, from 1804 sampling plots, distributed systematically every five kilometers in the study area.These plots were established by the National Forestry Commission (CONAFOR) for the National Forest and Soil Inventory 2004-2009 [25].These presence records were converted to geographical coordinates and finally a single database was generated by combining both subsets of data.

Obtaining Data
A total of 13 of the most representative species of the study area were analyzed, with eight of the genus Pinus (P.arizonica Engelm.var.stormiae Martínez, P. strobiformis Engelm, P. cembroides Zucc, P. cooperi C. E. Blanco, P. durangensis Martínez, P. engelmannii Carr, P. leiophylla Schiede ex It should be mentioned that the temperate forests of the SMO are under forest management and it is possible that anthropogenic activities have altered the natural distribution of the species studied.For example, some conifer species such as Pinus durangensis Martínez, P.cooperi C. E. Blanco, P. arizonica Engelm.var.stormiae, P. teocote Schiede ex Schltdl.et Cham and P. engelmannii Carr, are extracted more than Pinus strobiformis Engelm.or Pinus leiophylla Schiede ex Schltdl.et Cham.due to having desirable phenotypic characteristics for the forest industry [26].But since there is no precise quantification of the impact of anthropogenic activities, this factor was not incorporated into the models.

Distribution Modeling
The potential distribution of each species was modeled as a function of 19 environmental variables (Table 2) which were obtained from the WorldClim database (http://www.worldclim.org/version 2) [27,28], with a spatial resolution of 1 km.To accomplish this, we used the Maximum Entropy algorithm (MaxEnt) [15], following the methodology used by several authors [17,29,30].Potential distribution models for three periods of times were generated: (i) period of the most recent glaciation (21,000 years ago), (ii) present period and (iii) the future, corresponding to the period from 2080 to 2100.For this last period, the general circulation model used was the NIES 99 (http://www.ipcc.ch/ipccreports/sres/emission/index.php?idp=35) under an A2A scenario, a very severe scenario [31].The MaxEnt method has proven its effectiveness in making predictions based on information from presence data [32][33][34][35] and whose results express the value of habitat suitability for each species as a function of probability according to environmental variables.A high value of the distribution function in a given cell suggests very favorable conditions for the presence of an analyzed species [15].
The results obtained from the MaxEnt method, were changed to Boolean layers (presence-absence data) using ArcMap software package 10.1 (http://desktop.arcgis.com/en/arcmap/),considering a cutoff threshold equal to 10% omission errors [30,36].In each model, the variables that explained 70% or more of the data variability in the principal component analysis (PCA) were included [30,37].The PCA analyses were carried out with R software [38].
For each species, presence records were divided into 75% (randomly selected) for model training and 25% for model validation, using 100 replicates in 500 iterations with different random partitions (cross-validation method).
In order to measure the degree of association or similarity between the relative abundance of one species (presence-absence data per cell of 1 km 2 ) between one period of time and another, the Phi coefficient was used [39].Finally, to identify if there are significant differences between the potential areas of one period of time and another, the average area of the 13 species was compared, using the Kruskal-Wallis test, testing the hypothesis that the average area of a period of time is similar to the other periods.The level of significance used in all the analyses was 0.05.

Model Evaluation
Models were evaluated using the area under the curve (AUC) described by Phillips et al. [15], whose values are calculated from the Receiver Operating Characteristic (ROC).When AUC showed values equal to or greater than 0.9, we considered the models to be robust; values of AUC 0.7-0.9 were considered moderately robust; and values close to 0.5 were considered not robust [33].In all cases, the logistic output format was used [34].

Results
Observing the AUC values, the models for the scenarios of the three periods were consistent and robust showing values higher than 0.90 for all species, except for Quercus arizonica and Quercus sideroxyla (Table 1).The climatic variables used to model the potential niche of the species are listed in Table 2, which mainly includes average, maximum and minimum temperatures and rainfall in specific periods.The areas projected for the three periods of time, are shown in Figures 2 and 3.
The Phi coefficients revealed a high and positive correlation between presence/absence records of the contemporary and future period with an average of 0.81 for 92% of the species studied.The highest Phi coefficients were observed between the contemporary period and the future for P. durangensis and P. cooperi with a value of 0.87 for each species, followed by P. strobiformis with a coefficient of 0.85.In contrast, the smallest absolute coefficients (three of them negative) were observed between the period of the most recent glaciation and the contemporary period, and between the period of the most recent glaciation and the future period for Pinus cembroides and strobiformis, whose Phi coefficients were less than 0.05 (Table 3).When comparing the average values of the projected surfaces of the 13 species studied, no significant differences were found among the three periods of time, according to the Kruskal-Wallis test, with a significance level of 0.05.However, observing the areas of each species separately, several of them showed changes, as was the case of Pinus arizonica, P. leiophylla, P. cooperi and Quercus grisea, whose projected areas decreased as they passed from one period to another.This was more evident for the first two species between the period of the most recent glaciation and the contemporary period (Figure 4).The potential area of some species showed a change from one period to another.For example, P. cembroides projected an increase of 19,701 km 2 from the present to the future period and P. strobiformis projected an increase of 9633 km 2 (Figures 2 and 3).
In general, the projected potential areas were variable for each species, and although in most cases there were minor changes in the projected surface, slight displacements, a possible fragmentation of the habitats or discontinuity in the distribution of species were observed from one period of time to another (Figure 4).Habitat fragmentation was clearer for P. arizonica, P. leiophylla, Q. arizonica, Q. grisea and Q. magnolifolia (Figures 2 and 3).

Discussion
MaxEnt is a modeling instrument that has gained relevance in recent years for analyzing ecological characteristics of species using presence records, which are usually collected from specimen records in plant collection centers (e.g., MEXU) or field survey records.In this study, the areas projected (km 2 ) by models, apparently did not change from one period to another for most of  of means by the Kruskal-Wallis test, using the potential areas of the 13 species of each period; means sharing a letter are not significantly different (p < 0.05).sp1: P. arizonica, sp2: P. ayacahuite, sp3: P. cembroides, sp4: P. cooperi, sp5: P. duranguensis, sp6: P. engelmanni, sp7: P. leiophylla, sp8: P. teocote, sp9: In general, the projected potential areas were variable for each species, and although in most cases there were minor changes in the projected surface, slight displacements, a possible fragmentation of the habitats or discontinuity in the distribution of species were observed from one period of time to another (Figure 2).Habitat fragmentation was clearer for P. arizonica, P. leiophylla, Q. arizonica, Q. grisea and Q. magnolifolia (Figures 3 and 4). of means by the Kruskal-Wallis test, using the potential areas of the 13 species of each period; means sharing a letter are not significantly different (p < 0.05).sp1: P. arizonica, sp2: P. ayacahuite, sp3: P. cembroides, sp4: P. cooperi, sp5: P. duranguensis, sp6: P. engelmanni, sp7: P. leiophylla, sp8: P. teocote, sp9: Q. arizonica, sp10: Q. crassifolia, sp11: Q. grisea, sp12: Q. magnolifolia, sp13: Q. sideroxyla.

Discussion
MaxEnt is a modeling instrument that has gained relevance in recent years for analyzing ecological characteristics of species using presence records, which are usually collected from specimen records in plant collection centers (e.g., MEXU) or field survey records.In this study, the areas projected (km 2 ) by models, apparently did not change from one period to another for most of the species (Figures 2 and 3).Yet, the fragmentation or displacement of areas with a suitable climate was evident in the maps for several species, as in the case of Pinus arizonica, P. cembroides, P. engelmanni, P. leiophylla, Quercus arizonica, Q. magnolifolia and Q. sideroxila (Figures 3 and 4).Considering that the real interval of climate tolerance (optimal interval) is variable between one species and another, although they cohabit in the same region [26], it is possible that there were overestimates of the models.For this reason, the small difference between the areas projected does not necessarily suggest that the species have broad climatic tolerances or that the plants have a high capacity to adapt to contrasting climatic conditions.
The observations of which species show a fragmentation or discontinuity of the areas projected are valuable findings since they can be indicators of a high sensitivity to changes in climate values.In the end, a fragmentation can trigger a reduction in the real area of a species in the medium-or long-term, a problem that has been observed in the last decades in the temperate forests of the SMO [40] and in some Mexican terrestrial ecosystems [41,42].Moreover, in the SMO, there are no studies on the real anthropogenic effects, since they could accelerate the alteration of the natural distribution of some species.For example, during the harvesting of roundwood, regeneration and other herbaceous species are damaged in many zones of the study region due to the poor design of the road network used to transport the roundwood.If the quantitative effects of anthropogenic activity were known, a specific weight could be included and assigned when making the potential distribution models.
The Phi coefficients revealed different degrees of similarity between one period of time and another.For example, when comparing the corresponding values of the contemporary period and the future period of both Quercus magnolifolia and Q. grisea, a high correlation was observed (Phi = 0.87), suggesting a high similarity in the pattern of distribution of these species in these two periods of time.Conversely, a weak correlation was observed between the values of the past and present period of both Quercus crassifolia and Quercus arizonica (Table 3).
Potential distribution models have been useful for making important decisions in the area of ecology [14,43,44], although they should be used with caution as several authors have warned [45][46][47][48], because this type of modeling has limitations, which include the iterative procedure of the MaxEnt algorithm [49] as well as the possible errors of commission and omission when using climatic factors to characterize the habitat of a species [48].A commission error related to climatic variables could occur when the affinity of a species in a site is attributable to an event or situation that cannot be conditioned by climatic variables [48] whereas an error of omission is any error related to the underestimation of the climatic range [50] which can be translated as a false prediction of absence [48].
In some cases, the areas projected by the models are larger than the real distributions of the species studied, since the models do not consider biotic interactions [32] or the effects of human activities whose impact on the distribution pattern has increased significantly in recent years [2,51].Given that there is no precise quantification of the magnitude of the impact of human activities on the presence/absence of the species, this factor was not included in the model.Likewise, the models also exclude the dynamic nature of the individuals studied, the alteration by pests or diseases, or other situations that, for a short period, could drastically alter the distribution and abundance of forest species [52].
On the other hand, changes in the fundamental niche of a species can occur as a result of the plasticity of traits of morphological, physiological or behavioral type or by the classification of genotypes.For example, they could occur by the spatial segregation of individuals with certain functional features [53][54][55].Finally, the niche of a species can expand or move, after a long period of time, as a genetic response to the new environmental conditions [56,57].
In general, the statistical results suggest that the potential areas of most species are not experiencing drastic changes in the three modeled time periods.Nevertheless, establishing strategies for restoration, conservation, and reforestation in the medium-and long-term can guarantee not only the survival, but also the wide distribution, of the species of greatest economic and ecological interest, particularly those whose potential areas were reduced, displaced or fragmented, such as Pinus engelmanni, P. arizonica, P. leiophylla, P. cooperi and Quercus grisea, whose loss of potential area before the change from one period to another was evident (Figures 2 and 3).In these models, historical factors, biotic interactions or other limitations that affect the dispersion of the plants are not incorporated, since the magnitude of impact (quantitative values) of each of them on the distribution and abundance of organisms is unknown, and thus it is impossible to assign them an appropriate weight in the models.In future studies, we foresee incorporating physiographic variables, topographic variables and soil properties, as well as testing other analysis methods.

Conclusions
The projections reported here could be useful for decision makers.For example, for a preliminary assessment of risks of the decrease, displacement, or fragmentation of the spaces with ideal conditions for the species studied.The results could also be used to define the degree of change in the climatic domain for each of these species.Incorporating other variables like soil properties (pH, electrical conductivity, texture, etc.), physiographic variables, the annual deforestation index, or any variable that could modify the natural distribution of plants, could increase the certainty of the models.The potential areas for most of the species did not vary between one period and another.Still, establishing integral management plans including conservation and reforestation strategies, at a regional level can guarantee the continued wide distribution of the species studied.Likewise, it is highly recommended that strategies designed to minimize the impacts generated by the excessive extraction of some species, such as Pinus durangensis (cataloged as Near Threatened by the IUCN), are established.
Author Contributions: P.A. conceived and coordinated the research project, statistical analysis, and the primary writing of the text.M.E.S.-M.designed the models, analyzed the data, and helped in writing the text.C.V.-E.conducted formatting, and was the advisor and text editor.F.R.-A.: was the advisor and text editor.

Figure 1 .
Figure 1.Study area, located in northwest of Mexico.The axes are geographic coordinates.

Figure 3 .
Figure 3. Potential distribution areas of Pinus species for the period of the most recent glaciation (21,000 years ago), present period and into the future (2080 to 2100).

Figure 2 .
Figure 2. Potential distribution areas of Pinus species for the period of the most recent glaciation (21,000 years ago), present period and into the future (2080 to 2100).

Figure 4 .
Figure 4. Potential distribution areas of Quercus species for the period of the most recent glaciation (21,000 years ago), present period and into the future (2080 to 2100).

Figure 3 .
Figure 3. Potential distribution areas of Quercus species for the period of the most recent glaciation (21,000 years ago), present period and into the future (2080 to 2100).

Table 1 .
List of species studied and their records used in the potential distribution models.

Table 2 .
List of climatic variables used in the models and their respective acronyms.The variables that showed the highest and significant correlation coefficients (highlighted in bold) were used to model the potential areas.

Table 3 .
Phi correlation coefficients for the presence/absence values, comparing the degree of association among the three periods of time.