Introduction 

Natural forest habitats dominated by the genus Pinus exhibit wide biological diversity and provide ecological, economic, and social benefits, such as hydrological cycle regulation, water production, carbon sequestration, promotion of biodiversity, and scenic beauty ([6]). Many Pinus species are exploited for commercial purposes in Mexico and represent the most important source of wood, pulp, firewood, and resin, among other products ([60]).

Mexico has the second largest number of Pinus species worldwide ([22]), i.e., 52 out of the 111 known species (almost 50% - [49]). However, the majority of Pinus species in the country are restricted to very specific habitats and/or contrasting geographic environments. For example, Pinus caribaea var. hondurensis grows at sea level, while Pinus hartwegii is found up to 4000 m a.s.l., where it constitutes the upper timberline ([22]). The interaction with other species, including competition, contributes to determining the ecological distribution of P. greggii, as well as the climatic and edaphic conditions characterizing their growing sites ([8]).

Scientific evidence from Mexico indicates that the genus Pinus was always exposed to climatic changes throughout its evolutionary history. However, these changes have recently become faster due to anthropogenic activities, which caused an increase in the rate of change ([58]). Future projections foresee an increase of 2 °C in mean annual temperature by 2050, which will threaten global biodiversity ([32]). According to climate models, three different scenarios are anticipated for Pinus spp. populations in Mexico: (i) tolerate the climatic alterations through major adaptations; (ii) become locally or regionally extinct; (iii) undergo changes in their current distribution ([9], [57], [8]).

Pinus greggii Engelm. has a natural range restricted by the environmental conditions of the Sierra Madre Oriental. This species is socio-economically important in terms of firewood, fence posts, soil restoration, and the resin used to produce turpentine ([43]). In Mexico, P. greggii has been evaluated for reforestation with the aim of soil conservation and carbon sequestration ([47]). In particular, plantations of this species showed to grow well in semi-arid conditions on degraded soils ([35]). Under favorable conditions, P. greggii showed high growth rates ([59]), thus plantations of this species have also been established in other countries such as Argentina, Venezuela, South Africa, and Zimbabwe ([10]).

Knowledge of the ecological niche of P. greggii may allow environmental managers to distinguish different environmental patterns that contribute to establishment and distribution of the species, thereby obtaining useful information for conservation activities and management of genetic resources ([28]). Several tools are already available to this purpose ([11], [5]), such as cartographic representations displaying the capacity of a species to occupy a particular geographic area according to a set of variables, in conjunction with continuous or categorical characteristics of the region’s climatology, pedology, and topography ([24]).

MaxEnt is a spatial distribution algorithm widely used for assessing species’ habitat suitability at a geographic scale ([33]) and is considered a useful tool to predict the effect of climate change on future species distribution ([51], [17], [37]). Global circulation models are used to simulate future climate scenarios ([14]). Combining the current spatial distribution (ecological niche) and global circulation scenarios allows to generate a probabilistic map of future habitat suitability for a species, thus providing relevant information for its conservation and restoration programs ([57]). Indeed, MaxEnt has been successfully used to predict current and future spatial distributions of Pinaceae in Mexico ([52], [8], [17], [37]).

The negative effects of climate change (in particular drought) are expected to affect most species of the Pinaceae family in Mexico. In this context, species growing in arid or semi-arid regions, such as P. greggi, are likely to better face up the above effects and most areas of their current distribution to be preserved over time (niche conservatism theory - [62], [51]).

In this study, we analyzed geographic and environmental (climatic, topographic, and edaphic) records of P. greggii aiming to identify and delimit the most relevant environmental variables in current and future distributions, as well as to propose conservation areas within the current natural range of the species in Mexico. The specific goals were to: (i) estimate, identify and delimit the current range of P. greggii; (ii) identify the most relevant variables in the current and future (2041-2060) distribution; and (iii) propose conservation areas for P. greggii within its current natural range in Mexico.

  Methods 

Study area

