Study area

The natural pine forests in Mexico are located between the following extreme coordinates: in the north of the country: to the north: 32º43’06’’ N, 114º43’22” W; to the west: 29º03’54” N, 122°10’13” W; to the east 16°06'03” N, 91°20'06” W; and to the south 11°58'07" N, 94°26´02” W. They are present in almost all the Mexican states, except for Aguascalientes, Baja California Sur, Campeche, Quintana Roo, Tabasco and Yucatán (Figure 1). The most important orographic chains where these species are located are Sierra Madre Occidental, Sierra Madre Oriental, Neovolcánico Axis and Sierra del Sur (^{Inegi, 2016}).

Simbología = Symbology; Bosque de pino = Pine forest.

Figure 1 Location of Mexico’s pine forests. 

The dominant climates of pine forests are subhumid temperate [C (w), in 25 % of the total area of the country], subhumid semi-warm [(A) C, 24.5 %], semi-arid and arid (BS, 11.5 %), subhumid semi-cold (Cb, 10.5 %), warm subhumid (Aw, 6.0 %), and other diverse climates (22.5 %) (^{Inegi, 2016}).

The altitude where these forest masses grow varies from 200 m to 4 200 m above sea level (^{Inegi, 2016}). Extreme temperatures range from 0 to 40 °C; although, most are located in sites with records from 0 to 17 °C (^{Guevara and Arroyo, 2016}; ^{Inegi, 2016}). In Mexico, pine forests are, especially, on igneous rocks, also on a regular basis they are located on gneisses and shales, less frequently on marls, sandstones, lutites and limestones. Most pine trees tolerate acid soils (^{Rzedowski, 2006}). The types of soil where they mainly develop are Litosol, Eutric Regosol, Haplic Feozem, Eutric Cambisol, Chromic Luvisol, Chromic Cambisol, Humic Andosol, Dystric Cambisol, Humic Acrisol and Ocric Andosol (^{INIFAP-Conabio, 1995}).

Data of presence of pinyon pine species

Four species of pinyon pines at risk categories were chosen, according to the NOM-059-SEMARNAT-2010: Pinus culminicola Andresen & Beaman, Pinus johannis M. F. Robert, Pinus nelsonii Shaw and Pinus pinceana Gordon. A geographic database was prepared with a total of 884 records obtained from different herbaria: Herbario CHAP de la División de Ciencias Forestales de la Universidad Autónoma Chapingo (Division of Forest Sciences of Chapingo Autonomous University); Herbario Hortorio (CHAPA) del Colegio de Postgraduados (Hortorio Herbarium (CHAPA) of the Graduate College); Herbario Metropolitano (UAMIZ) de la Universidad Autónoma Metropolitana Unidad Iztapalapa (Metropolitan Herbarium (UAMIZ) of the Autonomous Metropolitan University Iztapalapa Campus); Herbario Nacional Forestal, “Biól. Luciano Vela Gálvez” (INIF) del Instituto Nacional de Investigación Forestal, Agrícola y Pecuaria INIFAP - Cenid-Comef (National Forest Herbarium, “Biól. Luciano Vela Gálvez ”(INIF) of the National Institute of Forestry, Agricultural and Livestock Research INIFAP - Cenid-Comef); Herbario CIIDIR del Centro Interdisciplinario de Investigación para el Desarrollo Integral Regional, Unidad Durango del Instituto Politécnico Nacional (CIIDIR Herbarium of the Interdisciplinary Research Center for Regional Integral Development, Durango Campus of the National Polytechnic Institute); Herbario CFNL de la Facultad de Ciencias Forestales de la Universidad Autónoma de Nuevo León (CFNL Herbarium of the School of Forest Sciences of the Autonomous University of Nuevo León); Herbario XAL of the Instituto de Ecología, A. C. (XAL Herbarium of the Institute of Ecology, A. C.) and Herbario ENCB of the Escuela Nacional de Ciencias Biológicas of the Instituto Politécnico Nacional (ENCB Herbarium of the National School of Biological Sciences of the National Polytechnic Institute). In addition, digital herbariums and outstanding electronic pages of large collections of plants such as the Herbario Nacional de México (MEXU) of the Universidad Nacional Autónoma de México (National Herbarium of Mexico (MEXU) of the National Autonomous University of Mexico); the Herbario Virtual de la Comisión Nacional para la Conservación y Uso de la Biodiversidad (HVC) (Virtual Herbarium of the National Commission for the Conservation and Use of Biodiversity (HVC)) and the Red de Herbarios del Noreste de México (Herbanwmex) (Network of Herbaria of the Northeast of Mexico (Herbanwmex)) were consulted. Other sources of data were the Inventario Nacional Forestal y Suelos (INFyS) (National Forest and Soil Inventory (INFyS)) of the Comisión Nacional Forestal (Conafor) (National Forestry Commission (Conafor)) and scientific literature.

