Climate Change Impact on the Habitat Suitability of Pseudotsuga menziesii Mirb. Franco in Mexico: An Approach for Its Conservation

: One of the conifers that survived after the last glaciation is the Pseudotsuga menziesii (Mirb.) Franco. Due to the gradual increase in temperature, this species was forced to move from the south to the north and to higher elevation, causing a fragmented and intermittent distribution in Mexico. The main objective of this study was to model and identify suitable areas for the future conservation of the P. menziesii in Mexico. The speciﬁc objectives of this research were: (i) to model the habitat suitability of P. menziesii in Mexico, (ii) to identify the most relevant environmental variables based on its current and future habitat suitability (2030, 2050, 2070 and 2090) and (iii) to suggest areas for the conservation of the species in Mexico. Records were compiled from different national and international sources. Climate and topographic variables were used. With MaxEnt software version 3.4.3 (Phillips, New York, NY, USA) 100 distribution models were obtained, where the model showed an area under the curve of 0.905 for training and 0.906 for validation and partial ROC of 1.95 and Z reliable ( p < 0.01), with TSS values > 0.80. The current area of the P. menziesii was 31,580.65 km 2 . The most important variables in the current and future distribution were maximum temperature of the hottest month, precipitation of the coldest trimester and average temperature of the coldest trimester. The percentage of permanence (resilience) for the 2030, 2050, 2070 and 2090 climate horizons was 49.79%, 25.14%, 17.45% and 16.46%, respectively, for the SSP 245 scenario. On the other hand, for the SSP 585 scenario and the analyzed horizons, the percentage resilience in areas of suitable habitat zones was 41.45%, 27.42%, 9.82% and 2.89%.


Introduction
Forests cover approximately 30% of the earth's land surface and provide a wide range of ecological, economic and social services, including the hydrological cycle, carbon capture and beautiful scenery [1,2]. In Mexico, its estimated that 5.4% of the national territory is coniferous forest, 4.9% oak forest and 0.4% mountain mesophyll forest [3]. Among the conifers, one species that survived in Mexico after the last glaciation is Pseudotsuga menziesii (Mirb.) Franco. Its distribution is fragmented and discontinuous. Due to the gradual increase in temperature, this species was forced to move from the south to the north of the country to sites with higher elevation [4].
The genus Pseudotsuga belongs to the Pinaceae, Coniferals order, Pinophyta class and Spermatophyta division. Martínez [5] suggests the existence of four Pseudotsuga species in Mexico: Pseudotsuga macrolepis Flous, P. rehderi Flous, P. flahualti Flous and P. guineri Flous. However, American taxonomists and those of other countries do not accept the existence of these species and only recognize two varieties: Pseudotsuga menziesii var. menziesii, known as fir Douglas, and Pseudotsuga menziesii var. Glauca [6]. The morphological differences of the Mexican populations are a consequence of climatic and ecotypic variation due to topographic and environmental variation where the species is located. Collections of Pseudotsuga in Mexico suggest that the Mexican populations are from the P. menziesii species [7,8].
In Mexico, the species is used for local construction, without abundant commercial production. Its seeds are used to establish commercial Christmas tree plantations. Its exploitation is limited because it is included in the list of special protections in the Mexican standard NOM-ECOL-059 [9]. However, the protections that have been implemented so far have not been sufficient for an adequate conservation of the species because in Mexican areas where populations are reduced and the anthropogenic effect is high, natural repopulation is low, which risks its permanence. Faced with this problem and the effects of climate change, it is important to implement conservation initiatives in the trees' natural habitat, which can be diversified according to the specific conditions of each population. However, for conservation measures to be successful, including the genetic component and the promotion of legal use, it is necessary to know how many trees exist, where they are located and under what conditions they develop.
Currently, the natural distribution of P. menziesii in Mexico includes the north part of eastern and western Sierra Madre with isolated and small groves to the east of the trans-Mexican volcanic belt, and localities in the South and North Sierras of Oaxaca [10,11]. They create forests with altitudes between 1500 and 3600 masl of cold and temperate climates, with deep slopes, and with north, northeast and northwest exposure, according to the greatest humidity amount, with little wind and the least insolation [12]. These climate requirements have not been clearly described, nor has the relevance of each variable in the actual distribution of P. menziesii been weighed.
Species distribution modelling (SDM) is one approach used to model the potential geographical distribution and ecological requirements of a species. This method analyzes those environmental conditions of a species' known occurrence to predict potentially suitable habitats [13]. The potential distribution model has recently become a very useful tool used by researchers to answer research questions and test hypotheses regarding biological, ecological, and evolutionary issues, as well as projections of the impact of climate change at different time scales and ecological niche reconstructions of the past [14][15][16]. Scientific evidence shows that past climate changes have modified the distribution and abundance of terrestrial biodiversity, including Mexican conifer forests [17,18]. Since the last great glaciation, about 21,000 years ago, the planet's climate has changed from very cold periods with a global mean annual temperature up to −8 • C to warm periods, known as interglacial, when the global mean annual temperature was up to 3 • C higher than today [19].
Such changes have been occurring with a greater periodicity due to anthropogenic activity, which is accelerating the rate of such changes [20]. Future climate projections predict an annual temperature increase of up to 2 • C by 2050 and 3.5 • C by 2100, which would become a global threat [21,22]. Faced with such changes in climate variability, contemporary biodiversity would have three possible options: tolerate climate alterations, local and/or regional extinction or change their distribution areas [23,24].
Knowledge of the habitat suitability of forest species helps researchers to understand the environmental patterns that contribute to their establishment, which allows the gathering of special information for conservation management of genetic resources [25]. MaxEnt is a useful tool to predict the effect of climate change on future species distribution [15,18]. Global circulation models are used to simulate future climate scenarios [26]. The combination of current spatial distribution (ecological niche) and global circulation scenarios allows the generation of a probabilistic map of future habitat suitability for a species, thus providing relevant information for its conservation and restoration programs [27]. Indeed, MaxEnt has been successfully used to predict current and future spatial distributions of conifers in Mexico [18,24,28,29].
P. menziesii is a relict species, whose populations were fragmented and isolated in Mexico after the last glacial age and whose environmental requirements have not been analyzed from an ecological and spatial perspective. Moreover, it is not known with certainty how climate change will affect its natural distribution area within the Mexican territory. For this reason, the main objective of this study was to model and identify suitable areas for future conservation of P. menziesii in Mexico. The specific objectives of this research were: (i) to model the habitat suitability of P. menziesii in Mexico, (ii) to identify the most relevant environmental variables according to its current and future habitat suitability (2030, 2050, 2070 and 2090) and (iii) to suggest conservation areas for this species in Mexico.