The study area includes the physiographic provinces of the Sierra Madre Oriental and a portion of the Neovolcanic Transversal Belt, specifically the subprovince Plains and Sierras of Queretaro and Hidalgo, and Lagoons and Volcanoes of Anahuac (Fig. 1), between 97°-105° W of longitude and 18°-30° N of latitude ([31]). We considered the whole area including the above provinces and subprovinces as the spatial area to be modelled, based on the presence of species records therein and biological and dispersal characteristics of the species ([61]). The highest elevation point in the region is the Pico de Orizaba (5610 m a.s.l.), while the lowest is at sea level on the coast of Veracruz ([30]). The precipitation ranges from 154 to 3866 mm, with an average of 685.09 mm and a mean annual temperatures range from -2 to 28 °C ([15]).

Fig. 1 - Geographical location of the study area.

  Enlarge/Shrink   Download   Full Width  Open in Viewer

Geographic records

We used four different sources of spatial records for P. greggii: (i) 262 records were retrieved from the Global Biodiversity Information Facility ([20] - ⇒ https:/­/­www.­gbif.­org/­species/­5285216); (ii) 40 records of were retrieved from the database of the National Herbarium of the Universidad Nacional Autonoma de Mexico ([42] - ⇒ http:/­/­www.­ib.­unam.­mx/­botanica/­herbario/­); (iii) 63 records were retrieved from the Global Network of Biodiversity Information by CONABIO, that contain different national and international collections ([55]); and (iv) 73 records were obtained through dendrochronological expeditions made by personnel of the National Dendrochronological Laboratory of the INIFAP CENID-RASPA in 2018-2019.

In total, 438 records were collected and cleaned up via the Niche ToolBox platform of the National Commission for the Knowledge and use of Biodiversity ([46]) in order to eliminate double records and sites closer than 1 km each other. This step helped to avoid the autocorrelation effect and subestimation of the distribution models ([50]). Overall, 250 spatial records of P. greggii were considered for modeling after the cleanup process.

Current and future climatic variables

Current climatic information was obtained from the 19 bioclimatic layers (Tab. 1) of the WorldClim database ver. 2.0 ([15]), which contain mean climatic global information from 1970 to 2000 with a spatial resolution of 30″ × 30″ (~1 km²). For future distribution analysis, we chose the Global Circulation Models (GCMs) CNRM-CM5 and HadGEM2-ES, which are two of the more recently used GCMs in Mexico ([37]) and were generated from the Project of Regional Models CMIP-5 ([7]) of the IPCC. The bioclimatic variables of these models were downloaded with two radiative forcings of 4.5 (constant trajectories of CO₂) and 8.5 (high trajectories of CO₂) for 2041-2060 (Tab. 1), with a spatial resolution of 30″ × 30″ (~1 km²). The list of BIOCLIM current and future layers downloaded from WorlClim ver. 2.0 is presented in Tab. 1.

Topographic and edaphic variables

Topographic information was obtained from a Digital Elevation Model (DEM) with a 30 m spatial resolution downloaded from the Mexican Continuous Elevation ver. 3.0 ([31]). Elevation (ELEV) and slope (SLO) variables were rescaled to a spatial resolution of 30″ × 30″ (~1 km²) as ASCII layers; the SLO layer was generated from the topographic information of the DEM, and the ELEV layer was obtained from the elevation data of the DEM after the process of filling empty spaces. In both cases, the software ArcMap® ver. 10.3 was used ([13]).