Climatic variables

Climatic, topographic and edaphological coverage were used, with a resolution of 926 m. 19 layers were used in raster format: temperature and precipitation of the 1902-2011 climate period, from WorldClim (Fernández et al., 2019). The altitude, exposure and slope coverage were obtained from the processing of the Continuous Model of Mexican Elevations 3.0 (CEM, for its acronym in Spanish), available on the official website of the National Institute of Statistics and Geography (^{Inegi, 2018}) and the building was generated from the vector layer of Edaphology of ^{INIFAP-Conabio (1995)}.

For the distribution model with climate change projections, the general circulation model (MCG, for its acronym in Spanish) GFDL-CM3 (Geophysical Fluid Dynamics Laboratory Climate Model version 3) of the United States of America was considered, under the scenario of greenhouse gas emissions (GHG) and Representative Concentration Pathways (CPR) 4.5 and 8.5, as well as the horizon of the medium future (2045-2069). CPR refers to radiative forcing expressed in W m^-2 and it is assumed that it varies with increasing GHG. These covers were obtained from the page of the Information Technology Unit for Atmospheric and Environmental Sciences of UNAM (^{Fernández et al., 2015}).

Modeling of the current potential distribution and under climate change scenarios

The current potential distribution models and under climate change scenarios of the stone species were generated with the MaxEnt 3.4.1 program (Maximum Entropy) (^{Phillips and Dudik, 2008}), software that estimates the probability of the presence of a species seeking the distribution of maximum entropy (as uniform as possible), based on environmental variables (covariates) and occurrence data (^{Phillips et al., 2006}). The probability is expressed in terms of "gain", which is calculated as the logarithm of the number of cells in the grid, minus the average of the negative logarithmic probabilities of the locations used to adjust the model (^{Phillips et al., 2006}; ^{Phillips and Dudík, 2008}).

During the execution of the models, the parameters predetermined by the MaxEnt program were maintained, as suggested by ^{Phillips et al. (2006)}. A replica of the model for each species was made, by random selection of 75 % of the occurrence locations for the training of the model and the remaining ones (25 %) as test data for validation.

Cut threshold for category definition

To generate the prediction models, the cut-off threshold was determined. The value established for the present study corresponds to the following categories: 0 as null, ≤ 0.30 low, 0.30 - 0.60 average and ≥ 0.60 high. The analysis and interpretation of the models generated by MaxEnt was executed with the ArcMap 10.4 software.

Evaluation of distribution models

The prediction quality of the distribution models was determined with the technique of evaluation of the characteristic operated by the receiver (ROC: Receiver Operating Characteristic), from the analysis of the area under the curve (AUC: Area Under the Curve). The AUC is the graphic representation of the discriminative capacity of a model for all possible cut-off points and can be interpreted as the highest probability that a randomly selected presence point is located in a raster cell, with a high probability value for the presence of the species.

The values of the area under the curve are categorized as follows: Excellent (AUC> 0.9), Good (0.8 <AUC <0.9), Acceptable (0.7 <AUC <0.8), Bad (0. 6 <AUC <0.7) and Invalid (0.5 <AUC <0.6) (^{Araujo et al., 2005}; ^{Elith et al., 2006}). In addition, the Confidence Interval was calculated for each of the estimates at a 95 % probability. Statistical analyzes were performed with software R version 3.4.2.

