Ecological Niches and Suitability Areas of Three Host Pine Species of Bark Beetle Dendroctonus mexicanus Hopkins

Bark beetles are a natural part of coniferous forests. Dendroctonus mexicanus Hopkins is the most widely distributed and most destructive bark beetle in Mexico, colonizing more than 21 pine species. The objectives of this study were to generate ecological niche models for D. mexicanus and three of its most important host species, to evaluate the overlap of climate suitability of the association Dendroctonus–Pinus, and to determine the possible expansion of the bark beetle. We used meticulously cleaned species occurrence records, 15 bioclimatic variables and ‘kuenm’, an R package that uses Maxent as a modeling algorithm. The Dendroctonus–Pinus ecological niches were compared using ordination methods and the kernel density function. We generated 1392 candidate models; not all were statistically significant (α = 0.05). The response type was quadratic; there is a positive correlation between suitability and precipitation, and negative with temperature, the latter determining climatic suitability of the studied species. Indeed, a single variable (Bio 1) contributed 93.9% to the model (Pinus leiophylla Schl. & Cham). The overlap of suitable areas for Dendroctonus–Pinus is 74.95% (P. leiophylla) and on average of 46.66% in ecological niches. It is observed that D. mexicanus begins to expand towards climates not currently occupied by the studied pine species.


Introduction
Temperate forests cover 15% of the world's land [1,2]. In Mexico, they cover 13% the country's territory [3]. Dominated by the genus Pinus, they are especially restricted to mountainous areas, between elevations from 1500 to 4000 m [4], having an affinity for cold temperate climates. These forests contain 49 of the 120 pine species described [5], reaching one of the highest diversities worldwide. Pine forests are of great economic importance for the country. From an ecological perspective, they contribute to regulate global climate acting as CO 2 sinks [6]. However, the accumulation of greenhouse effect gases in the atmosphere, caused mainly by human activity, has caused an increase of 0.87 • C in the last 10 years [7]. This trend is occurring more rapidly than had been predicted and could have significant effects on adaptations of living organisms.
Coexisting naturally with conifer species are species of the genus Dendroctonus Erichson, 1836, called bark beetles [8,9]. The genus includes 19 species and one subspecies, of which 17 are found in North and Central America and only two in Asia and Europe [8].
But only a small percentage of the more than 6000 bark beetles found worldwide are capable of causing significant economic damage [10]. Their name comes from Dendro (tree) and tonus (destroyer) [11]; they belong to the Curculionidae family and play a critical role in the dynamics of coniferous forests [12,13]. They vary in size from 0.1 to 6 mm and are endophytes; that is, they dig galleries under the bark of live trees to feed on the phloem [1].

Study Area
The study area differs for each of the species; in general, there are coniferous forests located in the main mountain systems of Mexico, such as the Sierra Madre Occidental (SMOc), Sierra Madre Oriental (SMOr), Trans-Mexico Volcanic Belt (TMVB), Sierra Madre del Sur, Sierra del Norte de Oaxaca, Sierras de Chiapas, and the extremes of Baja California. In these regions annual mean temperatures oscillate between 10 and 20 • C, and annual precipitation is 600 to 1000 mm [4]; altitudes are from 1600 m to a little more than 3000 m ( Figure 1).

Materials and Methods
The conifer species P. leiophylla, P. teocote and P. devoniana were selected because they have the highest percentual incidence of bark beetle D. mexicanus attack in Mexico and are the most susceptible species, with 35.6, 13.9 and 9.4%, respectively [14].

Study Area
The study area differs for each of the species; in general, there are coniferous forests located in the main mountain systems of Mexico, such as the Sierra Madre Occidental (SMOc), Sierra Madre Oriental (SMOr), Trans-Mexico Volcanic Belt (TMVB), Sierra Madre del Sur, Sierra del Norte de Oaxaca, Sierras de Chiapas, and the extremes of Baja California. In these regions annual mean temperatures oscillate between 10 and 20 °C, and annual precipitation is 600 to 1000 mm [4]; altitudes are from 1600 m to a little more than 3000 m (Figure 1).

