Introduction
There is strong evidence that anthropogenic climate change is generating negative impacts on natural systems, mainly changes in precipitation patterns and increases in the Earth's surface temperature (Intergovernmental Panel on Climate Change [IPCC], 2014). These climate modifications are increasing the frequency and intensity of extreme droughts, generating massive mortality of forest stands (Hammond et al., 2022), especially at the lower altitudinal limit of the natural distribution of populations (Mátyás, 2010). This is leading to a decoupling between local plant populations and the climatic conditions to which they have adapted over time (Chen et al., 2022). In this context, high mountain forests have become particularly sensitive to these effects (Brodersen et al., 2019; Rangwala & Miller, 2012) because the long-life cycle of trees does not allow for a rapid adaptation of populations to changing climatic conditions (Brodribb et al., 2020). On the other hand, populations at lower elevations are more sensitive to moisture availability (Hartl-Meier et al., 2014), and natural migration rate to higher elevations is very low compared to the speed required to cope with climate change (six to 11 times slower than necessary; Gómez-Pineda et al., 2020).
Sacred fir forests (Abies religiosa [Kunth] Schltdl. et Cham.) are considered forest relicts belonging to high mountain ecosystems (Araiza-Olivares, 2020). In Mexico, they occupy between 0.1 - 0.5 % of the forest area (Manzanilla-Quiñones et al., 2019; Pineda-López et al., 2013). Their populations are found in isolated and fragmented patches in locations with a humid and cold microclimate. They are often restricted to the top of mountains (between 2 400 and 3 600 m) in sites with high humidity (>1 000 mm·yr-1), mean annual temperatures of between 7 and 15 °C, and deep, moist, and well-drained soils (Rzedowski, 2006).
Populations of A. religiosa often exhibit low seed viability (de Pascual Pola et al., 2003), low germination capacity (Ortiz-Bibian et al., 2019), and, during the now warmer and drier March-to-May season, scarce natural regeneration (Guzmán-Aguilar et al., 2020), and defoliation due to drought stress (Sáenz-Romero et al., 2020). Habitat modeling studies by Gómez-Pineda et al. (2020) and Sáenz-Romero et al. (2012) have projected a drastic reduction in suitable habitat for the species between 77 % and 96.5 %, respectively, determining that the climatic conditions suitable for sacred fir will be found even higher up in the mountains. This reduction in suitable climatic habitat would not only affect the vegetation but also the migratory monarch butterfly (Danaus plexippus L.) colonies, which overwinter exclusively (November-March) in the dense stands of A. religiosa (Carlón-Allende et al., 2016). Therefore, assisted migration to higher elevations will be important to realign species with climates that will be favorable to them in the future (Sáenz-Romero et al., 2012).
This report presents a follow-up, at a later stage, of an experiment whose initial results were first reported by Cruzado-Vargas et al. (2021). The objective of the study was to determine the effect of the climate transfer distance (CTD = planting site climate - provenance climate) on survival ratio, height increase, and productivity index in A. religiosa seedlings. The effect was evaluated in seedlings from 11 provenances through reciprocal transplanting at trial sites with contrasting elevations: 3 000 m as the reference site, 3 400 m to assess the feasibility of assisted migration to higher elevations, and 2 600 m to estimate the impact of climate change in warmer and drier conditions. This was conducted to generate management guidelines for making management decisions to adapt to the effects of climate change.
Materials and Methods
A total of 11 seed collection sites of A. religiosa were selected, with an altitudinal gradient from 3 000 to 3 550 m, with a 50 m altitudinal difference, in the Monarch Butterfly Biosphere Reserve (MBBR) in Estado de México (Table 1). At each site, 10 trees were selected with a minimum linear distance of 30 m between them to reduce the risk of inbreeding, and five cones were collected from each tree.
Table 1.
| Elevation (m) | Latitude (N) | Longitude (W) |
|---|---|---|
| Seed provenances | ||
| 3 552 | 19° 33' 52.6'' | 100° 13' 47.4'' |
| 3 491 | 19° 34' 04.1'' | 100° 13' 59.5'' |
| 3 457 | 19° 34' 17.1'' | 100° 14' 08.2'' |
| 3 411 | 19° 34' 25.0'' | 100° 14' 13.8'' |
| 3 364 | 19° 34' 31.7'' | 100° 14' 03.8'' |
| 3 300 | 19° 34' 46.2'' | 100° 13' 53.6'' |
| 3 233 | 19° 34' 50.9'' | 100° 13' 26.6'' |
| 3 210 | 19° 34' 52.6'' | 100° 13' 15.5'' |
| 3 143 | 19° 34' 53.0'' | 100° 12' 53.1'' |
| 3 099 | 19° 35' 12.0'' | 100° 12' 52.9'' |
| 3 003 | 19° 35' 42.6'' | 100° 12' 37.5'' |
| Common garden experiment | ||
| 3 400 | 19° 34' 23.2'' | 100° 14' 08.4'' |
| 3 000 | 19° 35' 18.2'' | 100° 10' 46.0'' |
| 2 600 | 19° 48' 29.6'' | 100° 10' 52.7'' |
In the field, three trial sites were established in common gardens at contrasting elevations: 3 400 m (Llano Grande, Ejido La Mesa, Municipality of San José del Rincón, Estado de México; MBBR core zone), to assess the feasibility of conducting assisted migration; 3 000 m (Llano Grande, Ejido La Mesa, Municipality of San José del Rincón, Estado de México; MBBR buffer zone margin), as a reference site, similar to the current optimum for the species; and 2 600 m (Tlalpujahua de Rayón, Michoacán), to assess the impact of increased temperatures due to the predicted effects of climate change (Table 1).
At each trial site, three raised growing beds were established (wooden structure of 5 m in length x 1.5 m in width x 0.4 m in height. All beds at the three sites contained the same forest soil substrate obtained at 3 000 m). The common substrate and the dimensions of the raised bed were intended to prevent plant roots from coming into contact with the local soil, enabling assessment of the influence of climatic transfer distance alone on the growth and survival of the provenances. Subsequently, plants produced in a shade house at Instituto de Investigaciones sobre los Recursos Naturales (INIRENA), in Morelia, Michoacán, Mexico were taken to the field during the rainy season (July) of 2019 and used to establish six randomized complete blocks (two blocks per raised bed) at each field site: 11 provenances per block, five seedlings per provenance per plot in line, with a spacing of 0.2 x 0.2 m. At the head of each raised bed, seedlings from mixed provenances were established as border effect. The raised beds were covered with a 35 % shade net to protect from an excess of solar radiation and extreme low temperatures (Carbajal-Navarro et al., 2019).
In December 2019, initial plant height measurements were taken with a precision of ±1 mm. Basal diameter was measured at the root collar with a digital vernier. Subsequently, plant height and survival were measured bimonthly (for monitoring purposes), and basal diameter was recorded every six months. The final measurements were taken in November 2021 (end of the growing season, and approximately 2.5 years [28 months] after planting). The height increase was obtained from the difference between the final and initial measurements, and the survival ratio was obtained per provenance per block and then per provenance per site.
To determine the aerial biomass of each plant in the field (Ba, g), an estimation was made through a destructive analysis of five seedlings per provenance that had been kept in the shade house at INIRENA, i.e., 55 seedlings in total. Total height (Alt, cm) from the root collar to apical bud and basal diameter (Db, mm) at the root collar were obtained from each plant. The plants were weighed fresh with an analytical balance and placed in paper bags. The stem, branches and foliage of each plant were placed separately in each bag, put in a drying oven at 65 °C for 48 h (until constant weight) (Domec et al., 2010), and the dry weight was recorded. With this data, a fitted regression was performed using the PROC REG procedure of SAS (SAS Institute, 2014, 2015), obtaining the following parameters:
With the aboveground biomass and survival ratio data, an average was obtained per provenance per block per site, and a productivity index was calculated based on the methodology of Sáenz-Romero et al. (2021), in which: survival ratio per provenance per block (values from 0 to 1) was multiplied by the average aboveground biomass per provenance per block (in grams of dry weight).
For estimation of the climatic transfer distance, climate data for the provenances were obtained through spline models based on Sáenz-Romero et al. (2010) and available at the website https://charcoal2.cnre.vt.edu/climate/ (Moscow Forestry Sciences Laboratory, 2023). The climate of the provenance was considered as the historical climate in which the populations have developed, i.e., the reference period 1961-1990. Temperature data loggers were installed, and rainfall was measured monthly with rain traps strategically placed in each trial site to obtain the climate of the trial sites. Thirteen climate variables were estimated for the climates of provenance (1961-1990) and of the trial sites (2019-2021). Using these climate data, a series of mixed models were tested to determine which climate variables (Table 2) had the greatest impact on survival ratio, height increase, and productivity index. These models have been used similarly by Cruzado-Vargas et al. (2021). The structure of the mixed models was as follows:
where:
Y ijkl = survival ratio, height increase, or productivity index corresponding to the l -th plant of the j -th provenance in the k -th block nested within the i -th trial site
β 0 = intercept
β 1 - β4 = parameters of the fixed effects
β 5 - β8 = parameters of the random effects
CP j = value of the climate of the j -th provenance
CTD ij = climatic transfer distance, i.e., the difference between the value of the climatic variable at the trial site minus that of the climate of the provenance, for the j -th provenance in the i -th trial site
CP ij x CTD ij = interaction between the climate value of the j -th provenance and the CTD of the j -th provenance in the i -th trial site
S i = effect of the i -th trial site.
B k (S i ) = effect of the k -th block nested within the i -th trial site
P j = effect of the j -th provenance
(S i x P j ) = interaction between the i -th trial site and the j -th provenance
e ijkl = error term
Table 2.
| Code | Unit | Definition |
|---|---|---|
| MAT | °C | Mean annual temperature |
| MAP | mm | Mean annual precipitation |
| GSP | mm | Growing season precipitation (total precipitation April-September) |
| WDSP | mm | Total warm and dry season (March-May) precipitation |
| RSP | mm | Total rainy season (June-October) precipitation |
| CDSP | mm | Total cold and dry season (November-February) precipitation |
| MTCM | °C | Mean temperature in the coldest month |
| MMIN | °C | Mean minimum temperature in the coldest month |
| MTWM | °C | Mean temperature in the warmest month |
| MMAX | °C | Mean maximum temperature in the warmest month |
| DD5 | °C | Degree-days > 5 °C |
| WDSDD5 | °C | Warm and dry season (March-May) degree days > 5 °C |
| RSDD5 | °C | Rainy season (June-October) degree days > 5 °C |
| CDSDD5 | °C | Cold and dry season (November-February) degree days > 5 °C |
| AAI | index | Annual aridity index ( |
| GSAI | index | Growing season aridity index ( |
| WDSAI | index | Warm and dry season (March-May) aridity index ( |
| RSAI | index | Rainy season (June-October) aridity index ( |
| CDSAI | index | Cold and dry season (November-February) aridity index ( |
To select the mixed models that best described the response of the dependent variables, a series of steps were followed as described below (for more detail, see Sáenz-Romero et al 2017). 1) A Spearman correlation analysis was performed between the climate variables of the source of provenance and the mean value per provenance of the response variables, selecting those with the value closest to ±1. 2) A reduced model was fitted between the climate variables for the climatic transfer distance, i.e., the term for the climate of provenance was dropped. The climate variables from the CTD with a significant p-value and with the lowest Akaike Information Criterion (AIC) value were then selected. 3) Mixed models were run using the climate provenance and climatic transfer distance variables selected in the previous steps (1) and (2), and those with a negative and significant value of the quadratic term (Leites et al., 2012), as well as the lowest AIC, were selected. 4) A total of 75 models were run, of which one was selected for each response variable with the lowest AIC value.
The PROC MIXED procedure of SAS (SAS Institute, 2014) was used for the height increase and productivity index. In the case of survival ratio, mixed models for the survival ratio analysis were fitted in RStudio software version 2022.07.2 Build 576 (RStudio Team, 2020). The ‘glmer’ function of the lme4 package (Bates et al., 2015) was used and the "binomial" family was selected since the variable type is binary (0 = dead and 1= alive). Given the binomial ‘glmer’ predict in ‘logits’, we then transformed the response to exponential values following the formula
The final model was then used to predict and plot the response curve to the climatic transfer distance for: i) the average for the species (averaging the climate of origin across all populations), ii) the population collected at a higher elevation and iii) the population collected at a lower elevation. To build such response curves, the parameters of the fixed terms from the model were used, after choosing the combination of climate variables (of the climatic transfer distance and of the provenances) that yielded the best fit.
Results
For the three response variables analyzed (survival ratio, height increase, and productivity index) the response function best explained (based on their significance of the linear, quadratic, and lowest AIC terms) was the climatic transfer distance of the dry and cold season aridity index (CDSAI: November-February). For the climate of provenance variables, this was the warm and dry season aridity index from March to May (WDSAI) for the survival ratio (Table 3) and the rainy season aridity index (RSAI: June-October; Table 4) for height increase and productivity index.
Table 3.
| Parameter of the source of variation | Estimator | p |
|---|---|---|
| AIC | 223.32 | |
| Intercept | 6.1705 | 0.0022 |
| Warm and dry season aridity index (WDSAI) | -12.38 | 0.1944 |
| Cold and dry season aridity index (CDSAI) | -0.599 | 0.8940 |
| (CDSAI)2 | -6.99 | 0.0073 |
| Climatic seed source × CTD | -0.314 | 0.9826 |
Table 4.
| Parameter of the variation source | Height increase | Productivity index | ||||||
|---|---|---|---|---|---|---|---|---|
| Fixed effects | Estimator | p | Estimator | p | ||||
| AIC | 4432.7 | 566 | ||||||
| Intercept | 32.36 | 0.1728 | 4.60 | 0.2951 | ||||
| Rainy season aridity index (RSAI) | 759.39 | 0.0828 | 184.41 | 0.0465 | ||||
| Cold and dry season aridity index (CDSAI) | 276.62 | <0.0001 | 39.14 | <0.0001 | ||||
| (CDSAI)2 | -115.64 | <0.0001 | -26.36 | <0.0001 | ||||
| Interaction climatic seed source × CTD | -4243.83 | <0.0001 | -548.55 | 0.0081 | ||||
| Site | 0 | 0 | 1 | 0 | 0 | 1 | ||
| Provenance | 17.6 | 3.6 | 0.0864 | 0.2397 | 2.6 | 0.3569 | ||
| Block (Site) | 16.2 | 3.3 | 0.1254 | 1.97 | 21.2 | 0.0579 | ||
| Interaction site x provenance | 0 | 0 | 1 | 0.71 | 7.6 | 0.2117 | ||
| Error | 457.78 | 93.1 | <0.0001 | 6.39 | 68.6 | <0.0001 | ||
*Contribution to the total variance expressed as a percentage, where 100 % is the sum of the contribution to the total variance of the random terms.
Figure 1 shows the predicted response curves, constructed using the estimated parameters based on the fixed effects of the mixed model. The predicted response curve for the survival ratio (Figure 1a) indicated that survival was very high (0.95 = 95 %, range between 0 and 1) when plants were transplanted to the coldest and wettest site at 3 400 m (an elevation difference of 400 m corresponds to approximately -2 °C, based on a lapse rate of 0.5 °C per 100 m [Körner, 2007; Sáenz-Romero et al., 2010]). Similarly, a high survival rate (91.6 %) was observed at the site located at 3 000 m. In contrast, at the site at 2 600 m, which is warmer and drier, the survival rate (5.7 %) was ~16 times lower than at the other two sites. The flattened shape of the upper part of the survival curves is due to the values range between 0 and 1, given the binomial nature of the variable (0-1). The optimum survival value found lies between -0.5 and 0.6 of the x-axis climatic transfer distance values. In this range, the sites at 3 400 and 3 000 m showed the highest values, while the site at 2 600 m was located in the extreme lower part of the right-hand side of the curve (much drier), which explains the markedly low survival value. However, when comparing the mean annual temperature and mean annual precipitation values between the reference site at 3 000 m and the lower site at 2 600 m, it was observed that the latter is 1.8 °C warmer and drier (with 14 % less precipitation). These two variables are probably responsible for the high seedling mortality found at the 2 600 m site (94.3 %).
For plant height increase (Figure 1b), the response curve indicates that, at the 3 400 and 3 000 m sites, growth among the provenances per site was very homogeneous, with an average plant height of 31 cm at the 3 400 m site, and nearly three times that (87 cm) at the reference site (3 000 m). At the 2 600 m site, however, the (few) plants that survived presented heterogeneous growth, i.e., the maximum height recorded was 88 cm for the tallest plant and 29 cm for the shortest, with an overall average of 63 cm at the site.
Regarding the response function for the productivity index, which represents the balance between growth and survival, an average value of 6.4 was recorded at the site at the highest elevation, while the reference site had the highest productivity index (14.7). The site at the lowest elevation presented different and extreme values among the plants, but the average value was 5.2, which was the lowest value of the three trial sites (Figure 1c). Figure 2 shows the plant status (survival and growth) at the three trial sites.