Validation of the models

The generated models suggest a good prediction of the current geographic distribution patterns and under scenarios for the four taxa of pine trees under study. The AUC values of the training and validation data were estimated above 0.9 (Table 1). According to ^{Elith et al. (2006)}, models with AUC values of 0.75 are useful for predicting species distribution; in that sense, the AUC values greater than 0.91, obtained in this work, indicate that the models are highly described by the bioclimatic variables. The highest AUC (0.99) was recorded for the most geographically restricted species, P. nelsonii. The validity of the models was checked with the statistics used. In all cases, the AUC was above 90 % and there was significance (p-value <0.05) (Table 1).

Contribution of the variables to the modeling of the species

In the current models, the variables that most influenced the modeling of the presence of P. culminicola were temperature (58.1 %) and precipitation (26.3 %); while, in P. johannis, P. nelsonii and P. pinceana were precipitation (38.4 %, 61.1 % and 42.0 %, respectively) and altitude (29.9 %, 20.2 % and 32.2 %, respectively); the rest of the proportions corresponded to the variables, but to a lesser extent. With the RCP 4.5 scenario, the proportions were: in P. culminicola the temperature, with 54.3 %; and precipitation, with 17.6 %. In P. johannis, P. nelsonii and P. pinceana corresponded to altitude (50.7 %, 40.2 % and 53.2 %, respectively) and precipitation (27.7 %, 33.7 % and 24.7 %, respectively). With the RCP 8.5 scenario, for P. culminicola were: temperature (53.3 %) and precipitation (18.2 %), in P. johannis, P. nelsonii and P. pinceana were the altitude (51.2 %, 40.4 % and 53.5 %, respectively) and precipitation (27.0 %, 33.5 % and 25.2 %); the rest of the proportions fell on the other variables, but on smaller values.

The importance of the variables in the models is coincident with that indicated by ^{Rojas et al. (2012)} and ^{Mateo et al. (2013)}, who pointed out that the climate is fundamental and widely used in modeling potential geographic distributions.

Modeling of the current potential distribution

P. pinceana has a larger area of ideal high surface and, to a lesser extent, P. johannis, P. nelsonii and P. culminicola.

Table 2 shows the potential distribution area with category: High, Medium and Low of the four pinyon pines in Mexico.

The ideal areas for P. culminicola and P. nelsonii are mainly located in the states of Nuevo León, Coahuila, Tamaulipas, San Luis Potosí and minimum areas that are limited to Querétaro, Hidalgo, Tlaxcala and Puebla. It is estimated that P. johannis develops in Coahuila, Zacatecas, Guanajuato, Tamaulipas, Hidalgo, Nuevo León and Coahuila, and to a lesser extent, in other states. The distribution of P. pinceana recorded the greatest interval of all species, since it extends to Coahuila, Durango, San Luis Potosí, Zacatecas, Nuevo León, Guanajuato and Chihuahua and, to a lesser extent, in Querétaro and Hidalgo (Figure 2).

Categoría de Idoneidad = Suitability category; Alta = High; Media = Medium; Baja = Low.

Figure 2 Current potential distribution models of endemic pine trees in Mexico. 

P. pinceana has a little more than 2 000 ha of high potential (Table 2), an estimate that differs considerably from the results obtained by ^{Aguirre and Duivenvoorden (2010)}, who modeled the distribution of 56 species of pines in Mexico, including the species mentioned, and they used the MaxEnt Maximum Entropy algorithm. They obtained an area with high fitness of 15 910 750 ha, while ^{Treviño (2018)}, calculated 1 026 800 ha.

For P. nelsonii, P. culminicola and P. johannis the results obtained in this study were 200 000 ha, 126 000 ha and 495 000, respectively, which contrasts with the results of ^{Aguirre and Duivenvoorden (2010)} who reported 11 147 600, 679 160 and 26 559 170 ha. Comparative differences in the potential areas of the studies can be attributed to the method applied in each study.

