Study area

The Trans-Mexican Volcanic Belt is located at coordinates 17° 30’ and 20° 25’ N, and 96° 20’ and 105° 20’ W (^{Ferrusquía-Villafranca, 1998}) (Figure 1); it has elevations ranging from 256 to 5 650 m (^{Instituto Nacional de Estadística y Geografía [INEGI], 2018}). The area’s topography is rugged and steep in the mountains (^{Ferrusquía-Villafranca, 1998}). The dominant climate is temperate (humid and sub-humid) with an annual mean temperature of 12 to 18 °C and annual precipitation of 600 to 1 500 mm (^{García, 1998}).

Figure 1 Geographic location of the Trans-Mexican Volcanic Belt. 

Species description

Abies religiosa is a monoecious, evergreen tree, with a height of up to 60 m and maximum diameters of 1.80 m. The tree has the following characteristics: greyish, rough and cracked bark 18-25 mm thick; cross-shaped branches; simple and alternate leaves 15-35 mm long by 1.5 mm wide; and solitary cones 8-16 cm long and 4-6 cm wide (^{Protectora de Bosques [PROBOSQUE], 2007}; ^{Rzedowski, 2006}).

In Mexico, the species is located in Morelos, Estado de Mexico, Hidalgo, Puebla, Michoacán, Jalisco, Colima, Guerrero, Tlaxcala, Veracruz and Mexico City (^{Martínez-Méndez et al., 2016}) (Figure 2). In the Trans-Mexican Volcanic Belt, A. religiosa is distributed in an altitudinal gradient from 2 400 to 3 500 m and is the dominant species of the vegetation type known in Mexico as oyamel (sacred fir) forest (^{Madrigal, 1967}; ^{Rzedowski, 2006}).

Figure 2 Natural distribution of Abies religiosa in Mexico. 

Abies religiosa records

A total of 1 042 A. religiosa records were downloaded from the Niche Toolbox platform (^{Osorio-Olvera, Vijay, Narayani, Soberón, & Falconi, 2017}). In order to avoid spatial autocorrelation and overfitting of the models (^{Monterrubio-Rico et al., 2016}; ^{Peterson & Nakazawa, 2008}), duplicate, poorly georeferenced, and close coordinates (distance >1 km) were removed, leaving one record for each cell of 1 km². At the end of the depuration process, 341 records of A. religiosa distributed in the study area were obtained. These were used to generate the potential distribution models for the current and historical periods.

Current and historical bioclimatic information

Table 1 shows the 19 bioclimatic variables available in WorldClim version 2.0 for the current period (1970-2000) (^{Fick & Hijmans, 2017}) and the global circulation models (GCMs) CNRMCM5 (National Center for Meteorological Research of France) and MIROC_ESM (Institute of Oceanic and Atmospheric Research, National Institute of Environmental Studies and the Japanese Agency for Marine-Earth Science and Technology) for the Middle Holocene. The latter were generated from the CMIP5 (^{Coupled Model Intercomparison Project Phase 5, 2013}) Regional Models of the Intergovernmental Panel on Climate Change (IPCC). All variables had a spatial resolution of 1 km² and were in Geotiff (TIF) format. In addition, the variables soil type (SUE; scale 1: 250 000) and elevation (DEM; digital elevation model, 90-m resolution) were included, which were downloaded in vector and raster format from the INEGI platform (^{INEGI, 2014}, ²⁰¹⁸). At the end, all variables were homogenized to ASCII format with a scale of 1 km².

Selection of variables

Spatial autocorrelation was eliminated through a multicollinearity analysis (^{Monterrubio-Rico et al., 2016}; ^{Peterson & Nakazawa, 2008}), where those variables with a correlation coefficient ?#8805;0.85 (^{Manzanilla et al., 2019}; ^{Monterrubio-Rico et al., 2016}) were discarded to maximize the contribution of variables in the distribution models (^{Martínez-Méndez et al., 2016}; ^{Monterrubio-Rico et al., 2016}; ^{Peterson & Nakazawa, 2008}).

Delimitation of area M

Modeling area M is the environmental geographic space where a species has been reported or where it is supposed to be in accordance with the available biological knowledge (^{Martínez-Méndez et al., 2016}; ^{Soberón & Peterson, 2005}). To delimit area M, the physiographic subprovinces were used where A. religiosa records were located: Plains and Sierras of Querétaro and Hidalgo, Chapala, Mil Cumbres, Lakes and Volcanoes of Anáhuac, Neovolcanica Tarasca and Volcanoes of Colima (^{INEGI, 2001}). Area M was used to cut the variables to the same size and thus avoid the generation of overestimated distribution areas (^{Monterrubio-Rico et al., 2016}).

