Study area and species studied

Mexico is located in the northern hemisphere, has an area of 1 964 375 km² (^{INEGI, 2009}), is distributed on both sides of the Tropic of Cancer, with elevations ranging between 1500 and 5000 masl; sub-humid warm climate predominates ((Am (f)) with temperatures from 2 to 30 °C, and annual precipitation from 0 to 4,000 mm (^{García, 1998}).

The selected species (symbol - number of records) were (Table 1): Pinus arizonica Engelm. (Par - 541), P. ayacahuite Ehren. (Pay - 738), P. cembroides Zucc. (Pce - 1,244), P. devoniana Lindl. (Pde - 239), P. douglasiana Martínez (Pdo - 299), P. durangensis Martínez (Pdu - 1,217), P. engelmannii Carr. (Pen - 990), P. greggii Engelm. (Pgr- 33), P. hartwegii Lindl. (Pha - 64), P. herrerae Martínez (Phe - 338), P. lawsonii Roezl. (Pla - 61), P. leiophylla Schltdl. et Cham. (Ple - 906), P. lumholtzii Robins & Fern. (Plu- 812), P. maximinoi H. E. Moore (Pma- 76), P. montezumae Lamb. (Pmo- 127), P. oocarpa Schiede (Poo- 1,028), P. patula Schl. et Cham. (Ppa - 130), P. pringlei Shaw. (Ppr - 76), P. pseudostrobus Lindl. (Pps - 473) and P. teocote Schltdl. et Cham. (Pte - 830). The altitudinal distribution of many of these species ranges between 50 to 3000 masl (^{Eguiluz, 1982}), and some can even reach the upper tree vegetation limit of at 4,000 masl (^{Yeaton, 1982}). In the pine forests, the average annual temperature fluctuates between 6 and 28 °C and the precipitation ranges from 350 to more than 1000 mm year^-1 (García, 1988).

Databases: biological, topographical and climatic

From the WorldClim website, 19 geographic covers of climatic variables were obtained in 1 km² resolution raster format (Table 2), information derived from temperature and precipitation from 1950 to 2000 (^{Hijmans et al., 2005}). Altitude, exposure and slope were calculated in ArcMap 10.2 from the Mexican Elevations Continuum 3.0 (CEM 3.0) from INEGI. A total of 10 222 biological presence records (geographic coordinates) of the species were obtained from the National Forest and Soil Inventory database (INFyS 2009 - 2014) provided by the National Forestry Commission (CONAFOR). The data were verified so that there were no repeated records and that their distribution coincided with forested areas.

Modeling and species distribution validation

To generate the potential distribution models, MaxEnt version 3.3.3k (Maximum Entropy Species Distribution Modeling) was used (^{Phillips et al., 2006}), which requires only presence data of the species. In order to generate the each species model, all Bioclim and topographic layers were used, 500 iteration were set on the software on the logistics type output, this format can be interpreted as the probability estimate (from 0 to 1) of species conditioned to environmental variables (^{Phillips and Dudík, 2008}), 50 % of the total records were used for training and the rest for testing. The models were evaluated according to their response curve tests (omission / commission analysis) and sensitivity Recipient Operation Curve (ROC) by the - Area under the curve (AUC) (^{Elith et al., 2006}, ^{Aguirre and Duivenvoorden, 2010}) and the Jackknife tests provided by the program to assess the effect of each variable in the model (^{Hijmans et al., 2005}; ^{Ramos-Dorantes et al., 2017}). The models obtained via MaxEnt were reclassified on ArcMap 10.2, to obtain surfaces at different levels of probability and thereby obtaining the bioclimatic profile of each species. A principal component analysis was performed using the standardized contribution matrix of each bioclimatic and topographic variable to know its behavior according to Bioclim variables.

Evaluation of prediction models

The statistical test of the area under the ROC curve shows that the prediction models are sensitive for species with less than 340 records (Figures 1A, B, C, E, F, G, H, K, L and O), observing that the training and test curves have a high dispersion of the prediction threshold, that is, overestimate or underestimate the projections. For species with more than 470 only-presence records, the adjustment was good, since in the omission / commission analysis (Figures 1D, I, J, M, N, Ñ, P, Q, R and S) the omission rate is close to the estimated omission (diagonal line) of both, the training sample and the test sample. The evaluation of the predictions from the models using the ROC technique is accepted as a standard method and provides a simple measure of the model's performance (^{Benito de Pando and Peñas de Giles, 2007}, ^{Aguirre and Duivenvoorden, 2010}).

Figure 1 Analysis of the adjustment of the potential distribution models for 20 pine species in Mexico, using the ROC technique. 

The analysis of the results shows that both, the AUC value and the relative probability value (RPV), significantly decrease (p≤0.001) as the number of species records increases (Figure 2A and 2B). ^{Stockwell and Peterson (2002)} showed that the modeling of any species depends on the quantity and quality of the data; of the latter, the geographic precision of its location suggest that, for the prediction to be robust, a minimum of 50 records should be used. Although species modeling with less than six records was published (^{Aguirre and Duivenvoorden, 2010}), these are reported as non-significant. In our study, the records varied from 33 in P. greggii to 1,244 in P. cembroides (Figure 1), and it is clear that the modeling of the first species might not be robust due to having few records, and the rest of the species presents more than 60 data. It follows that the fit of the model increases with fewer registers, but underestimates and overestimates the predictions (Figure 2A). ^{Elith et al. (2006)} showed that parameter estimation in MaxEnt is critical when the sample size is small, and vice versa.

Figure 2 Relationship between values of the area under the curve, AUC (A), relative value of probability (B) and number of records. Trend of AUC values with respect to the number of records of each species (C). 

The average value of the AUC of all the species was 0.959, with P. pringlei 0.982, P. patula 0.980 and P. maximinoi 0.971 as the species with the highest adjustment, contrary to that observed with P. oocarpa 0.922, P. leiophylla 0.942 and P. teocote 0.945 (Figure 2C). Even the latter value is very acceptable, since ^{Phillips and Dudík (2008)} established that the closer to 1 the value of the AUC is the better model. In addition, this average exceeds the value 0.75 suggested by ^{Elith et al. (2006)} and ^{Aguirre and Duivenvoorden (2010)}, who assure that AUC values greater than 0.75 indicate that predictions based on only presence data are sufficiently precise to establish management plans.

The modeling of conifers potential distribution for this study using MaxEnt can be considered reliable. ^{Aguirre and Duivenvoorden, (2010)} reported an average AUC of 0.924, lower than the average of our study (0.959) but varying from 0.999 to 0.766. In P. pringlei, P. patula, P. maximinoi, P. teocote, P. leiophylla and P. oocarpa studied by these authors, the modeling of all of them improved in our study with AUC ranging from 0.928 to 0.967. With 391 records P. herrerae, ^{Ávila et al. (2014)} reported an AUC value of 0.973, similar to our study (Figure 1). These variations in model fit and prediction zones are due to the number of records used and their geographic location. In 11 modeled species by ^{Ramos-Dorantes et al. (2017)}, and considered in our study, the average AUC was 0.966, slightly higher than that in our study (0.959). In 54.5 % of them, our research improved modeling: P. ayacahuite, P. greggii, P. hartwegii, P. patula, P. pseudostrobus and P. teocote. The model adjustment is not always good, in Pinus oocarpa Schiede ex Schltdl. the AUC value was only 0.548, and 0.445 in Abies guatemalensis Rehd. (^{Cruz-Cárdenas et al., 2016}), but using other bioclimatic variables, a value close to 0.5 is really a random value (^{Felicísimo et al., 2012}).

Environmental variables that model the distribution of the species

When considering the first three variables in each species model, the ones that contributed the most to predicting their distribution are altitude, Bio10 and Bio1, in 10, 9 and 6 species, respectively (Table 3). However, a single variable (Bio4) can explain 57.4 % of the potential distribution of P. lawsonii. The bioclimatic profile is different between species (Table 3) and depends on its geographical distribution in the country. In this regard, ^{Eguiluz (1982)} associated the distribution of 77 taxa of the Pinus genus throughout Mexico, and showed that several of them depend on: 1) maximum temperature, 11 species where recorded up to 45 °C, including P. pinceana, P. pseudostrobus var. Stevezii Mart., P. greggii and P. oocarpa; 2) minimum temperature, up to -23 °C, where P. ayacahuite, P. brachyptera, P. lumholtzii, P. engelmannii and Pinus arizonica grow; 3) precipitation, P. radiata var. binata Lem and P. greggii in 100 and 2900 mm of annual precipitation; 4) altitude, P. muricata D. Don and P. hartwegii Lind. grow between 120 to 4000 masl.