Bioclimatic Variables and Selection
Because of the scale of the study area (>200 km), only bioclimatic variables were used (Table 1) [23,32]; these proposed by [33] were re-sampled at a resolution of ∼5 km 2 . Bio 8, Bio 9, Bio 18 and Bio 19 were excluded from the analysis since, by combining information on precipitation and temperature in the same layer, the resulting predictions were erratic and biased [34].
The selection of the bioclimatic variables was based on four criteria: (1) relative contribution of the variable to the bioclimatic profile of the species [35] was obtained through a principal components analysis (PCA) with the package 'FactoMiner' [36], previously extracting from each species occurrence records the value of the 15 bioclimatic variables, making the PCA to the standardized variables and selecting those that contributed most, (2) non-correlated variables (r < 0.8) were determined by a parametric correlation analysis (α = 0.05) of the variables transformed to natural logarithm [25,37,38], (3) variable frequency distribution using the Sturges [39] rule determined how the bioclimatic variable was distributed giving selection priority to those close to a normal distribution or a skewed distribution (left or right) [40], and (4) the predictive capacity of the variable consists of preliminary modeling with individual variables and transferring the model in time and space because it has been demonstrated that climate projections of precipitation of the General Circulation Models (GCMs) are biased at higher latitudes [33], resulting in an overestimation of climate suitability of a species in the same zones. The variables that least overestimate climate suitability of a species were selected. These procedures were carried out in each species.  14.87 (2) 9.82 (2) 12.36 (1,2,3) 20.08 (2) Bio 7 Annual Temperature Range ¶ ( • C) 8. 18 14.33 (1,3) 14.86 (1,2,3) 9.24 (1,2,3) Bio 10 Mean Temperature of Warmest Quarter ¶ ( • C) 9.85 (1,2,3) 9.83 (2) 14.42 (1,2) 11.40

Species Occurrence Records and Cleaning
Species occurrence records of each species were downloaded from the Global Biodiversity Information Facility (GBIF) and Red Mundial de Información sobre Biodiversidad (REMIB). Others were obtained by Inventario Nacional Forestal y de Suelos de México (INFyS), published articles and our own fieldwork.
Data cleaning of each specie consisted of eliminating records that were: (1) Outside the geographic range (latitude and longitude), (2) outside the altitudinal range (according to the descriptor of the species), allowing records between quantile 5 and 95 (and/or ±500 m), depending on the species, (3) not precise (equal to or less than three digits), (4) duplicates [17], (5) without author of identification, and (6) outside the ellipse at 99% of a PCA conducted with the package 'FactoMiner' [36], using 15 environmental variables and altitude. Spatial autocorrelation between records was then eliminated with the package 'spThin' [41], allowing a single record per pixel (∼5 km). If we did not correct for this characteristic, we would have incurred in a biased selection of variables or model coefficients [42].

Calibration Area
Areas where the species could be observed or that could be explored because of their biological capacity are denoted as «M» in the «BAM» diagram [29]. This was delimited preliminarily in ArcMap v.10.5, applying a 70 km buffer radius to each species occurrence records. In SDM and ENM, the calibration areas should include the complete distribution of the species [20]. If the extension of this area is smaller than the species range, the response function cannot have the predicted form, because of the niche theory [22]. The final delimitation of «M» was achieved with the cleansed records, that is, with those that passed all the cleaning criteria indicated in the previous section.

Model Calibration, Creation, and Evaluation
The model calibration, creation, and evaluation were done in 'kuenm', an R package that uses Maxent (maximum entropy) [43] as the modeling algorithm. Maxent has two main modifiable parameters: (1) "regularization multiplier" (β) and (2) "feature classes" such as linear (l), quadratic (q), product (p), threshold (t) and hinge (h). The first is a parameter that adds new constraints; it is a penalty imposed on the model, and the latter corresponds to a mathematical transformation of the different covariates used in the model to allow complex relationships to be modeled [37,44]. For each species, 16 regularization multiplicators were tested (0.1 to 1, 2 to 6, and 10), 29 response types (l, q, p, t, h, lq, lp, lt, lh, qp, qt, qh, pt, ph, th, lqp, lqt, lqh, lpt, lph, qpt, qph, qth, qth, lqpt, lqph, lqth, lpth and lqpth), and three different sets of environmental variables (optional), those that satisfied the selection criteria indicated in Section 2.2.
Modeling was carried out with approximately 70% of the records; with independent data (∼30%), the predictive capacity of the models was evaluated using cross validation [36]. This percentage depends on the availability of independent records of each species used for validation. The output format was of logistical type, which can be interpreted as probability of presence and is recommended if and only if the "regularization multiplier" and "feature classes" [45] are optimized. The resulting models (maps in raster format) represent species suitability values (0-1) [46]. The best fit model was selected according to: (1) the statistic partial ROC (Receptor Operated Curve) [47], (2) the rate of omission <0.05%, (3) the lowest value of the Akaike Information Criterion (AICc) [19,48,49], (4) species response curves to the environmental gradients [37], and (5) statistical significance of the model, p-values [19]. Here, the partial ROC was used instead of the area under the ROC curve (AUC) because the latter is not a good measure of fit in ENM [47,50]. Statistical significance was determined by a bootstrap resampling of 50% of testing data.
To reduce bias, the final selected model was represented by the mean of 10 repetitions, with which prediction uncertainty was obtained. The jackknife tests and the response curves to bioclimatic variables were implemented in 'kuenm' [43] to determine the variable's contribution to the model. The process for generating the climate suitability model of the species is shown in Figure 2.