Generation of distribution models

The distribution models were generated with the presence records of A. religiosa in CSV format and with the variables BIO1, BIO2, BIO3, BIO4, BIO12, BIO14, BIO15, BIO18, BIO19 in ASCII format (Table 1), soil type (SUE), and elevation (DEM).

The variables were introduced into the maximum entropy program (MaxEnt), where 75 % of the species records were used to perform the training test and the remaining 25 % for the validation test of the models (^{Martínez-Méndez et al., 2016}; ^{Monterrubio-Rico et al., 2016}; ^{Phillips, Anderson, & Schaphire, 2006}). An internal replication was applied by cross validation and the Extrapolate and Do clamping boxes were disabled to avoid overfitting of the models (^{Elith et al., 2011}). The outputs of the models were of the logistic type, which represents an environmental probability index with values from 0 to 1. Values close to 0 indicate unsuitable environmental conditions, while values close to 1 suggest excellent environmental conditions for the growth and development of the species (^{Coitiño, Montenegro, Fallabrino, González, & Hernández, 2013}; ^{Phillips et al., 2006}).

Five distribution models were generated and tested, in order to determine by means of ROC (analysis of the performance characteristics of the receiver of the area under the curve [AUC]), Partial Roc and Z tests, which model fitted the current distribution of the species (Table 2). The model with best statistical performance was overlapped on the federal protected natural areas shapefile decreed for the region (^{Comisión Nacional de Áreas Naturales Protegidas [CONANP], 2017}). The current sacred fir area within the polygonal of the protected natural areas was estimated with ArcGIS 10.3 software (^{Environmental Systems Research Institute [ESRI], 2014}).

The parameters of the current distribution model with the best statistical fit (^{Morrone & Escalante, 2016}) were transferred to the maximum entropy (MaxEnt) program to generate distribution models in the past.

Validation of the models

Models were evaluated using the AUC values of the ROC, where values between 0.7 to 0.9 are considered good and those greater than 0.9 are considered excellent (^{Coitiño et al., 2013}; ^{Peterson et al., 2011}). Nevertheless, this type of validation has been questioned for not considering true absences (^{Lobo, Jiménez, & Real, 2008}; ^{Peterson, Papes, & Soberón, 2008}), so it was necessary to perform a Partial Roc analysis in the Tool for Partial Roc version 1.0 software (^{Narayani, 2008}) to counteract AUC shortcomings (^{Peterson et al., 2008}).

The recommendations of ^{Peterson et al. (2008)} were followed by using 1 000 replicates per bootstrap between presence and habitat suitability files and establishing a 5 % omission error. The test generates values from 1 to 2, where a value with an average radius of 1.0 is equivalent to a random model (^{Garza-López et al., 2016}; ^{Lobo et al., 2008}; ^{Peterson et al., 2008}). To determine whether the models were statistically valid, a Z test was performed between the AUC ratios of Partial ROC (^{Martínez-Méndez et al., 2016}). The best model was chosen based on the highest value of the Partial ROC tests, lowest standard deviation, and a reliable Z value (P < 0.01).

The logistic output values of the models were reclassified into three habitat quality categories with equal intervals (low, moderate, and high), with the ArcMap 10.3 program reclassify tool (^{ESRI, 2014}). The values of the high habitat quality category were used to transform the continuous models to binary models (suitable-unsuitable). Based on the reclassification, the potential distribution area of A. religiosa in the Trans-Mexican Volcanic Belt was estimated for both analyzed periods.

Projections of the models with the best statistical performance for each period were visualized and drawn in potential distribution maps with the ArcMap 10.3 program (^{ESRI, 2014}). Finally, the most important variables of each analyzed period were determined with the Jackknife test, which allowed quantifying the contribution of the environmental variables in the distribution models (^{Phillips et al., 2006}).

Historical and current potential distribution

Given the recent diversification of the genus Abies in the world (^{Xiang et al., 2015}), it is very likely that climatic conditions of 6 000 to 12 000 years ago were warmer (+2 °C) in North America’s temperate and cold areas (^{Caballero et al., 2010}; ^{Svensson et al., 2008}). This possibly delimited the distribution of coniferous forests to the temperate and cold zones of the Trans-Mexican Volcanic Belt (^{Caballero et al., 2010}; ^{Madrigal, 1967}; ^{Martínez-Méndez et al., 2016}; ^{Rzedowski, 2006}).