Edaphic information was downloaded from the SoilGrids database (⇒ https:/­/­soilgrids.­org/­#!/­?layer=ORCDRC_M_sl2_250m&vector=1) at a spatial resolution of 250 m ([3]). This continuous edaphic information was developed in 2016 from the Global Soil Information Facilities (GSIF), which can be thought of as a spatial integration of a soil cartographic system at the global level ([25], [26]). SoilGrids data was generated by a model for the prediction of the physical and chemical characteristics of the soil worldwide. As for Mexico, INEGI provided series II soils profiles field data used as input for the models. Continuous variables such as coarse fragment volumetric, bulk density, absolute depth, pH, cation exchange capacity, and soil organic carbon content were extracted and adapted to an ASCII standard format with a spatial resolution of 30″ × 30″ (~1 km²).

Variable selection

For variable selection, a minimum convex polygon was generated according to the presence records of P. greggii in the study area ([16]). Later, 10.000 points of background were added and the climatic, topographic and edaphic information of each point was extracted. Environmental variables with correlation greater than r > 0.7 (p < 0.01) were eliminated to avoid the multicollinearity effect between variables ([41]). The selected environmental variables were rescaled at the same spatial resolution (30″ × 30″, ~1 km²) covering the whole study area, using the software ArcMap® ver. 10.3 ([13]).

Current distribution modeling

The MaxEnt ver. 3.4.1 algorithm was used for modelling current distribution to obtain the environmentally suitable area for the species ([52]). This algorithm was chosen because it is one of the most widely used methods for assessing species’ potential distribution and generates accurate geographic predictions based on presence records only ([11]). Seventy-five percent of records were used for training the model and 25% for the validation step. The BIO1, BIO7, BIO11, BIO15, BIO17, BIO19, ELEV, SLO, and pH variables were considered (Tab. 2).

The modeling criterion comprise internal replication by cross-validation, 1000 iterations, logistic output, 100 replicates, and a convergence threshold of 0.00001 ([52]). The “Extrapolate and Do” clamping options were deactivated, to avoid overestimation in the modeling prediction ([12]).

Model calibration was evaluated through the standardized coefficient of the Akaike information criterion (AICc), which provides model information, such as feature type and the regularization multiplier ([65]). The models showing the lowest AICc values were selected to generate the most accurate results. The calibration was carried out using the “ENMeval” library ([44]) in the R ver. 3.5.3 environment ([53]).

Modeling under future scenarios

To model P. greggii distribution under climate change scenarios, the calibration parameters and the model with the best statistical performances were transferred to the MaxEnt ver. 3.4.1 software ([45]). The estimated area (ha) of current and future distribution of P. greggii was obtained from the reclassification of the continuous values of both temporal projections (current and future) in three categories of suitability or habitat probability with equal intervals (low, medium, and high) using the “reclass” tool of ArcMap^® ver. 10.3 ([13]). The values of the high category were used as threshold cut to transform the continuous models to binary values (apt or non-apt) for every period ([37]). The conservation areas were identified using the “Intersect” tool of ArcMap, based on the simulation of current and future environmental conditions.

Model validation

The distribution models were evaluated through the statistical test of area under the curve (AUC) of the Receptor Operation Characteristics (ROC) analysis, which yields values in the range 0.0-1.0. Values from 0.7 to 0.9 indicate good model setting, while values above 0.9 indicate excellent setting ([51]). However, the utility of this analysis is strongly questioned as the algorithm uses only presence records, whereas this test requires true absences; for this reason, omission and commission errors are weighted evenly ([34]). It was necessary to perform a ROC partial test in the Niche ToolBox platform from the CONABIO ([46]) to counterbalance the AUC deficiencies. According to Peterson & Nakazawa ([50]), we generated 1000 replicates by bootstrapping ASCII data of every period and the presence records for the species, establishing a 5% omission error ([46]).

The ROC partial test generated values from 1 to 2, where a mean value of 1.0 indicates a random model ([34], [50]), [18]). A Z-test between the proportions of the AUC of partial ROC was performed to determine the statistical robustness of the models. The best model for each period was selected according to the highest value of the partial ROC, lower standard error, and statistically significant Z (p<0.01). Finally, the output of the selected models for every period were used to produce a distribution map in ArcMap ver. 10.3 software ([13]).

Relevant environmental variables

The weight of each environmental variables in the current and future modelled distribution of P. greggi over the study area were evaluated using the Jacknife test ([52]).

  Results 

Modeling current suitable area

