Abundance of Dendroctonus frontalis and D. mexicanus (Coleoptera: Scolytinae) along altitudinal transects in Mexico: Implications of climatic change for forest conservation

Bark beetle infestations have historically been primary drivers of stand thinning in Mexican pine forests. However, bark beetle impacts have become increasingly extensive and intense, apparently associated with climate change. Our objective was to describe the possible association between abundance of bark beetle flying populations and the occurrence of given value intervals of temperature, precipitation and their balance, in order to have a better comprehension of the climatic space that might trigger larger insect abundances, an issue relevant in the context of the ongoing climatic change. Here, we monitored the abundance of two of the most important bark beetle species in Mexico, Dendroctonus frontalis and D. mexicanus. We sampled 147 sites using pheromone-baited funnel traps along 24 altitudinal transects in 11 Mexican states, from northwestern Chihuahua to southeastern Chiapas, from 2015 to 2017. Through mixed model analysis, we found that the optimum Mean Annual Temperatures were 17°C–20°C for D. frontalis in low-elevation pine-oak forest, while D. mexicanus had two optimal intervals: 11–13°C and 15–18°C. Higher atmospheric Vapor Pressure Deficit (≥ 1.0) was correlated with higher D. frontalis abundances, indicating that warming-amplified drought stress intensifies trees’ vulnerability to beetle attack. As temperatures and drought stress increase further with projected future climatic changes, it is likely that these Dendroctonus species will increase tree damage at higher elevations. Pine forests in Mexico are an important source of livelihood for communities inhabiting those areas, so providing tools to tackle obstacles to forest growth and health posed by changing climate is imperative.