Materials and Methods
The study area comprised the places where P. menziesii grows in Mexico, with coordinates of 16 •  would become a global threat [21,22]. Faced with such changes in climate variability, contemporary biodiversity would have three possible options: tolerate climate alterations, local and/or regional extinction or change their distribution areas [23,24]. Knowledge of the habitat suitability of forest species helps researchers to understand the environmental patterns that contribute to their establishment, which allows the gathering of special information for conservation management of genetic resources [25]. MaxEnt is a useful tool to predict the effect of climate change on future species distribution [15,18]. Global circulation models are used to simulate future climate scenarios [26]. The combination of current spatial distribution (ecological niche) and global circulation scenarios allows the generation of a probabilistic map of future habitat suitability for a species, thus providing relevant information for its conservation and restoration programs [27]. Indeed, MaxEnt has been successfully used to predict current and future spatial distributions of conifers in Mexico [18,24,28,29].
P. menziesii is a relict species, whose populations were fragmented and isolated in Mexico after the last glacial age and whose environmental requirements have not been analyzed from an ecological and spatial perspective. Moreover, it is not known with certainty how climate change will affect its natural distribution area within the Mexican territory. For this reason, the main objective of this study was to model and identify suitable areas for future conservation of P. menziesii in Mexico. The specific objectives of this research were: (i) to model the habitat suitability of P. menziesii in Mexico, (ii) to identify the most relevant environmental variables according to its current and future habitat suitability (2030, 2050, 2070 and 2090) and (iii) to suggest conservation areas for this species in Mexico.