The AUC values of the 100 replicates varied from 0.879 to 0.886 for training and from 0.797 to 0.930 for validation datasets, indicating good model performances. For P. greggii in Mexico, the best model had a partial ROC value of 1.90 (Tab. 3), an AUC of 0.881, and 0.930 for training and validation steps, respectively. The results indicate a potential current distribution of P. greggii covering an area of 617,706.04 ha (Fig. 2) within the study area. The majority of the estimated area for P. greggii is located in the states of Nuevo Leon (260,028.94 ha - 42.1%) and Hidalgo (70,762.13 ha - 11.4%).

Fig. 2 - Potential current distribution model of P. greggii in the study area.

  Enlarge/Shrink   Download   Full Width  Open in Viewer

Modeling under climate change scenarios

The AUC values obtained from the ROC test for the CNRM-CM5 RCP 4.5 model ranged from 0.901 to 0.909 for the training dataset and from 0.809 to 0.953 for the validation dataset, while the values for the RCP 8.5 projection were from 0.885 to 0.894 for training and from 0.778 to 0.942 for validation.

The results obtained for the HadGEM2-ES RCP 4.5 model showed AUC values from 0.877 to 0.888 for the training step and from 0.781 to 0.932 for the validation step, while values for the RCP 8.5 model were from 0.874 to 0.881 for training and from 0.785 to 0.936 for validation. These results allowed to classify the models of the future distribution as very good.

Relevant variables in the current and future distribution

The most relevant variables in the current distribution were BIO1, BIO11, SLO, BIO19 and BIO7, which contributed to 81.2% of the model’s variability (Fig. 3a). The relevant variables for the 2041-2060 CNRM-CM5 RCP 4.5 model were BIO1, SLO, BIO19, BIO15 and BIO11, whereas the relevant variables for the RCP 8.5 model were BIO1, SLO, BIO15, BIO19 and BIO11, with contributions of 83.4% and 87.9%, respectively (Fig. 3b).

Fig. 3 - Percent contribution to model variability of the relevant environmental variables for: (a) current distribution models; (b) CNRM-CM5 RCP 4.5 and 8.5 models; (c) HadGEM2-ES RCP 4.5 and 8.5 models.

  Enlarge/Shrink   Download   Full Width  Open in Viewer

The most relevant variables of the HadGEM2-ES RCP 4.5 model for 2041-2060 were BIO1, SLO, BIO11, BIO15 and BIO19, with a contribution of 81.5%, while for the RCP 8.5 model the variables BIO1, BIO11, SLO, BIO7 and BIO17 had an overall contribution of 85.9% (Fig. 3c).

Current and future area of P. greggii in Mexico

Tab. 4presents the estimated current and future area of P. greggii under four climate change scenarios during 2041-2060 in the province of Sierra Madre Oriental and the subprovinces Plains, Sierras of Queretaro and Hidalgo, and Lagoons and Volcanoes of Anahuac inside the Neovolcanic Transversal Belt.

The CNRM-CM5 (RCP 4.5) model foresees an increase in mean annual temperature by 0.7 °C, which will reduce the species’ ecological niche by 21.8% (relative to current area). Similarly, the CNRM-CM5 (RCP 8.5) estimates a rise of 1.1 °C in mean annual temperature and an ecological niche reduction of 17.75%. Both scenarios predict a shrinkage of the natural range of P. greggi between 2041 and 2060.

The HadGEM2-ES (RCP 4.5) model anticipated an increase of 1.5 °C in the mean annual temperature, resulting in a reduction of the ecological niche of 11.2%. The HadGEM2-ES (RCP 8.5) model indicates a reduction of the ecological niche of 7.5% relative to the current area, with an estimated 2.1 °C increase of the mean annual temperature.

The four climate change scenarios considered predict from a slight reduction (HadGEM2-ES RCP 8.5: 7.8%) to a more extensive reduction (CNRM-CM5 RCP 4.5: 21.8%) in the ecological niche of P. greggii by 2041-2060. According to these models, the increase in mean annual temperature is the main responsible for the reduction of the ecological niche of P. greggii and the shrinkage of its future natural range.

Conservation of the ecological niche