Model Stratification
The final continuous model of probability (suitability) of each species was classified in three levels or strata: low, medium and high suitability. To this end, 5000 points were distributed randomly over the suitability model; from these points their value was extracted. Later, in R with the package 'stratifyR' [51], the thresholds of each stratum were calculated following the method of Khan [54]. This method determines the Optimum Strata Boundaries (OSB) and Optimum Sample Sizes (OSS) for the study variable, using the best-fit frequency distribution of a survey variable. It formulates the problem of determining the OSB as a mathematical programming problem, which is solved using a dynamic programming technique, using Neyman Allocation assuring minimum intra-stratum variance and maximum inter-stratum variance.

Bark Beetle-Free Areas
The final suitability models of each species were converted into binary maps in ArcMap v.10.5 to represent climate suitability-non-suitability. This was done by reclassifying the suitability to 1 and 0; the value of one (1) was assigned to suitability comprehended between the minimum value of the second stratum (calculated in the pre-

Model Stratification
The final continuous model of probability (suitability) of each species was classified in three levels or strata: low, medium and high suitability. To this end, 5000 points were distributed randomly over the suitability model; from these points their value was extracted. Later, in R with the package 'stratifyR' [51], the thresholds of each stratum were calculated following the method of Khan et al. (2002) [52], Khan et al. (2008) [53] and Khan et al. (2015) [54]. This method determines the Optimum Strata Boundaries (OSB) and Optimum Sample Sizes (OSS) for the study variable, using the best-fit frequency distribution of a survey variable. It formulates the problem of determining the OSB as a mathematical programming problem, which is solved using a dynamic programming technique, using Neyman Allocation assuring minimum intra-stratum variance and maximum inter-stratum variance.

Bark Beetle-Free Areas
The final suitability models of each species were converted into binary maps in Ar-cMap v.10.5 to represent climate suitability-non-suitability. This was done by reclassifying the suitability to 1 and 0; the value of one (1) was assigned to suitability comprehended between the minimum value of the second stratum (calculated in the previous section) and that of maximum suitability, while the value of zero (0) corresponded to the rest of the suitability. The binary maps were manipulated using 'raster algebra' of each pair of species (Pinus-D. mexicanus), and SAP free of SAD were calculated. The procedures based on spatial predictions of climate suitability (also called ENMs), besides estimating SAP-SAD [49], allow quantifying changes and niche overlap of two or more species in G space [20].

Quantifying Niche Similarity
Similarity between the Dendroctonus and Pinus niches was calculated with two indexes introduced by Warren [55] D and a measure derived from Hellinger's distance called I using ordination methods (PCA) and the 15 bioclimatic variables for each species [20] (Table 1). This method is the most precise [56] and the most recommended [20,57]. Both similarity measures range from 0, when species predicted environmental tolerances do not overlap at all, to 1, when overlap is total. Niche comparisons were made relative to the entire niche of the species, pooled from the two ranges [20,31], under the hypothesis of niche conservatism of the genus Pinus; it is known that it was established ∼145 million years ago in the lower Cretaceous [58], and the bark beetle as an invasive species. A kernel density function (standard smoothing parameters) was applied to determine the 'smoothed' density of occurrences in each cell in the environmental space for each dataset; the use of a kernel smoother makes the process of moving from G space to multivariate E space, independent of both sampling and resolution in environmental space. The shift of the middle position from centroid of the niche between the species of bark beetle with the pine species was calculated with the C metric, using Euclidian distance [20]. This analysis was performed with the package 'ecospat' [59]. All the packages mentioned in this study were run in R 3.6.3 [60].

Generalities
A total of 283, 3648, 2209 and 772 species occurrence records were collected for D. mexicanus, P. leiophylla, P. teocote and P. devoniana, subtracting for the modeling and validation 86, 900, 736 and 255 (30.39, 24.67, 33.32 and 33.03%); that is, up to 75.33% of the P. leiophylla records were eliminated, because they did not meet the established cleaning criteria (Section 2.3).
The PCA performed to select the variables by their contribution, explained 67.63, 73.48, 65.60 and 64.10% for the different species (Table 1). This analysis is statistically valid when the variables correlate with each other (>0.6). The Kaiser-Meyer-Olkin (KMO) index indicated that the global correlation of the PCA was 0.68, 0.71, 0.72, and 0.65, respectively. The results show that in the four species, the variables derived from temperature (Bio 1-Bio 11) are those that most contributed; variable 16 (Precipitation of Wettest Quarter, mm) is that which least contributed (Table 1). There are variables (e.g., Bio 12 and Bio 13, for the species D. mexicanus) that contribute significantly in the PCA (Table 1). However, they were not selected because they presented a bimodal distribution.

Generated Models and Their Statistics
A total of 1392 candidate models were generated for each species (Table 2). It is notable that not all the models were statistically significant (α = 0.05), registering 99.4% for D. mexicanus, 53.5% for P. leiophylla, 99.9% for P. teocote and only 9.5% for P. devoniana (Table 2), using the default parameters in Maxent does not necessarily produce the best model [61] (see Appendix A). The type of response that prevailed in the selected models was quadratic (Table 2). No linear response model was selected in any case. The algorithm implemented in 'kuenm' [43] selected the best regulating multiplier in each species varied from 2 (D. mexicanus and P. devoniana) to 5 (P. leiophylla). Fewer than 20% of the 1392 models generated passed the omission rate (false negative predictions) established (0.05%); indeed, none satisfied in the species P. leiophylla, having to increase the established value to 0.07% to select the model. Only one model (0.07% of the total) passed all the criteria established in Section 2.5 (Table 2) in D. mexicanus and P. devoniana.

Species Suitability Areas
In no case did climate suitability reach the maximum value (1), varying from 0.01 (P. teocote) to 0.82 (P. leiophylla). Suitability takes on different forms, triangular in D. mexicanus and P. devoniana with parameters a = 0.001 and 0.001; b = 0.732 and 0.795; c = 0.518 and 0.01, respectively. Gamma distribution in P. leiophylla and P. teocote had values of k = 0.702 and 0.526; λ = 3.082 and 3.728, respectively. The OSBs varied in each species and not in the same amplitude, showing that variance differs in all the spectrum of species climate suitability. For this reason, stratification should obey a statistical technique rather than stratify by proportions, assuring minimum and maximum intra-stratum variance among them.
It is estimated that the area of high suitability in Mexico is 234,649.1, 212,497.4, 177,904.8 and 159,630.4 km 2 for D. mexicanus, P. leiophylla, P. teocote and P. devoniana (Figure 3a-d). Except for P. teocote (Figure 3c), high suitability represents that largest part of «M», from 39.89% (P. devoniana, Figure 3d) to 47.85% (P. leiophylla, Figure 3b). The predicted area of suitability for a species is not dependent on the number of registers. Maxent uses presence records to model climate suitability. During this process, the algorithm generates "pseudo absences" where the species is not present. This helps to improve the predictions of the current distribution of the species.

Bioclimatic Profile
The 'kuenm' algorithm selected the set of variables ( Table 3) that showed the best predictive capacity, validated with the set of independent records of the species. The sets comprised different numbers of variables, from three (P. leiophylla) to seven (P. teocote). Like the PCA, the jackknife tests showed that the variables derived from temperature (Bio 1-Bio 11) contribute more than 80% to explain the bioclimatic profile of the species (Table 3), especially, the representatives of extreme values. In D. mexicanus and P. leiophylla a single variable contributes significantly to the species bioclimatic profile, 87.8% (Bio 10) and 93.9% (Bio 1), respectively. The average coefficient of variation of the variables that most contribute to the bioclimatic profile of each species (Bio 10, Bio 1, Bio 10 and Bio 11, Table 3) is 14.3%, indicating that the predictors selected by both PCA (preselection of variables) and the 'kuenm' algorithm adequately represent the species bioclimatic profile. Only one variable (Bio 6, in P. teocote) experienced high variability (198.9%) but contributes to the bioclimatic profile with only 15.6% (Table 3).

Bioclimatic Profile
The 'kuenm' algorithm selected the set of variables ( Table 3) that showed the best predictive capacity, validated with the set of independent records of the species. The sets comprised different numbers of variables, from three (P. leiophylla) to seven (P. teocote). Like the PCA, the jackknife tests showed that the variables derived from temperature (Bio 1-Bio 11) contribute more than 80% to explain the bioclimatic profile of the species (Table 3), especially, the representatives of extreme values. In D. mexicanus and P. leiophylla a single variable contributes significantly to the species bioclimatic profile, 87.8% (Bio 10) and 93.9% (Bio 1), respectively. The average coefficient of variation of the variables that most contribute to the bioclimatic profile of each species (Bio 10, Bio 1, Bio 10 and Bio 11, Table 3) is 14.3%, indicating that the predictors selected by both PCA (preselection of variables) and the 'kuenm' algorithm adequately represent the species bioclimatic profile. Only one variable (Bio 6, in P. teocote) experienced high variability (198.9%) but contributes to the bioclimatic profile with only 15.6% (Table 3).

Climatic Suitability of Pine Species, Free of Bark Beetle Suitability Areas
From the total suitable area predicted for P. leiophylla (444,100.4 km 2 ), P. teocote (729,358.3 km 2 ), and P. devoniana (400,142.8 km 2 ) (Figure 3b-d) and by obtaining SAP-free of SAD of each Pinus-Dendroctonus pair (Section 2.7), we found that only 92,995.2, 11,737.4 and 55,964.8 km 2 , respectively, are free of suitable bark beetle areas. For P. teocote, of the total of the suitable areas, only 3.02% is left.
Despite the wide distribution of P. leiophylla and P. teocote (Figure 3b-c), SAP-free of SAD are observed only in a northern part of the SMOc (Figure 4a,b), where there is a larger number (seven) of bark beetle species [13], in a compact form for the first pine species P. leiophylla and disperse for the second P. teocote, but inexistent in lower latitudes of the species distribution. In P. devoniana (Figure 4c), the SAP-free of SAD are observed discontinuously over the entire area of its distribution. In all cases, these areas are observed where there is high climatic suitability for the pine species and the bark beetle. Chihuahua and Sonora have 81.56% of the total SAD-free P. leiophylla suitable areas. Chihuahua and Durango contain 57.60% of the SAD-free P. teocote suitable areas, while Jalisco contributes 30.39% of the SAD-free P. devoniana suitable areas. It is possible that, in the future, these areas (SAP-free of SAD) may be susceptible to the bark beetle. MeanCI = confidence interval of mean, 0.05, 0.95 = quantiles of the bioclimatic variable, SD = standard deviation, CV = coefficient of variation (%); MAD = median absolute deviation, IQR = interquartile range.

Climatic Suitability of Pine Species, Free of Bark Beetle Suitability Areas
From the total suitable area predicted for P. leiophylla (444,100.4 km 2 ), P. teocote (729,358.3 km 2 ), and P. devoniana (400,142.8 km 2 ) (Figure 3b-d) and by obtaining SAPfree of SAD of each Pinus-Dendroctonus pair (Section 2.7), we found that only 92,995.2, 11,737.4 and 55,964.8 km 2 , respectively, are free of suitable bark beetle areas. For P. teocote, of the total of the suitable areas, only 3.02% is left.
Despite the wide distribution of P. leiophylla and P. teocote (Figure 3b-c), SAP-free of SAD are observed only in a northern part of the SMOc (Figure 4a,b), where there is a larger number (seven) of bark beetle species [13], in a compact form for the first pine species P. leiophylla and disperse for the second P. teocote, but inexistent in lower latitudes of the species distribution. In P. devoniana (Figure 4c), the SAP-free of SAD are observed discontinuously over the entire area of its distribution. In all cases, these areas are observed where there is high climatic suitability for the pine species and the bark beetle. Chihuahua and Sonora have 81.56% of the total SAD-free P. leiophylla suitable areas. Chihuahua and Durango contain 57.60% of the SAD-free P. teocote suitable areas, while Jalisco contributes 30.39% of the SAD-free P. devoniana suitable areas. It is possible that, in the future, these areas (SAP-free of SAD) may be susceptible to the bark beetle.
According to our results, uncertainty is not more than 30% of the coefficient of variation (Figure 4d-f). The lowest uncertainty (<0.15%) occurs in the areas of high suitability and the highest (>97%) in areas of low suitability.

Overlap of Suitability and Ecological Niches
The overlap of suitable areas of the bark beetle with those of P. leiophylla, P. teocote and P. devoniana in space G is 74.35, 96.98 and 82.44%, respectively (Figure 5a-c). The environmental space (Figure 5d-f) was built with Bio 5, Bio 6, and Bio 12 for the pur- According to our results, uncertainty is not more than 30% of the coefficient of variation (Figure 4d-f). The lowest uncertainty (<0.15%) occurs in the areas of high suitability and the highest (>97%) in areas of low suitability.

Overlap of Suitability and Ecological Niches
The overlap of suitable areas of the bark beetle with those of P. leiophylla, P. teocote and P. devoniana in space G is 74.35, 96.98 and 82.44%, respectively (Figure 5a-c). The environmental space (Figure 5d-f) was built with Bio 5, Bio 6, and Bio 12 for the purpose of comparing the fundamental niche and overlap between bark beetle and pine species. In Bio 5 (maximum temperature), the four species have the same tolerances, from 1 to 32 • C (Figure 5d-f); in Bio 6 (minimum temperature) the overlap is -2 to 11 • C, but D. mexicanus has broader tolerance (from −1 to 17 • C). The lowest occurs in P. leiophylla (from −4 to 11 • C). In Bio 12 (annual precipitation), the overlap occurs between 450 and 1755 mm. P. devoniana, the species of limited distribution, possesses the widest interval, from 429 to 2053 mm (Figure 4f); D. mexicanus tolerates a smaller interval (450 to 1755 mm). The fundamental niche of the pine species shows the same disposition in tridimensional space, but different from that of bark beetle species (Figure 5d-f), but in all niches unavailable climate (UC) and the existence of a potential niche (PN) are observed (Figure 5d-f).  (Figure 6a), P. teocote ( Figure 6b) and P. devoniana (Figure 6c). The variance explained by the first two principal components was 61.76% (P. devoniana) to 71.73% (P. leiophylla). In all cases, the variables derived from temperature contributed more in the PCA (PC1) and those of precipitation (PC2) contributed less (Figure 6d-f). Niche similarity between the 'invasive' species (D. mexicanus) and the pine species resulted in D = 0.48 and I = 0.67; D = 0.39 and I = 0.61; D = 0.53 and I = 0.69, for P. leiophylla (Figure 6a), P. teocote ( Figure 6b) and P. devoniana (Figure 6c). The variance explained by the first two principal components was 61.76% (P. devoniana) to 71.73% (P. leiophylla). In all cases, the variables derived from temperature contributed more in the PCA (PC1) and those of precipitation (PC2) contributed less (Figure 6d-f).
99.15, and 95.37%, respectively. The shift from niche centroid (C) of the 'invasive' species was more significant with P. leiophylla and P. teocote (Figure 6a-b), but in different directions, and almost the same centroid as observed with P. devoniana (Figure 6c). In the three cases, the shift moved over the temperature gradient. Under the similarity hypothesis of ecological niches of the species in niche conservatism (Pinus) and 'invasive' species (Dendroctonus), it is observed that both measures of niche similarity are significantly higher than what was expected of this null distribution, with p < 0.05 (Figure 6g-i). Therefore, this hypothesis is rejected, except in the case of D. mexicanus with P. devoniana (Figure 6g-i) where p > 0.05. Figure 6. Niche of the species along the two first axes of the PCA (a-c) of the native species and 'invasive' species. Green (Pinus) and blue (Dendroctonus mexicanus) (D), shading shows the density of the occurrences of the species by cell. The solid and dashed contour lines illustrate, respectively, 100% and 95% of the available (background) environment, E represents niche expansion and Up the non-overlapping niche between pine species and the 'invasive' species. The arrows represent the change of the niche centroid of the 'invasive' species, relative to the 'native' species. The contribution of the climatic variables on the two axes of the PCA and the percentage of inertia explained by the two axes (d-f). Histograms show the observed niche similarity I between the two ranges (lines with a diamond) and simulated niche similarity (grey bars). Pinus leiophylla (a,d,g), Pinus teocote (b,e,h) and Pinus devoniana (c,f,i). Figure 6. Niche of the species along the two first axes of the PCA (a-c) of the native species and 'invasive' species. Green (Pinus) and blue (Dendroctonus mexicanus) (D), shading shows the density of the occurrences of the species by cell. The solid and dashed contour lines illustrate, respectively, 100% and 95% of the available (background) environment, E represents niche expansion and Up the non-overlapping niche between pine species and the 'invasive' species. The arrows represent the change of the niche centroid of the 'invasive' species, relative to the 'native' species. The contribution of the climatic variables on the two axes of the PCA and the percentage of inertia explained by the two axes (d-f). Histograms show the observed niche similarity I between the two ranges (lines with a diamond) and simulated niche similarity (grey bars). Pinus leiophylla (a,d,g), Pinus teocote (b,e,h) and Pinus devoniana (c,f,i).
The proportion of the native niche (P. leiophylla, P. teocote and P. devoniana) nonoverlapping with the 'invasive' niche (D. mexicanus), denoted as Up (Unfilling) is 31.57, 19.44, and 21.15%; niche expansion E is 1.31, 0.84, and 4.62, while niche stability is 98. 68, 99.15, and 95.37%, respectively. The shift from niche centroid (C) of the 'invasive' species was more significant with P. leiophylla and P. teocote (Figure 6a-b), but in different directions, and almost the same centroid as observed with P. devoniana (Figure 6c). In the three cases, the shift moved over the temperature gradient. Under the similarity hypothesis of ecological niches of the species in niche conservatism (Pinus) and 'invasive' species (Dendroctonus), it is observed that both measures of niche similarity are significantly higher than what was expected of this null distribution, with p < 0.05 (Figure 6g-i). Therefore, this hypothesis is rejected, except in the case of D. mexicanus with P. devoniana (Figure 6g-i) where p > 0.05.

Species Occurrence Records in ENM
The used of reliable records is fundamental in ENM to avoid bias in prediction, mainly because the primary source of data is an opportunist sampling [22]. In general, 70% of the records were eliminated. The use of altitudinal limits of species distribution in records cleansing has not been documented. Here, it can be observed that altitude (included in the PCA with the Bios) was crucial to identifying atypical and erroneous records of the species, especially within «M». Much has been discussed over the number of species occurrence records to use in ENM; some show that 30 should be a minimum [25], while others indicate that 50 is sufficient [62]. Modeling D. mexicanus in this study was carried out with 86 records. Modeling the pine species surpassed 250 records. Actually, the number of observations is less important than adequately representing species prevalence distributed in the entire geographic and environmental space it occupies [22], for which it is of utmost importance that in SDM and ENM the records comprehend the complete distribution of the species (as was done here). Otherwise, the response of the variables would be wrong [20], as would be the ecological niche. It has been demonstrated that systematic sampling produces greater precision in species distribution models [22]. Some records used here (INFyS) come from this type of sampling, assuring more robust predictions. Moreover, this requirement was corrected by removing spatial autocorrelation.

Variable Importance in ENM
After Busby (1986) [26] defined bioclimatic variables from precipitation and temperature, a number of different predictors have been generated, but used indistinctly in ENM, despite the fact that some authors [22,29,32] have suggested their use depending on the study scale. Modeling the species of this study was carried out only with bioclimatic variables, [22,29,32]. The rigorous selection of predictors and the best configuration of the algorithm made it possible to choose a robust model to predict the climatic suitability of each species. Other researchers [19] demonstrated that if this process is appropriate the model will show minimum spatial correlation in their residuals. The use of dynamic (bioclimatic) variables for modeling allows determining the vulnerability of a species, in an unstable environment which does not occur with static variables (soil, slope, exposure and altitude) [63]. In 2017, Fourcade and other researchers surprised the scientific community when they demonstrated that pseudo-predictors derived from paints ('classical paintings') downloaded from Google Image ® , predicted species distribution even better than bioclimatic variables. Undoubtedly, such pseudo-predictors will not satisfy the criteria established in this study (Section 2.2), and so it is crucial to select the environmental predictors in ENM. To generate our models, besides considering these criteria, the 'kuenm' algorithm selected the best set of variables based on the validation of 1392 candidate models, using the set of independent records.
It has been demonstrated that using multiple variables brings problems of bias and uncertainty in the predictions [38,42] and a decrease in statistical power; complex models with a large number of parameters tend to overestimate the predictions [38,64] but this depends on the combination of the adjustment parameters of the "regularization multiplier" and "feature classes" [43]. However, this occur especially when the number of records is small [38,64]. Our models included no more than seven predictors avoiding these problems. Other authors argue that a large number of predictors resulted in problems of collinearity [25,35,38], which tend to inflate the variance of response variables and parameters [42], and an erroneous representation of species distribution [16,40], moreover, biological interpretation of the model is complex or null. Austin (2002) [65] argues that species responses are often nonlinear; likewise, ecological theory suggests response curves are often unimodal [66], and hence, quadratic features may be appropriate. None of our generated models for the studied species was linear; most were quadratic (Table 2), in accord with theory these authors, which demonstrates the correct selection of predictors, quality of records and modeling. Studies reveal that pine species show non-linear responses to climate suitability [67], like what occurs in bark beetle species [18].
Our findings indicate that regardless of the genus and species, the variables of temperature are those that most contributed to the bioclimatic profile, especially seasonal temperature means and extremes (e.g., Bio 10 and Bio 11, Table 3), as occurs in most studies of ENM [67][68][69], even when topographic [28], soil and vegetation [67] variables are included in modeling. For bark beetles, the variables that best predict climatic suitability are temperature [12,70]; e.g., Bio 1 in Dendroctonus rizophagus Thomas & Bright [35]; Bio 5 in Dendroctonus valens Le Conte [71] and Bio 10 in D. mexicanus [18], and in this last species Bio 7, Bio 8, and Bio 10 [28]. It has been found that D. mexicanus possesses a more ample bioclimatic profile than that of several other bark beetle species [18], making it possible to occupy larger geographic areas and probably adapt to new host pine species.
It is surprising that a single variable, e.g., Bio 1 in P. leiophylla and Bio 10 in D. mexicanus, contributes more than 90 and 87%, respectively (Table 3). Like what we found here, it has been demonstrated [71,72] that a single variable predicts plant species distribution very well. In contrast, other authors [23] report models of up to 38 predictors. The number of variables that make up a model (as long as they have been correctly selected) possess a fundamental interpretation. Our analysis shows that this could determine species vulnerability. When the contribution of a single variable is high, there is a risk that if it changes (increases/decreases) and varies over time, it will have important effects on predictions of climate suitability, making the species vulnerable in the same proportions. In contrast, if the model is composed of multiple variables, their contribution would be shared, increasing the possibility that not all are being modified at the same rhythm, even if they do not.
Specifically, in the study area, Bio 1 has increased significantly in recent years [7], a fact that we have verified through climate projections of the GCM's [33], with an expected increase of 2 to 3.5 • C by the year 2050, making P. leiophylla and D. mexicanus highly vulnerable to climate change. The latter is even more vulnerable because of its high dependence and sensitivity to temperature [10], while P. teocote and P. devoniana are less vulnerable. Previous studies on ecological niche models [42] demonstrate that P. leiophylla is highly vulnerable to climate change, as found here.

Climate Suitability and Niche Overlap
The stratification of species climate suitability can determine priority areas for conservation or management, especially those of high suitability. The areas of low suitability represent areas that in the past were of medium of high suitability, but nothing can be done to remedy this condition. Low suitability for D. mexicanus ( Figure 3a) and P. devoniana (Figure 3d) is observed only at higher latitudes and almost never at lower latitudes, challenging distributional ecology where maximum suitability of the species occurs toward the center of spaces E and G [24]. Analogously, it has been found [13] that greater suitability of D. mexicanus is registered in small portions of the TMVB and the SMOc at higher altitudes and where Pinus-Dendroctonus host diversity is high but discontinuous [14,49]. For Pinus, suitability also occurs at higher altitudes [42], possibly following regimes of higher moisture and lower temperatures.
Superimposing climate suitability (Dendroctonus-Pinus) resulted in very few discontinuous SAP free of SAD (between 30.03 and 25.65%), fortunately, where suitability for pine species is high, in contrast with [49], where less than 1% of all the Pinus host species is D. rhizophagus-free, these differences are explained by the cutoff threshold selected to obtain the overlap. It has been demonstrated that the analysis of suitability overlap in G is problematic because it will depend on the extension and distribution of the environmental gradients in the study area [56]. Despite the importance of identifying SAP-free of SAD, only the study of Smith et al. (2013) [49] is known.
The overlap (D) of the 'invasive' species niche with the 'conservationism' species with the highest percentage of incidence (P. leiophylla) is almost 50%, but a change in niche of D. mexicanus, relative to the three pine species has been observed (Figure 6a-c). The niche overlap in similar taxonomic groups of pine species is not very high, on average D = 0.20 [67]. Studies reveal that one of the two species (with high plant prevalence) shows a niche shift [20]; this number is quite high, especially in exotic species, it is possible that the studies reported a shift when there was actually none. Less than 1% of the studies show niche conservationism. Some authors [73] suggest that dramatic niche changes found should be carefully interpreted since they are dependent on the methods and data used. The kernel density method in ecological niche studies (used here) is one of the most robust and produces the best results [56,57].

Conclusions
Not all the models generated in Maxent were statistically significant (α = 0.05). In 'kuenm' it is possible to generate n candidate models and to select a robust model rather than a random model. The statistical procedure (PCA) was a crucial tool for preselecting climate predictors. By including altitude in the analysis (PCA), it was possible to identify atypical and erroneous records within the calibration area, which is not possible through bivariate environmental space. The variables representing extreme temperatures play an important role in defining species climate suitability; they are also indicators of climate change and thus evidence that this will affect species distribution, proportional to their rate of change. The areas of overlapping climate suitability in geographic space and of niches in environmental space average 84.59 and 46.66%, which indicates small areas of Pinus species free of the bark beetle that are isolated in the distribution area. The ordination methods show that pine species have a broader ecological niche, but P. leiophylla and D. mexicanus were identified as highly vulnerable to climate change. In addition, expansion of the bark beetle toward new climates is observed and, consequently, toward new geographic areas following its climate preferences. The areas of high suitability for Pinus, especially those areas free of suitability areas for D. mexicanus, should be prioritized for conservation. The redistribution of the bark beetle species is highly probable in the coming years, consequently, fewer suitable areas for the pine species, free of bark beetle. It is proposed that, when generating ecological niche models, robust methodologies be used, considering the association Dendroctonus-Pinus. This study enriches previous knowledge of the species, improving the delineation of geographic distribution of their ecological niches and specific climate tolerances, contributing tools for the timely implementation of actions and strategies for managing the country's forests for species conservation and preservation.