According to the Jackknife test, the most important limiting environmental variables in the distribution of A. religiosa in the Trans-Mexican Volcanic Belt were the elevation (DEM with 56.5 %), cumulative annual precipitation (BIO12 with 18.4 %) and precipitation of the warmest quarter (BIO18 with 8.3 %). These three variables explained 83.2 % of the model variability. These results are similar to those found by ^{Farjon and Filer (2013)}, ^{Madrigal (1967)}, and ^{Rzedowski (2006)}, who indicate that constant annual humidity and abundant rainfall in summer are essential for the distribution of this fir. The present study found that the ideal altitudinal gradient for the growth and development of A. religiosa was from 2 129 to 3 687 m. These results are similar to those reported in the literature, since previous studies mention that the distribution of A. religiosa ranges from 2 400 to 3 660 m (^{Calderón de Rzedowski & Rzedowski, 2001}; ^{Madrigal, 1967}; ^{Nieto de Pascual-Pola, 1995}, ^{Rzedowski, 2006}). It should be noted that, according to the results of this study, the distribution of the sacred fir forest can be found 271 m below the minimum elevation and 27 m above the maximum elevation reported in the literature. Regarding cumulative annual precipitation, the current models estimated an average of 1 125 mm. This value was slightly higher than that reported by ^{Madrigal (1967)} and similar to that mentioned by ^{Rzedowski (2006)} (1 000 mm and >1 000 mm, respectively). ^{Martínez-Méndez et al. (2016)} and ^{Sáenz-Romero et al. (2012)} mention that precipitation during the summer months is a determining factor in the distribution of the species, which agrees with what was obtained in the present study. However, dendrochronological studies carried out for the species within the study area indicate that winter-spring precipitation is of high relevance for the growth of A. religiosa (^{Cerano-Paredes et al., 2014}; ^{Huante, Rincón, & Swetnam, 1991}). This differs from what was obtained in this study, since the evaluated models did not consider the winter-spring precipitation as a key factor in sacred fir distribution. However, this variable could be of greater relevance than precipitation in the summer months during the radial growth stage of the species.

On the other hand, the most important variables in the prediction of historical distribution (MIROC_ESM model) were annual mean temperature (BIO1 with 65.9 %), cumulative annual precipitation (BIO12 with 18.6 %), and diurnal temperature oscillation (BIO2 with 4.6 %); these three variables accounted for 89.1 % of the model’s variation. For the CNRMCM5 model, the most important variables were annual mean temperature (BIO1 with 62.4 %), cumulative annual precipitation (BIO12 with 16.2 %), and diurnal temperature range (BIO2 with 6.5 %), which contributed with 85.1 % of the model. The results of the present study agree with those of ^{Madrigal (1967)} and ^{Rzedowski (2006)}, who mention that diurnal temperature range (11 to 16 °C) is a determining factor for the growth of A. religiosa. Therefore, it was corroborated that this variable has been a limiting factor for 6 000 years (10.6 to 17.2 °C) in the distribution of sacred fir in the Trans-Mexican Volcanic Belt. With respect to the annual mean temperature, the values between 4.6 and 15.7 °C are congruent with those reported by ^{Rzedowski (2006)}, who mentions that they fluctuate between 7 and 15 °C. Otherwise, according to the climatic conditions estimated by the global circulation models, annual precipitation during the middle Holocene was 80 to 224 mm higher than today’s and the annual mean temperature was 1 °C colder and not warmer as mentioned by ^{Svensson et al. (2008)} and ^{Caballero et al. (2010)}.

Because the estimated areas and environmental conditions described in the literature and in this study showed similarities, it can be assumed that the potential distribution of A. religiosa has remained stable over time (according to the theory of ecological niche conservatism; ^{Peterson, 2011b}). Thus, the evidence found in this study seems to indicate that the distribution of the species in the Trans-Mexican Volcanic Belt has not changed in 6,000 years. This would prove (with the limitation that there are no palynological records) that the ecological niche of A. religiosa has been maintained during the evaluated periods. It should be noted that, added to the fact that the current and historical distribution did not change significantly, it was possible to determine that sites such as Nevado de Toluca and Mexico City had larger sacred fir areas during the past: +2 658.4 ha and +20 129.7 ha, respectively (Table 4; Figure 5A). Consequently, it can be inferred that the climatic conditions of 6,000 years ago were colder and wetter in the temperate and cold mountainous areas of the Trans-Mexican Volcanic Belt, which caused the sacred fir distribution to remain there until the present (^{Madrigal, 1967}; ^{Martínez-Méndez et al., 2016}; ^{Ramírez-Barahona & Eguiarte, 2013}; ^{Rzedowski, 2006}).

Figure 5 Potential distribution (current and historical) of the sacred fir (Abies religiosa) in Nevado de Toluca (A) and Mexico City (B).