The estimated area for the ecological niche conservation of P. greggii for 2041-2060 in Mexico, according to the CNRM-CM5RCP 4.5 model was 392,923.28 ha and 366,697.07 ha with the RCP 8.5 model. Areas of 467,108.76 and 464,252.59 ha were estimated according to the HadGEM2-ES RCP 4.5 and 8.5, respectively for 2041-2060 (Fig. 4). The interpretation of these results pointed to the HadGEM2-ES RCP 4.5 as the model with the largest area of ecological niche conservation in 2041-2060, with a 75.6% increase relative to the current estimated area. The conservation of P. greggi populations in these geographical areas is crucial for in situ restoration activities of the species in the next future.

Fig. 4 - Conservation areas of ecological niche of P. greggii in the study area under climatic change scenario models for the period 2041-2060. (A) CNRM-CM5 RCP 4.5; (B) CNRM-CM RCP 8.5; (C) HadGEM2-ES RCP 4.5; (D) HadGEM2-ES RCP 8.5.

  Enlarge/Shrink   Download   Full Width  Open in Viewer

The majority of the conservation areas for P. greggii are located in the states of Nuevo Leon (223,589.71 ha) and Hidalgo (46,341.99 ha). According to the results of the current and future model HadGEM2-ES RCP 4.5, two suitable areas for conservation and restoration activities of the species are proposed: the first is located in the north, in the state of Nuevo Leon; the second is located in the center of the country, in the state of Hidalgo (Fig. 5).

Fig. 5 - Proposal of suitable areas for in situ conservation and restoration activities of P. greggii in the study area.

  Enlarge/Shrink   Download   Full Width  Open in Viewer

  Discussion 

Spatial modeling

Current and future distribution models developed using the MaxEnt algorithm in this study showed good performances, as indicated by the AUC test for the training dataset (0.88) and for the validation dataset (0.93), as well as an excellent adjustment in the partial ROC (1.85 to 1.94) and significant values of Z (p<0.01). Peterson ([51]) reports that AUC values between 0.7 and 0.9 indicate a good performance of the model, and values close to 2.0 of partial ROC are adequate with no random effects ([50], [18]).

The present results are based on 250 spatial records of P. greggii distributed over the study area. Stockwell & Peterson ([63]) suggested a minimum of 50 records to develop the species distribution analysis. Aceves-Rangel et al. ([1]) modelled the potential distribution of Pinus species using only 33 records of P. greggii with an AUC of 0.95, which is higher than that obtained in the present study. Such results could be due to the smaller number of records, the lack of debugging and calibration analysis, and the inclusion of records from the state of Chiapas, though P. greggii is endemic of the Sierra Madre Oriental and the eastern part of the Neovolcanic Transversal Belt ([54]).

A potential current area of 617,706.04 ha has been estimated for P. greggii in this study. Contrastingly, Aceves-Rangel et al. ([1]) estimated a smaller potential suitable area for the same species (550,300 ha). Such difference is remarkable (about 67,400 ha) and could be attributed to the fact that the previous study did not include edaphic variables in the model, which are considered important in modelling species potential distribution ([8], [37]).

Relevant environmental variables

The most important variable for P. greggii current distribution was BIO1, i.e., the mean annual temperature. This coincides with previous studies on Pinus habitat, which showed that this variable is a driving factor for at least ten different species ([1]). Furthermore, it confirms previous findings on the crucial role of temperature in the establishment and growth of conifer species ([64]), which has been corroborated by the association of drought index in arid zones ([36]). Indeed, elevated temperatures promote an increase in evapotranspiration and consequently a metabolic alteration that impacts the assimilation of photosynthates ([23]), leading to a reduction in tree growth ([21]). This should be taken into account especially in northern Mexico, where significant increases in temperature are forecast by the climatic models ([40]). Studies like Martínez-Méndez et al. ([39]) and Manzanilla-Quiñones et al. ([37]) found that an increase in temperature may lead to a decrease in the suitable area of Abies, Quercus, and Pinus genera.