The results of modeling coincide with those of ^{Favela (2010)}, who described that P. culminicola is located at high elevations of 3 000-3 650 masl; also, ^{Contreras et al. (2012)} confirms that the climates where it grows are cold, with an average monthly temperature of 10 °C during the warmer months and an annual rainfall of 401.3 mm. ^{Romero et al. (2000)} indicate that P. johannis thrives at altitudes of 2 400 to 2 800 m with annual rainfall of 300 to 400 mm. According to Suzan et al. (2002), P. nelsonii is distributed between 1 900 and 2 300 masl, in areas with annual rainfall between 300 and 600 mm; it is located on mountain feet from arid to semi-arid environments (^{Perry, 1991}; ^{Farjon and Styles, 1997}). ^{Estrada et al. (2014)} indicated that P. pinceana is located in an altitudinal range of 1 500 to 2 300 m, in places with an annual rainfall between 350 and 600 mm (^{Martínez et al., 2015}). According to ^{Barrera et al. (2018)} the pine trees have wide latitudinal and altitudinal distribution, and occupy different climatic and edaphic regimes due to the interaction of different environmental factors.

Modeling of the potential distribution under climate change scenarios

The GFDL-CM3 model estimated increases in temperature ranges and decreased precipitation, relative to the current ranges in the potential distribution of pinyon pine trees due to the phenomenon of climate change (Table 3). According to the results, in the RCP 4.5 scenario, the simulated temperature range for the pine trees was 10 to 25 °C, while for the RCP 8.5 between 11 and 26 °C. The precipitation in RCP 4.5 was 153 to 2 584 mm and in RCP 8.5 125 to 2 605 mm. Regarding the humidity conditions, P. culminicola and P. pinceana are the species with the distribution in areas of less precipitation (Table 3).

Tables 4 and 5, respectively, show the potential distribution of the four pinion species with the GFDL-CM3 model with the climate change scenarios for the mid-future horizon (2045-2069) RCP 4.5 and RCP 8.5. In the first, in descending order of greater potential high and medium areas, P. pinceana, P. nelsonii, P. johannis and P. culminicola stand out. The first species in Coahuila, Nuevo León, Zacatecas, Guanajuato, San Luis Potosí, Tamaulipas, Hidalgo, Querétaro and Puebla. The second in Nuevo León, Coahuila, Zacatecas, Tamaulipas, Guanajuato, San Luis Potosí, Chihuahua and Hidalgo. The third in Coahuila, Zacatecas, Nuevo León and San Luis Potosí, and the fourth, in Nuevo León, Coahuila, State of Mexico, Puebla and Zacatecas.

The spatial distribution of the pinyon pines with this scenario (RCP 4.5), for the mid-2045-2069 horizon, together, showed significant surface decreases with respect to the current one. For example, P. pinceana would be reduced to 85.3 % and P. johannis 77.4 % (Figure 3). These data are consistent with the results of the work of ^{Gómez and Arriaga (2007)}, who indicate that they are species vulnerable to climate change.

Porcentaje de distribución = Percentage distribution

Figure 3 Predicted changes of the species of pine trees with reference to the current potential distribution in the 2045-2069 horizon considering the RCP scenarios 4.5 and 8.5. 

In the 8.5 RCP scenario in descending order of higher potential high and medium areas, P. pinceana, P. nelsonii, P. johannis and P. culminicola stand out. The first species would be located mainly, by its surface, in Coahuila, Nuevo León, Zacatecas, Guanjuato, San Luis Potosí, Tamaulipas, Hidalgo, Puebla and Querétaro. The second in Coahuila, Nuevo León, Zacatecas, Tamaulipas, Guanajuato, Chihuahua, San Luis Potosí, Hidalgo and Durango. The third in Coahuila, Nuevo León, Zacatecas, Puebla and San Luis Potosí; and the fourth, in Nuevo León, Coahuila, State of Mexico, Puebla, Zacatecas and Baja California.

Modeling with the RCP 8.5 emission scenario suggests losses similar to the RCP 4.5 scenario. As a whole, there were significant reductions with respect to the current potential of P. pinceana would be reduced to 85.0 % and P. johannis, 76.8 % (Figure 3).