Figure 3 represents the change in mean annual temperature and mean annual precipitation of the three sites within the historical period 1961-1990: a) the site at the highest elevation recorded a decrease in temperature of 2.4 °C, thus being even colder than the historical period, with an increase in mean annual precipitation of 15 %; b) the reference site at 3 000 m recorded an increase of 1.7 °C compared to the historical value; that is, it has become warmer despite an 8.5 % increase in precipitation; c) the site at 2 600 m has become ~2.1 °C warmer and also slightly drier (considering only the change in precipitation, but not the change in evapotranspiration due to the increased temperatures), with a loss of 1.25 % in the annual precipitation compared to the historical value. August recorded the highest precipitation levels for all three sites.

Discussion
The balance between temperature and precipitation in the cold and dry season (November - February); i.e., the CDSAI, proved to be the most significant climate variable of climatic transfer distance for the three response variables (survival ratio, height increase, and productivity index). Winter-spring precipitation is important for the survival, growth, and productivity of several conifer species, including Pinus cembroides Zucc. (Carlón Allende et al., 2018a b), Pinus pseudostrobus Lindl., and A. religiosa, since the warm temperatures influenced by the El Niño phenomenon lead to reduced precipitation during this season in central Mexico (Carlón Allende et al., 2016). Murray-Tortarolo (2021) notes that changes in the timing and intensity of precipitation have caused seasonal extremes, with the wet season becoming wetter and the dry season drier.
At the site at 3 400 m, a high survival rate (95 %) was recorded. This was probably related to the absence of extreme heat and the 51 % increase in precipitation during the cold and dry season (November - February) compared to the historical period (1961-1990). This availability of moisture most likely conferred a greater capacity for survival on the plants (Allen et al., 2015; Mátyás et al., 2018). At the site at 3 000 m, the 22 % decrease in precipitation during the same November - February period, compared to the historical period of 1961-1990, had no effect on plant survival, possibly because precipitation during the rainy season (June - October) of the 2019-2021 trial period was higher (>15 %) compared to that of the historical period 1961-1990, thus providing the plants with sufficient soil moisture to survive the cold and dry season. In sharp contrast, at the 2 600 m site during the November - February cold and dry season, there was an increase in temperature of ~3.5 °C and a decrease in precipitation of 8 % (compared to the historical 1961-1990 values), causing survival to be ~16.7 times lower than at the higher elevation site and ~16.1 times lower than at the 3 000 m site.
It is important to highlight the value of the comparison of the effect of the balance between temperature and precipitation at the reference site at 3 000 m, which presented the best performance of the plants in terms of growth, as well as high survival (91.6 %), compared to the site at 2 600 m. Although the latter site at a lower elevation is outside the distribution of the origin of the provenances (>3 000 m), the fact that it has both higher temperatures (1.8 °C) and lower precipitation (-14 %) than the reference site at 3 000 m, allows us to visualize the magnitude of the climatic impact to be expected with a temperature increase of less than 2 °C and a decrease in precipitation of only 14 %. These results indicate a highly probable scenario for Mexico between the 2060s and 2090s (Sáenz-Romero et al., 2010); such changes in climatic variables will be sufficient to induce massive seedling mortality, as observed at the 2 600 m site (94.3 % mortality).
For plant height increase, contrasting results were found among the trial sites. For the site at the highest elevation (3 400 m), height growth values were the lowest in the entire experiment. This result is consistent with the fact that such elevation is near the natural upper altitudinal limit of A. religiosa in the studied region. This is due to the negative correlation, observed in conifers, between high growth rates and cold resistance (Rehfeldt et al., 2018); in other words, growth decreases as cold increases, and the growing season shortens with increasing elevation (Liu & El-Kassaby, 2018). Furthermore, it has been shown that temperatures below 5 °C limit root growth because the transport of the H⁺-ATPase enzyme decreases, thereby inhibiting several functions essential for cellular growth (Alvarez-Uria & Körner, 2007), thus growing degree days (DD5; daily temperature above 5 °C) have a significant influence on plant height growth. During the cold and dry season from November to February, DD5 values differed significantly among the sites: 11 DD5 were recorded at 3 400 m, 816 DD5 at 3 000 m, and 1 032 DD5 at 2 600 m. This considerable difference in DD5 may restrict growth at the highest-altitude site while benefiting plants at the lower sites. It was also observed that, at the 2 600 m site, the few plants that survived until the next rainy season exhibited the second-highest average height growth among all sites. It is likely that populations found in the driest sites (located at lower elevations within the natural distribution of the species) have a greater capacity for growth recovery following a period of drought stress compared to populations from more wetter regions (higher elevations). This has also been observed in populations of Pinus pinaster Ait. in the western Mediterranean region (Sánchez-Salguero et al., 2018), and in some forest species more frequently found in North America and Europe, such as Pinus ponderosa Douglas ex Lawson et C. Lawson, Pseudotsuga menziesii (Mirb.) Franco, Quercus alba L., Picea glauca (Moench.) Voss., Tsuga canadenis (L.) Carr., and Taxodium distichum (L.) Rich. (Gazol et al., 2017).
Regarding the relationship between growth and survival, defined as the productivity index, it can be assumed that the plants established on the site at 3 400 m grew less in part because they designated a significant amount of resources obtained from photosynthesis for storage and defense against the low temperatures (Lusk & Jorgensen, 2013; Qin et al., 2022). This is a common response in shade-tolerant and slow-growing species such as A. religiosa, which allows them to survive longer but not increase their formation of biomass (Harsch & Bader, 2011). At the 3 000 m site, a high productivity index was recorded, because the most favorable climatic conditions for the performance of A. religiosa trees are found between 3 050 and 3 250 m: Ortiz-Bibian et al. (2019) and Guzmán-Aguilar et al. (2020) respectively showed that the highest number of viable seeds and recruitment were found in the intermediate part of the natural distribution of this species. Similarly, Musule et al. (2016), reported that between 3 100 and 3 200 m, the anatomical structure of sacred fir wood features a lower content of hemicellulose (related to frost or drought resistance). As for the site at 2 600 m, the productivity index was lower. Although the surviving plants at this site showed a good performance in terms of height increase, survival was too low and was insufficient to present a high value for the productivity index. A study by Dixit et al. (2021) presented similar results, in which the survival of P. ponderosa seedlings was the poorest recorded at the lowest elevation site.
The results of this study seem to indicate that at 3 000 m is currently the site of ecological optimum for the provenances tested. The ecological optimum refers to the site where the species occurs most frequently and can be competitively exclusive (Rehfeldt et al., 2018). However, given the effects of climate change, it has been predicted that soon (by the decade centered on the year 2060), the optimal habitat for A. religiosa will be 500 m higher in elevation, reaching 4 000 m at its upper altitudinal limit (Gómez-Pineda et al., 2020). Considering that the altitudinal limit of the sacred fir in the MBBR is approximately 3 550 m, assisted migration of populations of the species to other higher peaks of the Trans-Mexican Volcanic Belt, such as Nevado de Toluca, in the Estado de México, should be regarded as a vital conservation strategy (Sáenz-Romero et al., 2012). This measure would also contribute to the adaptation and mitigation of climate change impacts.
Based on the results obtained in this study, it is recommended to implement assisted migration to higher elevation sites as a conservation measure for the species and as part of the effort to mitigate the expected effects of climate change.
Conclusions
The altitudinal movement of Abies religiosa provenances toward cooler and more wetter sites (up to ~2 °C cooler= 400 m altitudinal difference upwards) will allow the balance between temperature and available moisture to remain within the range necessary for plant survival. However, at planting elevations such as 3 400 m, growth rates will be limited-at least in the short term-because cold temperatures will act as a limiting factor. Populations of A. religiosa located at the lower limit of their distribution will increasingly suffer from water stress due to increasing temperatures. This is likely to result in mass seedling mortality, apparently induced by an increase in mean annual temperature of 1.8 °C and a decrease in precipitation of 14 % at the 2 600 m site, compared to the reference site at 3 000 m.