The principal component analysis shows two clearly defined groups: a compact group of species (Pinus cembroides, P. durangensis, P. arizonica, P. lumholtzii, P. leiophylla, P. herrerae, P. leiophylla and P. ayacahuite) were especially distributed to Bio1, Bio11, Bio18 and Bio19; and a second group highlighting P. devoniana, P. douglasiana and P. oocarpa, depending on Bioclim variables derived from precipitation, Bio13, Bio15 and Bio16 (Figure 3A). Pinus greggii showed a different response to other species, and its distribution depends on the driest month's precipitation (Bio14) and exposure. When considering a greater than 10 % contribution, altitude appears in the model structures in more than 55 % of the species, followed by Bio10 (50 %) and Bio1 (35 %); the variables that do not explain more than that value are Bio6, Bio9, Bio15, Bio17, Bio19 and exposure (Figure 3B).

Figure 3 Principal component analysis of the modeling contribution of bioclimatic and topographic variables in 20 Mexican conifers species (A) and frequency in the models (B). The symbols are explained in Materials and Methods. 

Similar to our study and derived from the Jackknife test, ^{Téllez et al. (2005)} reported that the variables that explain P. arizonica, P. devoniana, P. durangensis and P. pseudostrobus distribution are Bio1 and Bio12. In 12 conifers species evaluated by ^{Ramos-Dorantes et al. (2017)}, Bio1 was the most important variable, 11 of these species included in this same study responding to this same bioclimatic variable. For Pinus herrerae Martínez, the maximum temperature and altitude (1985 - 2227 masl) determine its distribution (^{Ávila et al., 2014}), as well as in our results, which does not explain conifer species distribution. Research by ^{Cruz-Cárdenas et al. (2016)} at a regional scale, on eight species included in our study, contrasts with our data: Bio6 explains 92.3 % of the distribution of P. douglasiana, 76.7 % of P. lawsonii, 92.4 of P. pseudostrobus and 92.3 % of P. teocote; and agree with P. montezumae (48.9 %), when observing that altitude determines its distribution.