Introduction
Mexico is considered the center of origin and diversification for Pinus tree species [1], and there is evidence of a parallel, co-evolved, regional diversification of the bark beetle genus Dendroctonus (Coleoptera: Curculionidae: Scolytinae) [2,3]. Multiple species of Dendroctonus bark beetles are aggressive tree-killing pests of pine forests in North and Central America [4][5][6]. In Mexico the genus Dendroctonus contains abundant and diverse bark beetle species, with 13 beetle species that colonize 24 of Mexico's 47 pine species [3,7]. An assessment of pine vulnerability to bark beetles in Mexico's forests shows that high pine diversity is a major determining risk factor at the regional level; the most affected regions are the Transverse Volcanic Belt, followed by the Sierra Madre Occidental and Sierra Madre del Sur [5]. At the local level, anthropic and natural disturbance are the key risk factors behind vulnerability [8].
Native bark beetles contribute to the maintenance of ecologically healthy and vigorous forests by removing old, sick, and stressed individuals [9,10]. In this role, bark beetles represent a natural and constant disturbance force in the forest [11,12]. Stand factors such as age structure, species distribution, basal area, site index, etc., influence the risk of large bark beetle populations by acting on host tree resource availability and even on production and mobility of aggregation pheromone plumes that concentrate beetle mass-attack colonization of vulnerable trees [9,11,13]. Successful colonization and reproduction of bark beetles in living trees requires the release of enough aggregation pheromone to ensure that sufficient conspecifics are attracted to the host tree to overwhelm its defenses [14]. However, many other intrinsic biological processes of the beetles are involved. For example, the particular timing of beetle flight and the distance they can cover are related to lipid reserves on emergence and lipid metabolism in flight [15,16].
An increasing number of studies have linked climate change to an increase of bark beetle outbreaks, especially heat and drought, which are the climate factors most clearly causally linked to host tree stress [11,[17][18][19][20][21][22]. Host tree stress makes them more vulnerable to insect pests, particularly bark beetles [13,[23][24][25][26]. There is evidence that climate change is now altering the life cycles and distributions of multiple bark beetle species around the world [27][28][29]. For Mexico studies are scarcer, but bark beetle outbreaks appear to have worsened since 1970, and especially strongly since about 2010 [21,36]. For example, between 2011 and 2013, there were unusually large bark beetle outbreaks of D. mexicanus Hopkins, D. frontalis, and Ips lecontei Swaine (D. Cibrián-Tovar, Universidad Autónoma Chapingo, México, pers. com.). During 2010-2012 the state of Durango lost 25% of its forest production and between 2012 and 2013, 43,859 hectares of forest in Chihuahua were damaged by bark beetles [37]. Similar situations have occurred in Honduras and Guatemala with D. frontalis [38]. After a report of increased bark beetle incidence in Michoacán in west-central Mexico, Rubin-Aguirre et al.
[39] studied Scolytinae communities and found greater beetle prevalence at lower elevations, correlated with higher maximum temperatures.
Despite correlations with temperature, the role of climate in bark beetle outbreaks is not yet fully understood [40]. Turchin et al. [41] analyzed 30 years of D. frontalis population data in East Texas, revealing that these bark beetles respond to a delayed density-dependent process. However, Friedenberg et al. [33], working with models that combine both density-dependent population regulation and climate variables over a large area across 48 years, found that population growth declined with the number of days exceeding 32˚C and that density dependence also decreased with the number of extremely cold (-4.32˚C) days. These results suggest that although increasingly warm temperatures (due to climatic change) could favor larger bark beetle populations, there is an upper limit of favorable temperatures, and importantly, extreme events (e.g., heatwaves) may have a negative effect [42].
Monitoring bark beetle populations while they are flying has been done successfully using traps baited with their own aggregation pheromones [43]. These pheromone traps are routinely used in government programs in Mexico and Central America [44]. Based on that previous experience, in 2014 the Mexican National Forest Commission (CONAFOR, for its acronym in Spanish), along with the Mexican Council of Science and Technology (CONA-CyT), promoted and funded studies on the abundance of the three most important pine barkbeetle species in the country-D. frontalis, D. mexicanus and D. adjunctus Blandford [3,45]using pheromone-baited traps along altitudinal gradients [46]. Initial findings of this project showed that the most critical months for the abundance of flying Dendroctonus spp. are from March to May, corresponding with the warm dry season [47][48][49][50]. Altitudinal patterns of several Dendroctonus species are also known from entomological collections across the country [8,51], where D. frontalis is primarily distributed at lower elevations (below 2,000 m a.s.l.) than D. mexicanus. This finding is consistent with studies using pheromone-baited traps at different elevations in Querétaro [46,52].
Although bark beetles outbreaks are an increasing management problem in Mexican temperate forests, there are no previous nation-wide studies that simultaneously coordinate trapping of flying bark beetle populations and measure climate variables right at the trapping sites. In order to be able to link the abundance of these populations at a given climatic intervals, a large scale analysis of a possible association between the two was conducted for this work. We aimed to contribute to fill the knowledge gaps, especially related to the association between bark beetle abundance and climate, a critical information under the context of climatic change. We explored the limits of the climatic space occupied by each bark beetle species by conducting a trapping along altitudinal gradients. Given the strong association between temperatures and elevation on the complex Mexican mountain ranges (see Fig 2 on [53]), we developed a sampling strategy that provides an array of a diverse occurrence of temperatures and precipitations on each of the studied locations. This sampling approach along elevational gradients, could be applied outside México, like in Guatemala, Belize, Honduras and Nicaragua, that not only share the same bark beetle problems, but also share several keystone pine species with México at low and high elevations, like Pinus oocarpa and P. hartwegii, respectively, as well as Dendroctonus species like Dendroctonus frontalis, D. adjunctus, D. mesoamericanus, D. valens and D. aproximatus, among others [7,54,55].
Thus, we present a nation-wide analysis of the CONAFOR-CONACyT project on the abundances of D. frontalis and D. mexicanus across elevational gradients and the relationship between beetle abundance and local measurements of air temperature and relative humidity and estimates of precipitation and evapotranspiration. The present work aims to determine the climatic envelopes occupied by each of these bark beetle species at the broadest scale across Mexico, rather than building site-based predictive models that would require more detailed parameterization of site-specific conditions (e.g., stand characteristics, insect populations). After determining these national-scale climatic envelopes, we discuss the possible consequences for bark beetle abundance in Mexico under ongoing and projected future climate change conditions.

Study sites
The study was conducted in the temperate pine and pine-oak forests of 11 Mexican States distributed in the northern, central and southern parts of the country (Fig 1). The aim was to sample populations of Dendroctonus frontalis and D. mexicanus species along altitudinal forest More details of transects on S1 to S11 Tables in S1 File. The shapefile was obtained from the as TIFF with a 15m resolution. Source: INEGI, Continuo de Elevaciones Mexicano 3.0 (CEM 3.0), Ver: 202303311214 (https://www.inegi.org.mx/app/geo2/elevacionesmex/). The figure and contour lines were drawn by E.Gómez-Pineda in QGIS 3.10 (https://qgis.org/es/site/). This figure is to be published under the CC-BY 4.0 license. gradients. The range of altitudes covered by the one to three transects per state are summarized in Table 1, and more details on each transect are provided in S1 to S11 Tables in S1 File.
The transects are located in all the four Mexican geographic regions where the occurrence of D. mexicanus was primary studied by Pérez-Miranda et al., [56]: the Sierra Madre Occidental, Sierra Madre Oriental, Transverse Volcanic Belt and Sierra del Sur.
The transects cover the two pine species that conform the extremes of the pine and pine oak forest in México: The transects at Sierra de Arteaga, Coahuila state, includes pine forests dominated by Pinus hartwegii, a high altitude pine that conforms the timberline at México, that reach 4,000 m of altitude and dominate several Natural Protected Areas on the highest mountains of México [57]. On the other elevational extreme, transects at Chiapas state covers a pine forest dominated by P. oocarpa, that neighbor with tropical dry forest, and is also the most important resin producer in México [58]. Pinus leiophylla, a pine widely used also as resin producer [59], is also covered on transects at Durango. Transects at forest dominated by P. durangensis and P. pseudostrobus, at Chihuahua and Michoacán states, respectively, covers two economically very important pine species for commercial timber exploitation: P. durangensis for its large distribution at Sierra Madre Occidental, and P. pseudostrobus for its also wide distribution at Transverse Volcanic Belt at central México and on Sierra Madre del Sur (the last covered by transects at Oaxaca state), a species characterized by a fast growing rate and excellent stem form [60]. Forest dominated by P. greggi, a species widely used in México on degraded pine forest ecological restoration projects [61] is covered by transects on Querétaro state. Forest dominated by P. teocote, that grows frequently on sites relatively marginal in moisture [62], is covered by transects at Nuevo León and Jalisco.
The enlisted pine species covered by the transects, were dominant in the stands where the traps were placed (Table 1). However, there are more pine species present on the sampling sites. It is important to highlight that Mexico is a megadiverse country, center of pine speciation, with about 47 species of the genus Pinus which represent about 42% of the World total number of pine species [1,63,64]. Historically, México has been a country donor of pine species germoplasm as source for pine plantations as exotics, like P. patula in South Africa, to say just one example [65]. Because of that, the Mexican pine forests can be seen as a Worldwide natural reservoir of germplasm for contemporary and potential use in ecological restoration and commercial forest plantations, an asset eventually endangered by the ongoing climatic change [22].

Sampling
The elevation gradient included in the sampling sites covers the main distribution ranges of each bark beetle species: 1500-2000 m a.s.l. for D. frontalis and 2000-2500 m a.s.l. for D. mexicanus [8,51]. A total of 147 sampling sites were established on 24 altitudinal gradients. The lowest and highest elevation at which traps were placed was 1342 and 3611 m a.s.l., respectively. Each altitudinal transect had between 4 and 10 traps, placed sequentially every 100 m of altitudinal change (see S01 to S11 Tables in S1 File for details). The length of the altitudinal transects also aimed to cover the full natural altitudinal distribution of the most economically and/or ecologically important pine species in each studied region. For example, in the case of the State of Michoacán, two transects were placed to cover the natural altitudinal distribution of Pinus pseudostrobus Lind., the most economically and ecologically important pine species of the Mexican Transvolcanic Belt in that State. On each of those two transects, eight trapping sites were placed, covering an altitudinal range of 700 m of elevational difference. The slope of the transects are illustrated for two states with zoom panels in Fig 1. At each of the 147 sampling sites, 12-unit Lindgren funnel traps were placed and baited with a mixture of frontalin, racemic endo-brevicomin and turpentine from Synergy Semiochemicals Corporation, Canada. Release devices for each semiochemical were as follows: frontalin was dispensed from an Eppendorf tube inside of a plastic sachet (release rate 12 mg / day, chemical purity racemic); endo-brevicomin was dispensed from a flexlure (release rate 0.15 mg/ day, chemical purity racemic); and turpentine was dispensed from a single sealed plastic bag (8 × 19 cm) (release rate 3.6 g / day, chemical purity 70:30 mixture of alpha-and betapinene). Semiochemicals were replaced every 45 days in the field. Each Lindgren trap had a collecting cup filled with antifreeze solution (PRESTONA AF EX composed mainly of ethylene glycol) to preserve the beetles [66].
Every two weeks traps were visited, antifreeze was replaced, and trapped insects were collected and stored in 70% ethanol. In Michoacán and Chiapas, only "year 2015" was considered. Referring to the "year" as March through February of the following year corresponds more closely to the actual seasons of ecological activity in Mexican pine and pine-oak temperate forests than the calendar year. The warm dry season occurs from March through May, the rainy season from June to October, and the cold dry season from November to February. Seasonal and annual climate estimates were calculated for these periods. Our operative definition of seasons and start and end of year is in general coincidental with the definition of climatic regions of México by Vidal-Zepeda [67], for the corresponding climatic regions of: North, Northeast, Center, Balsas Depression and Oaxaca Valleys and Southeast. In all those regions, in general, March is a month where temperature increases as a transition from colder winter temperatures to the warmer temperatures of Spring and Summer; also, June to October are months that are included in the rainy season of such regions.
All individuals of D. frontalis and D. mexicanus were identified and counted to estimate the abundance of bark beetle flight activity as a proxy of the abundance of those insects. Taxonomical identification was done following Armendáriz-Toledano et al. [7]. Specimens were collected in accordance with all applicable Mexican regulations, and a sample of the specimens from each locality were deposited at the Colección de Insectos of the Facultad de Biología, Universidad Autónoma de Querétaro. Field experiments were approved by the National Science Research Council (CONACYT) and the Comisión Nacional Forestal (CONAFOR) (Permit number CONACYT-CONAFOR-2014-C01-234547) and by the Dirección General de Vida Silvestre SGPA/DGVS/02594/16. Since we did not worked or interviewed people we did not ask for a informed consent.

In situ climate recording
Air temperature and relative humidity were recorded hourly with data loggers (Hobo Microbaq EL-USB-2) over the same periods as the insect trapping. One data logger was placed adjacent to each trap. Daily average, minimum, and maximum temperature were calculated, which were used to calculate the average, minimum, and maximum monthly values for air temperature and relative humidity, using SAS statistical software [68]. Additional derived variables, such as Mean Annual Temperature (MAT), Mean Temperature of the Coldest Month (MTCM) and of the Warmest Month (MTWM) were estimated in order to represent the climate of each site. We also calculated the mean temperatures over the March-May warm dry season because that was the peak period of bark beetle collection in many traps.

Climate variables from other sources
Available soil moisture is a critical environmental factor influencing tree growth and vigor, and therefore affects resistance to herbivory, especially from bark beetles [13,25,69,70]. Also, relative humidity (atmospheric moisture) is known to affect insect performance by accelerating or delaying metamorphosis and growth rate [71]. Our in situ data collection did not include measurements of precipitation. We therefore complemented our climatic data with additional information from the TerraClimate website (https://doi.org/10.1038/sdata.2017. 191), a globally gridded climate and hydroclimate product that provide (additional than other web sites resources) climatic variables useful to explore the environmental stress induced by hotter droughts, which is a feature characteristic of the ongoing climatic change [22,72]. From that source, we added the following spatially-interpolated climate variables to our analyses of each trap site for each month: vapor pressure deficit (VPD, effectively a measure of atmospheric drought, or the air's demand for moisture); climatic water deficit (CWD, the difference between precipitation and potential evapotranspiration); soil moisture (SOIL_M, mm of soil extractable water); precipitation (PPT); the Palmer Drought Severity Index (PDSI, a long-term index of soil drought); and the Standardized Precipitation-Evapotranspiration Index (SPEI). From the monthly values per trap, we estimated annual means and means for the warm dry season in March-May. The main climate variables used in this study are described in Table 2.
In order to explore if the studied years (2015-2016) were somehow atypical or not in the context of a larger time period, we examined the Mean Annual Temperature from the center of each transect for each year from 1958-2021 with data obtained from TerraClimate website.

Statistical analyses
To visualize how the Dendroctonus species are distributed across the climate space of our study, we constructed histograms of the sum of insects per species per trap, and then summed those values within each one degree Celsius of Mean Annual Temperature.
Next, we fit a response function curve for each species using abundance of the Dendroctonus species per year per trap as the dependent variable, as a function of each climatic variable enlisted on Table 2, either temperatures recorded in situ (obtained by averaging first by day and then by month the raw temperature data of the field dataloggers) or climate variables related to moisture obtained from TerraClimate (as described earlier). We used a quadratic mixed model (Proc Mixed of SAS, 2021), following a similar statistical analysis to the one conducted by Leites et al. [73,74] and Sáenz-Romero et al. [75,76]. The underlying assumption of the quadratic model was that Dendroctonus species would have an optimum climate interval where it would be most abundant, with lower abundances at either colder or warmer temperatures than the optimum [77,78].
The model was: Where: Y ijkl = observation of number of insects trapped of a given species. β 0 . . . β 4 = Regression parameters. C i = i th climate variable. S j = j th State. T(S) k(j) = k th Transect nested in the j th State, and e ijkl = error term. The climate variable (C i ) was a fixed term, while state and transect nested in state were the random terms.
Year of observation (March 2015 to February 2016, March 2016 to February 2017) was considered as a repetition and its effect was also included in the error term.
The best model (i.e., the set of climate variables that best fit the response function curve), was chosen based on a combination of three criteria: (a) the lowest value of the Akaike Information Criterion (AIC, [79]), (b) the climatic quadratic term needed to be negative, and (c) at least one of the climatic terms should be statistically significant (P � 0.05). Table 2. Climatic variables used in this study. These climate variables also were estimated for the dry and warm season of March-May. The first three variables were estimated from data collected with dataloggers in situ, and the others from TerraClimate (https://doi.org/10.1038/sdata.2017.191). After a preliminary analysis, we found that for each species there were many traps with zero individuals or very low abundance (< 50 individuals) of a given Dendroctonus species, likely because the geographic and altitudinal location of such traps were outside the main distribution area of the given species. Thus, to avoid an undesirable bias of fit on response functions due to excessive weight of the zero or very low abundance values, we decided (based on local expertise about each species) to set a minimum abundance of captured insects per year per trap to be included in the mixed-model analyses as follows: 50 individuals for D. mexicanus and 500 for D. frontalis; such criteria was applied when doing the analysis per Dendroctonus species.

Distribution of Dendroctonus spp. abundance along the climatic space
We found the greatest abundance of D. frontalis between 17 and 20˚C MAT, with an optimum at 18˚C. The number of individuals trapped between 18.0 and 19.0˚C in one year reached 146,388, by far the highest abundance found during our study (Fig 2; for a histogram by species see S1a Fig).
We observed a less clear distribution pattern for D. mexicanus in relation to temperature. The highest abundances ranged between 11 to 13˚C of MAT, with an optimum at 12˚C (which mostly corresponds to conifer forests). A second interval of high abundances appeared at warmer sites, between 15˚C and 18˚C of MAT (corresponding to pine-oak forest), with an optimum at 16˚C. Note that the second distribution interval of D. mexicanus widely overlaps with that of D. frontalis (Figs 2 and S1b).

Climatic response functions
When fitting the mixed model of the bark beetle abundance per species captured per year per trap against the climatic variables, we found that the climatic variable that best fit the quadratic model for D. frontalis was Vapor Pressure Deficit (VPD, Table 3, Fig 3). Larger values of VPD were associated with greater abundances of D. frontalis per trap per year, until declining somewhat at the highest VPD values (Fig 3). Since VPD increases nonlinearly with warmer temperatures, these data show how D. frontalis is thus highly temperature dependent.
In contrast, for D. mexicanus, VPD was not the most important climatic variable explaining its abundance. The best mixed model considered Mean Annual Maximum Temperature and the highest peak of D. mexicanus abundance occurred between 20˚C and 26˚C (Fig 4, Table 4). Maximum temperatures either below 20˚C or above 26˚C had lower abundance of D. mexicanus.
State and Transect nested in State as sources of variation and as random effects were not significant (P > 0.05) for both Dendroctonus species (Tables 3 and 4).
We also explored grouping States as Regions (as random effect), where (same grouping as in Table 1

Relationship between climate variables and bark beetle abundance
The abundance of bark beetle flying populations was related to climate variables: Vapor Pressure Deficit and Mean Annual Temperature for D. frontalis, and Mean Maximum Annual Temperature and Mean Annual Temperature for D. mexicanus. The D. frontalis abundance pattern clearly shows that abundance of this beetle is highest at sites with a Mean Annual Temperature (MAT) of 18˚C (Fig 2) and Vapor Pressure Deficit � 1 (Fig 3). Vapor Pressure Deficit can be seen as a proxy indicator of atmospheric drought [72]. Rising VPD has been linked to declining (more negative) water potential inside trees [80], mechanistically linking increasing VPD with physiological drought stress in plants [81]. Those two variables (MAT and VPD) that promote drought are also known to be important for bark beetle outbreaks in xeric habitats [82,83], although such association with specifically VPD have not been proven previously in Mexico. During the warm dry season (March-May) lower altitudes experience higher temperatures in the center of the Transverse Volcanic Belt in Mexico-the same region and elevation (between 600 and 1500 m a.s.l.) where the main distribution of D. frontalis occurs in the country [8]. Our results are similar to findings in the Sierra Gorda mountains of Querétaro state by López-Gómez et al. [84], which found that D. frontalis abundance was positively correlated with temperatures and temperature/precipitation ratio (a climate variable that, like VPD, expresses the balance between temperature and precipitation). For D. mexicanus, the highest relative abundances occur at a lower Mean Annual Temperature than for D. frontalis (Fig 2). This is consistent with the fact that D. mexicanus occurs at higher altitudes, with its main distribution range along the Transverse Volcanic Belt and between 2000 and 2500 m a.s.l. [8]. Our data indicate that a Mean Maximum Annual Temperature between 20 and 26˚C correlates with higher abundances of D. mexicanus, as shown in Fig 4 (notice that D. mexicanus occupy a climatic space colder than D. frontalis, because Mean Maximum Annual Temperatures are not fully comparable with the Mean Annual Temperatures, since the former have overall by far higher values than the later). Our results are also similar to findings by López-Gómez et al., [84], where higher temperatures were associated  Table 2.
https://doi.org/10.1371/journal.pone.0288067.g003 with a higher abundance of flying D. mexicanus in the mountains of Querétaro. VPD was not the most important climatic variable for D. mexicanus abundance, perhaps because this species occurs at higher elevations than D. frontalis, where there is a narrower range of VPD values. Note that another study with a bark beetle species occurring in cool boreal forests found that they are not responsive to climatic moisture deficit [85]; this hypothesis warrants further investigation. Also, it is worth considering that the semiochemicals used in this study were designed  Table 3.
https://doi.org/10.1371/journal.pone.0288067.g004 Table 4. Mixed model analysis for total year abundance of Dendroctonus mexicanus per trap. Akaike Information Criterion (AIC, Akaike, 1973), estimated parameters, contribution to total variance (of random terms) and significance for the selected mixed model are indicated.

PLOS ONE
and tested for monitoring D. frontalis [86]. We included D. mexicanus in our analyses because it consistently appeared in our traps, but we cannot ensure that the semiochemical combination and concentrations were optimal for D. mexicanus. Our conclusions in this species should therefore be considered preliminary, and further detailed investigation is needed on the level of attraction produced in this species by the semiochemicals used. Our findings of a strong relationship between weather and Dendroctonus species distribution are consistent with the results of Moser et al., (2005) from pheromone-baited traps for Dendroctonus species in Arizona (the northern distribution limit for D. mexicanus). They reported that maximum temperatures during bark beetle flight fluctuate from 19˚C in winter to approximately 38˚C in summer. In our study, captures of D. mexicanus individuals dropped at temperatures below 16˚C, and even D. frontalis individuals were less active at cold temperatures [87].
The It is important to stress that our study provide evidence of the climatic intervals on which flying populations of Dendroctonus frontalis and D. mexicanus occur at locations southern that the ones reported for most of the comparable available literature for the genus Dendroctonus, generated mostly at USA and Canada.
The fact that both sources of variation, State or Transect nested in Sate, were not significant, indicates that most of the climatic effects expressed at State and Transect were likely absorbed first by the climatic fixed terms of the model. In other words, State and Transects nested in States should be viewed as sources of variations that contain factors of variation other than climatic effects. The large significance of the climatic fixed terms and the lack of significance of States and Transects as random effects, indicate that the abundance of the two species of Dendroctonus is influenced by some specific intervals of climatic values, independently if those climates occurring at the Northern Mexican states (as Chihuahua, Durango, Coahuila or Nuevo León), or at the Southern states (as Oaxaca or Chiapas).

Possible consequences under future climatic change scenarios
Our results indicate that if temperature continues to increase due to climate change, abundances of D. frontalis will likely continue to increase in numbers at lower elevations (at least until temperatures surpass their optimum) and/or move to higher elevations as their optimal conditions are reached at higher altitudes that were previously too cold for them, as reported for other species under climate change conditions [32,88]. Consequently, these species may pose an increasingly serious threat to the health of Mexican pine forests, especially at their low-elevation xeric limits (generally the extreme lower limit of the natural distribution of each tree species, sensu Mátyás [89]). From our data, we infer that since D. mexicanus has a large area of sympatry with D. frontalis, it would likely also move to higher elevations. Our data do not rule out the possibility that warming temperatures might exceed the optimum values for these insects, and thus their populations might decline beyond that warmer temperature threshold; however, those higher temperatures could also be deleterious to their host tree populations. For example, other studies of insect-plant interactions have documented that temperatures above herbivores' thermal limits can decrease their populations directly, at the same time as extreme temperatures may decrease plant quality for herbivores and therefore decrease herbivore performance indirectly as well [90].
The amount of moisture available for the trees also seems to play a role in the abundance of bark beetles in general. As atmospheric drought stress intensifies with warming temperatures and rising VPD, tree water potentials decline, weakening the trees and inhibiting their defense responses [15,25,81]. This may be linked to an increase in emission of volatile organic compounds, which are known to play important roles in signaling and aggregation of insects [91,92]. Likely the trees are weakened and thus easier attacked by the bark beetles [81].
For both bark beetle species, it is likely that as temperature increases and precipitation decreases due to climatic change [93], they will follow their optimal temperature ranges, shifting their most serious outbreaks to higher-elevation forests than presently observed. Several studies have also found that Scolytinae abundance is related to higher ambient temperatures, which coincides with the lower altitude interval of its distribution [39, 94,95].
As stated by Biederman et al. [28], the drivers of bark beetle population irruptions and crashes are still not fully understood. These beetle population dynamics are apparently driven by complex interactions among tree host abundance and susceptibility, beetle density, weather, symbiotic associations, and natural enemies, with all these interactions affecting the potential of bark beetles to function as major agents of widespread forest disturbance [11]. Thus, further research is needed on these biotic interactions with bark beetle populations, in addition to associations with climatic drivers.
To extrapolate these correlation patterns between bark beetle abundance and temperature, we should also consider how pine hosts may or may not shift their current distributions in future climates [24]. One scenario is that these pine populations stay where they currently are distributed, and therefore host pines would be more attacked in sites where temperature is optimal for bark beetles. If these attacks are sufficiently severe and recurrent, forest populations may tend towards local extinction at those sites. Also, it has been reported that naïve pine populations (not exposed to bark beetles previously) are less able to defend themselves, and therefore may be more susceptible to Dendroctonus attack (reviewed by Raffa et al. [96]). A hypothetical alternative situation is that host pine species start to shift their distributions in response to climate change and bark beetles follow them without causing major distress. However, this seems unlikely because the natural migration speed of forest tree species is much slower than the migration rate projected to be necessary to keep these tree species coupled to the climate conditions to which they are adapted [97][98][99]. A third scenario is that once preferred pines disappear from a particular site, the resident bark beetles will colonize new pine species. This is more complicated and uncertain, since the ability to colonize a specific pine host has been linked with bark beetle fungal symbionts [100,101]. Of course, these hypothetical scenarios need to be taken with caution, since more sophisticated, predictive, site-specific models considering age, species composition and structure of the stands, are needed for a more robust prediction.

Limitations of this research
There is a large body of biological and ecological information on D. frontalis (Coulson and Klepzig [102] and references therein). This includes a particularly wide breadth of knowledge on its pheromonal communication system. For D. mexicanus, on the other hand, information is extremely limited [44]. The pheromone blend used in our work was designed for D. frontali [44,86]. However, it is also reported to attract individuals of D. mexicanus in Mexico and the United States [46-50,52, 87,103]. Therefore, we are confident that the conclusions we have reached about the relative abundance of the trapped D. frontalis numbers in this study are solid, but we are somewhat more cautious in our interpretations concerning D. mexicanus. Although more confident interpretations for trap data on D. mexicanus relative abundance will be obtained someday when an appropriate species-specific pheromone blend becomes available and when the basic ecology of this species is better described, the nationwide study reported here certainly advances our knowledge of D. mexicanus distribution and abundance.
This determination of how climate contributes to the dynamics of these two important Dendroctonus species is a useful advance but represents only an initial understanding of what are likely very complex drivers of both the distribution and abundance of Dendroctonus populations as climate change unfolds. Further developments will require more sophisticated models that consider specific characteristics of current forest stands (e.g., pine species present, basal area, and recent history of growth trends and beetle outbreaks) and account for potential shifts in host tree distributions [104][105][106][107]. However, such finely parameterized models may only prove effective at local or regional scales, which highlights the importance of less complex, broad-scale empirical models like the one used in this study.
Monitoring bark beetle populations using pheromone-baited traps (as per Hayes et al. [108]) proved to be an effective approach to collect beetle population data for this study, as this approach provided reliable information on variation in species presence and abundance across elevational gradients using trap transects that sampled flying beetles. However, Hayes et al. [108] also indicate an inconsistent relationship between flying beetle population data and concurrent host tree mortality. In other words, sampling flying beetles at the same time as tree mortality may underestimate the impacts of beetles on trees because processes like tree mortality often lag significantly behind beetle flights and infestations. Thus, long-term monitoring like the successful D. frontalis monitoring system in the Southern United States may be required [43]. In our sites, the correlation between bark beetle abundance and tree infestation was not recorded; these relationships need to be addressed in future studies.

Future research needs
From a different perspective, bark beetles can be considered an ecosystem component that promotes the renewal of aging stands, but this role needs to be better understood in the context of climatic change. The number and intensity of outbreaks are undoubtedly increasing due to debilitated trees because hotter droughts are occurring due to climate change [21,22]. Unfortunately, at the same time, in the case of the main mountain systems in México, recruitment via natural regeneration has been more difficult because the March-May warm dry season is becoming ever warmer and drier due to climate change [109]. Such a combination of processes poses an important challenge for both commercial forest management and for the conservation of Natural Protected Areas to achieve stand renewal after bark beetle outbreaks.
There is an important body of data accumulating in recent years, that indicates an association of increasing annual temperatures and hotter drought periods with a larger number of bark beetles outbreaks that need a sanitary intervention by logging. In our view, this is clear evidence of the association between the ongoing climatic change and the increase of bark beetle populations [21]. However, a nation-wide monitoring effort is needed to trap bark beetle on sites with disparate climatic characteristics, as we have presented here, but for longer periods of time. That in order to have a stronger support of the association between increasing temperatures and drought periods with the occurrence of more abundant bark beetle flying populations.

Conclusions
1. Bark beetles appear to have optimal temperature conditions: between 17˚C and 20˚C of Mean Annual Temperature for Dendroctonus frontalis and two optimal intervals for D. mexicanus: between 11˚C and 13˚C and between 15˚C and 18˚C.
2. Drought conditions seem to favor increases in bark beetle populations, as suggested by the association between higher D. frontalis abundance and higher Vapor Pressure Deficit values.
3. Because environmental conditions are shifting due to human-caused climate change, and historical climate for a given site is progressively occurring at higher altitudes in mountain regions, it might be possible that bark beetles will likely follow their optimal climate conditions by shifting to higher altitudes.
4. Given the projected increase of temperatures and decrease of precipitation for Mexican continental regions due to climatic change, and considering the association between high values of vapor pressure deficit and values of mean maximum annual temperatures with the abundance of D. frontalis and D. mexicanus, respectively, it seems likely that bark beetle populations will increase along with the increase of temperatures and decrease of precipitation, which may lead to worsening infestations and affect pine fitness. (TIF) S1 File. S1 to S11 Tables (one Table per