In the study by Aceves-Rangel et al. ([1]), the variable BIO11 (mean temperature of coldest quarter) was found to be relevant for P. arizonica, which is ecologically associated with P. greggii and can be found in similar climatic conditions ([35]). Analogously, García-Aranda et al. ([17]) found that BIO11 plays an important role in the distribution of P. nelsonii in northwestern Mexico, which is also associated with P. greggii in the Sierra Madre Oriental and the Neovolcanic Transversal Belt. The mean value of the variable BIO11 was 4.3 °C in the present study, which is very similar to the mean value found for P. nelsonii (4.6 °C).

García-Aranda et al. ([17]) found that the SLO (slope) variable has great relevance in three Pinus species with distribution restricted to northeastern Mexico (P. cembroides, P. culminicola, and P. nelsonii), accounting for 21.1% of the model variability, which is close to the contribution value of SLO found in this study for current distribution model of P. greggii (18.3%). Indeed, Muñoz et al. ([43]) mentioned that P. greggii is present on slopes up to a 5%; the mean slope value found in this study is 8%.

The BIO19 variable (precipitation of the coldest quarter) was considered an important variable for at least seven Abies species in Mexico ([39]). According to the study of Aceves-Rangel et al. ([1]) BIO19 was relevant for P. lumholtzii at a national level with a lowest value of contribution to model variability of 8.2%, in comparison to the 12.8% found for P. greggii. In this study, the BIO19 averaged over the current modelled distribution was 427 mm. However, P. greggii is adapted to zones with low precipitation ranging from 293 to 747 mm ([54]).

The annual oscillation of the temperature (BIO7) for the current distribution of P. greggii in the present study had a mean value of 25.4 °C. This variable has been reported to be relevant in other similar studies as well. Hernández et al. ([27]) found that BIO7 had an importance factor of 12% in the distribution of Cedrela odorata in Mexico, slightly less than our 15.5% in the present study for P. greggii. Martínez-Méndez et al. ([39]) considered BIO7 a relevant environmental variable in the distribution of four species of Abies at the national level in Mexico.

Future scenarios

Several studies about climate change scenarios have been carried out in Mexico for Pinaceae, most of which are focused on temperate and cold climates. These studies agree with the hypothesis of a significant reduction in the natural range of Pinaceae by 2050 ([57], [8], [37]). However, this type of study has not been widely applied to Pinus species growing in arid and semiarid regions of the country.

According to the increase in temperature predicted by the CNRM-CM5 and HadGEM2-ES models with two radiative forcings (RCP 4.5 and 8.5) for 2041-2060, the ecological niche of P. greggii will decrease between 7.8% and 21.8% within its current endemic zone, but with a tendency to modify their distribution as mentioned by Gavilán ([19]).

In a modeling study of pinyon pines under climate change scenarios in Mexico, Pérez et al. ([48]) found that P. culminicola, P. johannis, and P. pinceana will undergo a decrease in the current species’ range, with a larger area predicted using RCP 8.5 compared to RCP 4.5. This situation is similar to that of P. greggii with respect to scenarios using both constant and increasing trajectories of CO₂ concentration in the atmosphere, likely because all the above species grow in similar environmental conditions, i.e., arid and semiarid regions located at the bottom of the mountains ([49]).

By modelling the future distribution of P. arizonica and P. cembroides (species ecologically similar to P. greggii), Romero-Sánchez et al. ([56]) predicted by 2050 an increase in the suitable area of 52.29% and 45.95%, respectively, compared to the current area in Sierra de Zapaliname, Coahuila. While the study of Romero-Sánchez et al. ([56]) was regionalized, the anticipated increase in suitable areas in arid climates highlights favorable effects on species such as P. greggii, P. arizonica and P. cembroides. This suggests that, despite the global effect of climate change may be negative, some species could benefit of changing environmental conditions at local level.

Conservation areas of niche