Surface and potential distribution maps

The three species with the largest potential distribution area estimated by MaxEnt according to a probability greater than or equal to 0.70 are: P. montezumae, P. devoniana and P. pseudostrobus with 14,744.8, 14,436.1 and 11,594.8 km² (Figures 4 and 5O, L and Ñ), while those with the smaller potential area were P. durangensis, P. lumholtzii and P. engelmannii with 0.1, 275.9 and 819.1 km² (Figures 4 and 5N, J and M). Oaxaca (+) and Durango (*) are the states that register ideal area for nine species (Figure 4), Chihuahua for seven, Guerrero six, Michoacán and Puebla five (Figure 4), whereas seven states record conditions only for one species, each different. Of the registered species 50 % has a probability value greater than 0.9 (Figures 5A, B, C, E, F, G, H, L, O and S), and only P. durangensis does not have a surface within this range (Figure 5M).

Figure 4 Estimated area by species in a probability threshold greater than 0.70 and the three states of the country with highest conditions for the suitability of each species. 

Figure 5 Potential distribution maps of twenty Pinus species in Mexico. 

The comparison of the potential distribution areas of coniferous species in Mexico, between studies, would be ambiguous because they may represent different scales of information and spatial analysis (^{Ávila et al., 2014}; ^{Cruz-Cárdenas et al., 2016}; ^{Ramos-Dorantes et al., 2017}). Besides, ^{Téllez et al. (2005)} and ^{Aguirre and Duivenvoorden (2010)} report distribution areas at the species and taxa level without indicating records number, and present areas in different probability ranges. Of these latter authors, the bioclimatic profile and geographic distribution is similar to that found here for P. arizonica, P. devoniana, P. durangensis and P. pseudostrobus. ^{Stockwell and Peterson (2002)} agree that the number of records is relevant in the quality of the prediction models and, therefore, the estimation of distribution areas is important, since ^{Cruz-Cárdenas et al. (2016)} indicate that species such as P. leiophylla and P. teocote under different climate change scenarios will reduce their distribution area. The delineation of the geographical distribution of the species studied here perfectly includes the five geographic zones reported by ^{Sánchez-González (2008)}, and are useful in the assessment of the species distribution due to the effects of climate change (^{Phillips and Dudík, 2008}).