Materials and Methods
The study area comprised the places where P. menziesii grows in Mexico, with coordinates of 16°45′ to 29°58′ LN and 96°38′ to 108°50′ LW (Figure 1

Geographic Records of Presence
There were 491 P. menziesii records from two information sources; the first is the Global Biodiversity Information Facility [31] and the second is the national forest inventory [32]. To avoid the effect of spatial autocorrelation, and therefore an over estimation of the habitat suitability models, a debugging analysis was developed [33]. This process consists of deleting double records and leaving only one record for 2.5 lineal minutes (~19.87 km 2 ) through the routine of the Niche Toolbox program of the National Biodiversity use and knowledge Committee [34]. After debugging, 126 P. menziesii records were considered for the model process.

Current and Future Climatic Variables
Current climatic information was obtained from the 19 bioclimatic layers of WorldClim version 2.1 (Fick, S.E. and R.J. Hijmans, Environment Institute Sweden, Stockholm, Sweden) ( Table 1), at a special resolution of 2.5 min (~19.87 km 2 ), based on the average climatic information from 1970 to 2000 [30]. For the future habitat suitability analysis, the CNRM-CM6 Global Circulation Model (GCM) [35], which is one of the most recently used in Mexico, was considered [18,26]. The GCM was generated from the Coupled Model Intercomparison

Topographic Variables
The topographic information was obtained from the Mexican digital elevation model 3.0 (DEM) with a resolution of 30m per pixel [36]. The elevation information (ELE) was changed to a spatial resolution of 2.5 min (~19.87 km 2 ) and the slope (PEN) and orientation (ORI) variables were generated with the ArcMap 10.3 software (ESRI, Redlands, CA, USA) [37].

Environmental Variable Selection and Model Calibration
The variable selection that represents the greatest influence on the species habitat was developed through a minimum convex polygon (MCP) with the post-debugging records. This process consisted of placing 10,000 background centroids in order to extract information for each variable [38]. To avoid the multicollinearity effect between variables, a correlation matrix was generated and those with r > 0.70 (p < 0.05) were eliminated [39]. The selected variables were delimited for Mexico-viz., the zone that represents the model area (M), which consists of the geographical space where the species has been registered according to its biological capacity [14]. The model calibration was developed through the standardized coefficient criteria of information Akaike AICc, a process that was done to reduce the model complexity through the PMC variables [40]. The calibration consisted of selecting the lowest ACIc value, which was generated through R version 3.5.3 environment [41], with the ENMeval library [42].

Current Habitat Suitability Model
The habitat suitability model of P. menziesii was generated using the MaxEnt software version 3.4.3 (Phillips, New York, NY, USA) [28]. This software is widely used to determine potential distribution areas of species by occurrence records [43]. The variables used after debugging were maximum temperature of the hottest month, precipitation of the coldest trimester, average temperature of the coldest trimester, average temperature of the most humid trimester, precipitation of the driest month, precipitation of the driest trimester, orientation and elevation. The model parameters were internally replicated by cross validation, 1000 iterations, logistic type exit, 100 replicas and convergence threshold of 0.0001, and 75% of the record was selected for the model and the remaining 25% for validation [44].

Future Habitat Suitability Model and Conservation Zones
Future habitat suitability models were made by transferring the parameter from the current calibration and statistical results of the best model generated in MaxEnt [45].
To determine habitat suitability zones, a suitable-non-suitable threshold was generated through reclassification of the same raster model into three classifications (low, medium and high), and the high classification was used to discriminate zones and transform the continuous models to binary [46] through the ArcMap software version 10.3 (ESRI, Redlands, CA, USA) [37]. Under the climate change scenario, the area occupied by the species and the percentage of change in the current model and the future total model were estimated, as well as each polygon of the Natural Protected Areas in the state and national category for the most critical climate change scenario in Mexico. The areas that could be suitable for species conservation according to their environmental requirements through the intersection between polygons from the current model and the areas of the future model have a 2090 horizon. These areas allow the identification of zones with permissible conditions to generate conservation, reforestation and restoration activities [26].

Validation of Habitat Suitability Models
The current and future habitat suitability models of P. menziesii were evaluated using the area below the curve test (AUC) [47]. The validation test shows a range from 0.0 to 1.0, where values from 0.7 to 0.9 classify the model as good and values above 0.9 classify it as excellent [15]. Due to the questioning of the use of presence data algorithms, by similarly weighting omission and commission errors, the models were additionally evaluated with the partial ROC analysis through the Niche Toolbox platform [34]. The criterion of 1000 replicas by bootstrap with an omission mistake of 5% followed this [48]. To evaluate the performance of the models, z tests (p < 0.01) were generated with a high partial ROC value and a lower standard error. The models were also evaluated by the True Skill Statistics (TSS = sensitivity + specificity − 1) [34].

Contribution of the Environmental Variables
The contribution of the variables used in the model generation were evaluated through the Jackknife test [49]. The proportion of each variable allowed the determination of the relative importance in the development of the P. menziesii habitat suitability for the current and future periods [28].

Model of the Current Habitat Suitability of the P. menziesii
The results of the current suitability models of the 100 replicas of P. menziesii showed the AUC of ROC standard statistical values from 0.900 to 0.905 for training and 0.866 to 0.906 for validation. For the best model, a partial ROC value of 1.95 was presented (S.D. = 0.02), with AUC from 0.905 to 0.906 for training and validation, respectively. The TSS value for the best model was 0.85. The current area of P. menziesii in Mexico was estimated as 31,580.95km 2 (Figure 2). The areas identified as having the highest habitat suitability in Mexico were located mainly in Durango (26.21%), Chihuahua (13.25%), Estado de México (10.75%), Puebla (10.03%), Nuevo León (9.42%) and Veracruz (6.52%). The refuge areas within the polygons of the protected natural areas of P. menziesii were located in the northeast of Mexico, in the areas known as (1)

Contribution of the Environmental Variables
The contribution of the variables used in the model generation were evaluated through the Jackknife test [49]. The proportion of each variable allowed the determination of the relative importance in the development of the P. menziesii habitat suitability for the current and future periods [28].

Model of the Current Habitat Suitability of the P. menziesii
The results of the current suitability models of the 100 replicas of P. menziesii showed the AUC of ROC standard statistical values from 0.900 to 0.905 for training and 0.866 to 0.906 for validation. For the best model, a partial ROC value of 1.95 was presented (S.D. = 0.02), with AUC from 0.905 to 0.906 for training and validation, respectively. The TSS value for the best model was 0.85. The current area of P. menziesii in Mexico was estimated as 31,580.95km 2 (Figure 2). The areas identified as having the highest habitat suitability in Mexico were located mainly in Durango (26.21%), Chihuahua (13.25%), Estado de México (10.75%), Puebla (10.03%), Nuevo León (9.42%) and Veracruz (6.52%). The refuge areas within the polygons of the protected natural areas of P. menziesii were located in the northeast of Mexico, in the areas known as (1) Bajo San Juan River (1383.64 km 2 ) and (2)

Habitat Suitability Model under the Climate Change Scenarios
The maximum and minimum values obtained from the AUC of standard ROC test for training and validation and the CNRM-CM6 model under SSP 245 and 585 scenarios, projected to the 2030, 2050, 2070 and 2090 climate horizons, are shown in Table 2. The partial ROC value for the best model under climate change circumstances is shown in Table 3. The TSS values for all the future models were >0.80.

Habitat Suitability Model under the Climate Change Scenarios
The maximum and minimum values obtained from the AUC of standard ROC test for training and validation and the CNRM-CM6 model under SSP 245 and 585 scenarios, projected to the 2030, 2050, 2070 and 2090 climate horizons, are shown in Table 2. The partial ROC value for the best model under climate change circumstances is shown in Table 3. The TSS values for all the future models were >0.80.

Relevant Variables in the Habitat Suitability of Current and Future Climate
The most important variables for habitat suitability for P. menziesii (Mirb.) Franco identification in Mexico are shown in Table 4.

Future P. menziesii Area in México
The estimated area for P. menziesii habitat suitability under climate change conditions is described in Table 5. The habitat suitability maps for P. menziesii under 2030, 2050, 2070 and 2090 climate horizons under the SSP 245 scenario are found in Figure 3. In addition, maps for the SSP 585 climate scenario and the horizons analyzed are found in Figure 4.
Maximum temperature of the hottest month, Precipitation of the coldest trimester, Precipitation of the driest trimester, Average temperature of the coldest trimester and Precipitation of the driest month 97.2

Future P. menziesii area in México
The estimated area for P. menziesii habitat suitability under climate change conditions is described in Table 5. The habitat suitability maps for P. menziesii under 2030, 2050, 2070 and 2090 climate horizons under the SSP 245 scenario are found in Figure 3. In addition, maps for the SSP 585 climate scenario and the horizons analyzed are found in Figure 4.      The results obtained in the models under the climate change scenario suggest a trend of reduction in the area of P. menziesii in Mexico. The percentage of permanency when considered by the total area (100%)-the one generated by the habitat suitability model with current climate for the climate horizons 2030, 2050, 2070 and 209-was 49.79%, 25.14%, 17.45% and 16.46%, respectively, for SSP 245 scenario. On the other hand, for the SSP 585 scenario and the horizons analyzed, the percentage of resilience in the habitat suitability zones was 41.45%, 27.42%, 9.82% and 2.89%.
The effects of climate change will subject the species to modifications in its current distribution, causing a trend of expansion and contraction in its distribution at different latitudes and altitudes in Mexico, as well as a habitat permanence over time ( Figure 5A). Changes in distribution are predicted, with large contractions in 2030, 2050, 2070 and 2090, for 245 and 585 scenarios, mainly in latitudes higher than 20 • and elevations from 2,000 to 2500 masl ( Figure 5B-I). Given the most critical climate change conditions (2090, SSP 585), it is interesting to see how the species is predicted to expand at latitudes below 17 • with an elevation range of 2400 to 3200 masl ( Figure 5I). Given the most critical climate change outlook, which belongs to the 2090 horizon of the SSP 245 and SSP 585 scenarios, the projected areas that could prevail within the polygons of the federal and state protected natural areas with the largest surface area are described in Table 6.  Given the most critical climate change outlook, which belongs to the 2090 horizon of the SSP 245 and SSP 585 scenarios, the projected areas that could prevail within the polygons of the federal and state protected natural areas with the largest surface area are described in Table 6.

Suggested P. menziesii Conservation Areas in México
According to the results in areas destined for the conservation of P. menziesii, for the 2090 horizon based on the CNRM-CM6 climate model and scenario 245 (conservative scenario), it was possible to identify 4991.96 km 2 with the ideal environmental conditions for the growth and development of the species in the future. Through this analysis, it was possible to identify five important regions located in northern and central Mexico, in the states of Puebla (1201.11 km 2 ), Estado de Mexico (895.61 km 2 ), Veracruz (886.66 km 2 ), Coahuila (569.93 km 2 ) and Nuevo León (518.50 km 2 ) ( Figure 6A). In the most critical scenario, SSP 585 (extreme scenario), only 857.11 km 2 of conservation areas are suitable for the species under the assumption of permanence of an ideal climate in the future. The zones with the largest ideal areas were located mainly in central Mexico in the states of Veracruz (335.47 km 2 ), Puebla (303.97 km 2 ) and Estado de Mexico (151.32 km 2 ) ( Figure 6B).

Spatial P. menziesii Model
The habitat suitability model generated in this study for the current and future periods showed good performance of the AUC of ROC standard statistic in the training (0.905) and validation (0.906) tests. AUC statistical values of 0.7 to 0.9 are considered good model performance, and partial ROC values close to 2 are considered reliable at >95% [15,50]. For Mexico, a previous study by Ramos-Dorantes et al. [51] analyzed the potential distribution of P. menziesii in the state of Puebla, where AUC performance (0.90) was very similar to that obtained in this study.
The total area occupied by P. menziesii in Mexico is currently unknown, as the populations are fragmented and distributed in patches combined and dominated by other species. According to the results of the current habitat suitability model, the presence of the species coincides with different studies that locate it in the northern region of Mexico in the states of Chihuahua, Sonora, Durango, Zacatecas, Nuevo León, Coahuila and Tamau-

Spatial P. menziesii Model
The habitat suitability model generated in this study for the current and future periods showed good performance of the AUC of ROC standard statistic in the training (0.905) and validation (0.906) tests. AUC statistical values of 0.7 to 0.9 are considered good model performance, and partial ROC values close to 2 are considered reliable at >95% [15,50]. For Mexico, a previous study by Ramos-Dorantes et al. [51] analyzed the potential distribution of P. menziesii in the state of Puebla, where AUC performance (0.90) was very similar to that obtained in this study.
The total area occupied by P. menziesii in Mexico is currently unknown, as the populations are fragmented and distributed in patches combined and dominated by other species. According to the results of the current habitat suitability model, the presence of the species coincides with different studies that locate it in the northern region of Mexico in the states of Chihuahua, Sonora, Durango, Zacatecas, Nuevo León, Coahuila and Tamaulipas [7,11,52], as well as the middle regions Hidalgo, Puebla, Tlaxcala and Veracruz [53,54], and more recently the state of Guanajuato [55]. In the southern region of the country, its presence has been reported in the northern and southern Sierra of Oaxaca [11,56].
According to palaeobotanical studies, the species is considered to have a wider distribution, with presence in the Mexico Valley and Chiapas [57], which presumes that at present, relict populations may exist in non-explored areas-an assumption that is strengthened through the results of the current spatial model, where areas with ideal climate and topographic characteristics for the species were located in Mexico State and Chiapas.
The results of the areas suitable for the species under the current climate allow the identification of most of the regions where populations are found in Mexico, but it is important to highlight that a few areas of presence were not identified, as in the case of the state of Querétaro, where the species has been recorded, but even with the good performance of the models of this study, only a reduced area of presence was identified.

Relevant Environmental Variables
The most important variable in this study was the highest temperature of the hottest month, coinciding with different studies related to the habitat suitability of the conifer species, genus Pinus, such as P. lumholtzii, P. teocote, P. pseudostrobus and P. patula in Mexico [58] using only bioclimatic variables. Thus, in studies of the combination of climate, soil chemistry and orographic variables, there is an important contribution of P. hartwegii, P. leiophylla, P. oocarpa and P. rzedowskii species [24]. Therefore, a very high variation in this variable directly affects the ecophysiology of conifers [59].
Research such as that of [16,20,26] predicts that the increase in temperature due to climate change will reduce the distribution areas for tree species of genus Abies, Pinus and Taxodium in Mexico. This study shows that, with a maximum temperature of 22 • C in July, the presence of the species begins to decrease, which is confirmed in the species index cards [60].
Precipitation in the coldest trimester explained 18.4% of the total variability of the habitat suitability model of P. menziesii. In the same way, it was also considered an important variable for habitat suitability models in at least seven Abies species in Mexico, according to the study by Martínez-Méndez et al. [16]. The study of Aceves-Rangel et al. [58] is also relevant for P. lumholtzii at a national level, with an 8.2% lower value compared to that found in this study for P. menziesii, but close to 14.6% found for Taxodium mucronatum [26]. In this analysis, the precipitation of the coldest trimester shows a mean of 125 mm, although P. menziesii is adapted to temperate zones with precipitation up to ±1200 mm throughout the year [60].
Although precipitation in this period may be less than 10% of the annual total [61], low temperatures cause less evaporation. Rainfall in this period is stored in the soil's outline and is used by the species at the beginning of the growth stage, promoting its presence in places where the winter rain is significant [62].
The average temperature in the coldest trimester is a common denominator in the spatial pattern of terrestrial vascular plant species in Mexico [18,58]. This trend was attributed to high temperatures, which promote an increase in evapotranspiration, and consequently a metabolic alteration involving photosynthate assimilation, especially in species of Holarctic origin [63]. Temperature is an important component in global climate change and has shown changes in its variability [64], especially in northern Mexico, where the climate models predict a significant increase [65]. This suggests that growth may occur earlier than usual but with a reduction in the radial increment and biomass production, because the P. menziesii species is susceptible to temperature change [66].
The average temperature in the most humid trimester is a variable with wide participation in the analysis of habitat suitability in conifers such as P. herrerae, P. durangensis, P. engelmani and P. leiophylla [58] and species such as Cedrela odorata in Mexico under current and future conditions [67]. According to the average temperature of the most humid trimester, the average was 16.82 • C, and the average temperature of the coldest quarter is in the optimal range for P. menziesii in Mexico for the month of September. The low temperatures of this season (August-September) and the short photoperiod cause an initial quiescence transition to intense dormancy in P. menziesii [68].
The precipitation of the driest month has been reported as a determinant variable in species with restricted distribution, such as Taxus globosa in Mexico [58] with a participation in the variability of the models of 6.2% [69]. This is similar to what was found in this study (6.9%). Precipitation in the winter season in northern Mexico has a strong influence on the radial growth of P. menziesii [70] because it infiltrates the soil's outline and is used by the deep roots of the species that depend on the stored moisture [71].
In this study, topographic variables (elevation and orientation) represented 2.0% of the variability of the model; however, these types of variables are important in the habitat suitability of forest species in Mexico [72].
The elevation variable indicated a habitat suitability of the species in a range of 2800 to 3200 masl with a value of 0.6%, a range that is consistent with the optimum of 2800 to 3600 masl suitable for its development. Elevation is an important variable for the presence of diverse pine species. It can explain about 30% in the case of Pinus lawsonii [58] and 29.0% for Pinus montezumae [18], as well as for species of the Taxodium and Abies genera in Mexico [16,26].

Climate Change Scenarios
Several studies on climate change scenarios for tree species have been carried out in Mexico. Most focused on cold and temperate climate species, which confirm the theory that there will be an important reduction in the distribution areas for those taxonomies in the future [14,18,27,46]. However, for P. menziesii, there is a lack of spatial information under the climate change context that could help generate strategies in important areas for the conservation of this species in Mexico.
The chosen GCM has been used in the habitat suitability analysis of tree species because it shows the characteristic climate conditions of the geographic area of this study [46]. According to the temperature increases predicted by the CNRM-CM6 model for the horizons and scenarios analyzed, the ecological niche would show a decrease from 16.4 to 2.8% within the natural distribution areas of P. menziesii in Mexico. This is due to the high susceptibility to temperature change of the species [59]. The results obtained by the GCM show similar tends reported by [46], who found ecological niche reductions for the future of P. greggii of 21.8% with the CNRM-CM5 model under climate change scenarios.
The combination of high temperatures and a lower availability of the hydric resource outside of precipitation would be the main causes of the reduction of the habitat of P. menziesii in its natural range in Mexico. However, although there is not a baseline study on the distribution of the species in the country, the pattern of reduction shows a similarity with other forest species with fragmented distribution, such as Taxus globosa [29].
The analysis under climate change scenarios suggests the resilience of P. menziesii in latitudes lower than 21 • N-that is, the creation of a refuge zone in central Mexico. This trend is confirmed by [73], who foresee the disappearance of tree and shrub species in northern Mexico during this century.

Preservation Areas
Studies by [74] show that, out of 70% of the 56 species of the genus Pinus modeled in Mexico, only 10% have their distribution areas protected within a natural protected area. In this study, something similar is shown since it is estimated that between 5 and 28% of Research developed on the conservation of the ecological niche mentions that the species that tend to maintain their niche adapt better to climate change over time or will most likely move to colonize new geographic areas that have similar characteristics to their ecological niche [15].
Ecological niche modelling estimates a probabilistic index of environmental suitability over large-scale areas based on environmental variables of the sites where the species currently grow. Such an index provides detailed information about the niche components because it applies the niche theory as a multidimensional hypervolume proposed by [75], which has been complemented by geographical analysis [14].
According to the relevant variables in the future and current distribution of P. menziesii, it was possible to define the conservation areas of the ecological niche suitable for the species conservation and restoration of the species within its natural range in Mexico. These areas are located mainly in central Mexico, in the states of Puebla, Estado de Mexico and Veracruz and are similar to other conservation proposals made for P. greggii [46] and P. hartwegii [18]. Thus, the environmental and topographic characteristics of the Trans-Mexican volcanic belt would act as climate refuge areas for these species, which could have occupied these niches since the last great glaciation, as occurred with Pinus spp [17] and during the mid-Holocene period (6000 years ago) with Abies religiosa [17]. Nowadays, due to the accelerating processes of modification in the patterns and dynamics of global climate, Mexican conifer species would have to move 100 m in altitude per year to reach a 0.5 • C increase to survive the negative effects of the climate change [21]. Therefore, it is important to support for this type of research that provides different perspectives for decision-making processes in the face of adverse climate scenarios [76] and distinguish between the management and conservation strategies for the country's forest resources.

Conclusions
This study allowed the delimitation and estimation of the current suitability of the natural habitat of Pseudotsuga menziesii (Mirb.) Franco in Mexico within the natural protected areas. Additionally, suitability areas were defined and estimated under two climate change scenarios-SSP 245 and 585-for 2030, 2050, 5070 and 2090 horizons. Currently, there are no studies related to this ecologically and economically important species. The most important climate variable determining ecological niche suitability was the maximum temperature of the hottest month, being a component of interest due to the increase of the variable predicted in the future due to climate change. In this case, a significant reduction in area is predicted in all future projections, but with a contraction tendency towards mid-latitudes, i.e., the disappearance of this species in northern Mexico is suggested. The identification of the presence of P. menziesii in natural protected areas allows researchers to understand the areas of interest to increase the conservation strategies in current climate conditions and the remaining areas that may remain under climate change, specifically in the states of Puebla, Estado de México, Veracruz, Coahuila and Nuevo León.

Conflicts of Interest:
The authors declare no conflict of interest.