Studies conducted on the conservation of ecological niche mention that, in order to preserve their niche, species may either adapt to climatic changes through time or move to colonize new geographical areas with characteristics similar to those of their original niche ([51]). According to Booth et al. ([4]) the realized niche and parts of the fundamental niche could be measured through species distribution modeling, and the use of simulation models can assist understanding how the climatic change will affect species distribution in the future. Ecological niche modelling estimates a probabilistic index of environmental suitability over large areas for the species analysed, based on environmental variables of the sites where the species currently grow. Such index provides detailed information about the niche components, because it applies the niche theory as a multidimensional hypervolume proposed by Hutchinson ([29]) which has been complemented with the geographical analysis by Soberón & Peterson ([61]), where some or most of environmental components are preserved allowing the species to persist at different temporal scales ([51]). However, these types of studies are very scarce so far for Mexican conifers. Martínez-Méndez et al. ([39]) mention (without testing the hypothesis) that the ecological niche (as determined for Mexico) of the genus Abies remained stable over time. Manzanilla-Quiñones et al. ([37]) tested this hypothesis on Abies religiosa [Kuth] Schltdl & Cham, finding that the ecological niche of this species has been preserved since 6000 years ago in the high and humid parts of the Neovolcanic Transversal Belt. In this study, we estimated 75.6% of the conservation niche using the model HadGEM2-ES RCP 4.5, indicating a smaller conservation niche area in comparison to Manzanilla-Quiñones et al. ([37]) for A. religiosa. This could be due to a reduction in the amount of moisture and an increase in temperature of the ecological niche of P. greggii.

According to the analysis of the interaction of relevant variables in the current and future distributions of Pinus greggi, it was possible to delimit two conservation niche zones suitable for conservation and restoration activities inside the current natural distribution in Mexico. Aguirre & Duivenvoorden ([2]) modelled the current and future distribution of 56 Pinus species in Mexico, and proposed to establish new protection areas for several species of this genus in the Sierra Madre Sur, where few areas that are under regulatory protection already exist. However, their study was focused on Pinus species growing in more temperate climates, which are different from the arid to semiarid zones where P. greggii can be found. Manzanilla et al. ([38]) proposed two zones of conservation and seed production for P. pseudostrobus and P. montezumae. Although the latter species have different environmental requirements than P. greggii, the proposed zones of conservation and seed production are similar, with two overlapping zones from both studies.

According the objectives of this study, it was possible to estimate and delimit the current natural distribution of P. greggii in Mexico and the most relevant environmental variables were identified. Moreover, from the niche conservation analysis between the current and future distributions, it was possible to propose conservation areas which will likely maintain similar environmental conditions in the future.

  Conclusions 

Our results allowed to delimit the natural geographic distribution of P. greggii in the Sierra Madre Oriental and Neovolcanic Transversal Belt of Mexico. The prediction and mapping of the species distribution under four scenarios of climate change from 2041 to 2060 was also carried out.

The most important variable affecting current and future ecological niches of P. greggii in Mexico was the mean annual temperature, which is expected to increase due to climate change. A decrease in the natural range of P. greggii is predicted according to the four projected scenarios (two with constant and two with increased concentrations of greenhouse gas emissions). However, a slight increase of the suitable area outside its current range is expected by the modeled scenarios (RCP 4.5 to RCP 8.5), which will likely favor the species in the future. The regions most affected by climate change are expected to be the states of Hidalgo and Puebla, according to the four projected scenarios for P. greggii during the period of 2041-2060, though this do not imply local or regional extinction of the species therein.

Finally, the niche conservation analysis of P. greggii allowed to identify and delimit areas under similar environmental conditions (Nuevo Leon in the north, Hidalgo in Central Mexico) that could be used for in situ conservation, restoration, and forest propagation purposes.

  Acknowledgments 

This study was carried out through funds provided by CONACyT (Consejo Nacional de Ciencia y Tecnología, DF, Mexico), with the project no. 283134 “Red dendrocronológica mexicana: aplicaciones hidroclimáticas y ecológicas”.

  Author Contributions 

ARMS designed the study, carried out the analysis, and wrote the paper; JVD provided a thorough review of the paper; UMQ participated in spatial analysis and editing of the paper; JLBL and JAHH provided a thorough review and edited the paper; JEA and AHVP participated in the study design and conceptualization.

  Conflict of interest 

The authors declare no conflict of interest.

  References