Highlights:
Climate, fire, tree density, and topography influence Dendroctonus spp. infestations.
Current and future (2015-2039) spatial distribution of bark beetle infestation risk was modeled.
High risk of D. mexicanus infestation increased from 3.9 % (current scenario) to 5.0 % (future scenario).
High risk of D. frontalis infestation decreased from 10.8 % (current scenario) to 9.6 % (future scenario).
The northeast of Michoacán has the highest risk of bark beetle infestation.
Introduction
Conifer bark beetle insects generate changes in forest ecosystems when they reach epidemic levels leading to massive tree death (Krist, Sapio, & Tracz, 2007; Morris et al., 2018). Under such conditions, insects affect wood supply and ecosystem services for human well-being (Dhar, Parrott, & Heckbert, 2018).
In Western United States and Canada, infestations by Dendroctonus ponderosae Hopkins were recorded for more than a decade on about 28 million hectares (Dhar et al., 2018; United States Department of Agriculture [USDA] Forest Service, 2018). In Mexico, conifer bark beetle insects (Dendroctonus spp. Erichson, Coleoptera: Curculionidae: Scolytinae) are considered the most damaging forest pest because they cause tree death (Soto-Correa, Avilés-Carrillo, Giron-Gutiérrez, & Cambrón-Sandoval, 2019). Nevertheless, these insects are considered biotic regulators of forests, as they play important roles in forest stand regeneration (Dhar et al., 2018; Raffa, Grégoire, & Lindgren, 2015) and in determining its structure and composition (Raffa et al., 2008).
Dendroctonus spp. populations are favored by low precipitation (Chapman, Veblen, & Schoennagel, 2012) and high temperatures (Raffa et al., 2015; Soto-Correa et al., 2019), as well as by occurrence of forest fires (Billings et al., 2004; Fonseca, Llanderal, Cibrián, Equihua, & de los Santos, 2009). On the other hand, large, homogeneous and mature stands are also more susceptible to insect attack (Raffa et al., 2015).
Seidl et al. (2017) believe that, in the future, forests will be at greater risk and frequency of disturbance due to climate change, with the possibility of increases in bark beetle insect populations beyond their historical range of spatial variability (Bentz et al., 2019). Therefore, early detection of forest areas susceptible to bark beetle insect attack would benefit mitigation of effects of these potential infestations (Bone, Wulder, White, Robertson, & Nelson, 2013).
Recently, temperate forests in Michoacán have been attacked by Dendroctonus mexicanus Hopkins and Dendroctonus frontalis Zimmermann bark beetle insects (Rubin-Aguirre et al., 2015; Salinas-Moreno et al., 2004). About 58 191 m3 of wood were affected in 2017, which represented 14.78 % of the total pine wood production in the state.
The objective of this research was to model the potential spatial distribution under current climatic conditions and an assumed scenario (2015-2039) of climate change and thus generate a tool that can be used to determine probability of epidemic outbreaks of D. mexicanus and D. frontalis in temperate forests of Michoacán.
Materials and Methods
Study area
The study area comprises temperate coniferous forests of Michoacán, Mexico, covering approximately 1 124 973 ha (Comisión Forestal del Estado de Michoacán [COFOM], 2015). These forests are part of the physiographic provinces Trans-Mexican Volcanic Belt (TVB) and Sierra Madre del Sur (SMS). The region has temperate climates, with mean annual temperatures ranging from 8 to 30 °C, and annual precipitation from 400 to 2 000 mm (Instituto Nacional de Estadística y Geografía [INEGI], 2018). The altitude varies from 104 to 3 641 m (INEGI, 2019).
Selection and analysis of variables associated with bark beetle outbreaks
Within the framework of the analytic hierarchy process (AHP; Saaty, 1980), criteria (variables) were organized, analyzed and combined to estimate the risk of bark beetle outbreaks. The AHP comprises three phases (Malczewski, 1999): 1) construction of a hierarchy problem based on relevant criteria and sub-criteria, 2) analyzing criteria and sub-criteria (variables) by experts in the field, and 3) estimation of a consistency value associated with analyzing criteria and sub-criteria to ensure that weights assigned are not the product of chance, but of a rational and correct analysis.
Phase 1
Relevant criteria and sub-criteria for the outbreak of D. frontalis and D. mexicanus were identified and defined based on a literature review and consultation with seven Mexican forest pest experts (Malczewski, 1999). The literature review provided information on geographic distribution (Salinas-Moreno et al., 2010), physiographic characteristics (Armendáriz-Toledano, Zuñiga, García-Roman, Valerio-Mendoza, & García-Navarrete, 2018; Salinas-Moreno et al., 2004), climate conditions (Morales-Rangel, Cambrón-Sandoval, Soto-Correa, Jones, & Obregón-Zuñiga, 2018; Raffa et al., 2015) and tree-size measurements of the sites where bark beetle outbreaks occur (Krist et al., 2007; Negrón, Anhold, & Munson, 2001; Raffa et al., 2015). The criteria identified as relevant were climate, occurrence of wildfire, stand tree density (basal area), and topography. The sub-criteria total annual precipitation, mean annual temperature, exposure, altitude and slope were identified. The hierarchy problem (Figure 1) was constructed from these criteria and sub-criteria (Malczewski, 1999).
Phase 2
A total of seven surveys were completed by experts with knowledge on forest pests to provide, by means of paired comparison matrices, each criterion and sub-criterion identified as relevant. Respondents were instructed to use the 1 to 9 rating scale proposed by Saaty (1980), which specifies that by assigning a higher number to a factor/criterion, it represents more importance than another for the purpose being evaluated. The matrices were solved in the WEIGHT module of TerrSet® version 18.31 (Eastman, 2016) to calculate the importance value (weighting) of each criterion and sub-criterion.
Phase 3
The consistency ratio (CR) defined by Saaty (1980) was calculated to evaluate the weighting confidence. A CR < 0.10 indicates that consistency is acceptable. If CR > 0.10, inconsistency in the comparisons is assumed and should be corrected or discarded (Malczewski, 1999; Saaty, 1980). CR value for each survey was calculated in the WEIGHT module of TerrSet® version 18.31 (Eastman, 2016). Of the seven surveys conducted, only three that were consistent were considered and, from these, the arithmetic mean was calculated as the weighting value for each of the criteria and sub-criteria. These values were incorporated as multiplicative vectors in the multi-criteria evaluation process to estimate risk of outbreak.
Spatial databases
Previously generated maps of total annual precipitation and mean annual temperature (~1 km2 spatial resolution) were used for the current scenario (Cuervo-Robayo et al., 2014) and a future scenario (Fernández, Zavala, Romero, Conde, & Trejo, 2015). The maps for the current scenario were constructed with average monthly climate data for the period 1910-2009. Those corresponding to the future scenario, called representative concentration trajectory (RCP 4.5) with near horizon 2015-2039, resulted from modeling a low greenhouse gas concentration trajectory (580-650 ppm CO2eq) that could be reached by the year 2100 and in which temperature changes of approximately 2.3 to 2.6 °C may occur. The RCP 4.5 scenario is one of the four greenhouse gas concentration trajectories adopted by the Intergovernmental Panel on Climate Change (IPCC, 2014) in climate change scenarios and which has more coincidence with cumulative emissions in the short term (Hausfather & Peters, 2020). On the other hand, geo-referenced data of forest fires recorded by the National Forestry Commission and the Comisión Nacional Forestal and the Sistema Nacional de Información y Gestión Forestal (CONAFOR & SNIGF, 2020) were obtained for the period 2010 to 2019.
Tree density criterion was represented by a basal area map (m2∙ha-1), based on data from the Inventario Nacional Forestal y de Suelos (INFyS) (CONAFOR & SNIGF, 2020b). The existing basal area at INFyS sites was determined with ratio-mean estimators (Šmelko & Merganič, 2008); subsequently, basal area was mapped over the entire area comprising the state's coniferous forests using the weighted inverse distance interpolation method (Chirinos & Mallqui, 2016). Finally, slope, aspect and altitude maps were generated using the digital elevation model with 30 m spatial resolution (INEGI, 2019). Maps were constructed in ArcGis® 10.5, homogenizing pixel size to 30 m, and cartographic projection to geographic coordinate system and the North American Datum of 1927 (NAD 27).
Standardization of sub-criteria
Criteria and sub-criteria identified as relevant for bark beetle infestation are recorded in different units of measurement, so it was necessary to standardize them to a common unit, i.e., to values on a scale from 0 to 1 (Gómez & Barredo, 2005). For this purpose, fuzzy membership functions were used (Malczewski, 1999). First, inflection points of the functions were defined (Table 1), which were used as a basis for defining the degree of membership (suitability) of each pixel that makes up the study area. The degree of membership is higher when the standardized value is close to one and lower when it is close to zero. Membership functions can have varied shapes, the most common are linear, sigmoidal and “J” shaped with increasing, decreasing or symmetric direction (Eastman, 2016). The “fire” criterion was standardized using binary values (0 or 1), because information regarding the degree of intensity was not available. Areas with and without fire presence were classified with value equal to one and zero, respectively. Standardized values of the sub-criteria were generated in the Fuzzy module of TerrSet® version 18.31 (Eastman, 2016).
Risk estimation for current infestation and under a climate change scenario
The risk of outbreak for both species of bark beetles was estimated considering four criteria and five sub-criteria (Figure 1). After sub-criteria were standardized using membership functions (Table 1), they were weighted by multiplying by the weight assigned by the experts (Figure 2). Subsequently, the weighted sub-criteria maps were summed to give rise to criterion maps. Using a similar procedure, the weighted criterion maps were generated and multiplied by the weight assigned by the experts. Finally, criterion maps were summed to create the map that estimates the risk of bark beetle outbreak at the pixel level (Figure 2). For easier interpretation, the calculated values were grouped into four risk classes: low (0-0.2), medium (0.2-0.4), high (0.4-0.6), very high (>0.6).
The procedure to generate the future risk map was similar to that explained in the previous paragraph to estimate the risk of infestation by Dendroctonus under current conditions; the only difference was to replace the climate criterion maps (sub-criteria total annual precipitation and mean annual temperature) with maps that record the projection of these sub-criteria considering the climate change scenario RCP 4.5. This shows, in general, increase in temperature and decrease in precipitation (Fernandez et al., 2015). The rest of the criteria remained constant.
The risk map generated for each bark beetle species was compared with the geographic coordinates recorded for bark beetle outbreaks in the period 2010 to 2019.
Comparison of current and future spatial distribution of infestation risk
Infestation risk maps for the current scenario and the future scenario (climate change) were compared using the CrossTab module of TerrSet® version 18.31 (Eastman, 2016). This module generated a cross-classification matrix that identified changes in estimated infestation risk values for the two modeled scenarios: gains, losses, permanence, and total change in hectares, for each of the infestation risk categories defined for the two Dendroctonus species (Martínez, Martín, & Román, 2000). In addition, CrossTab provided the Kappa index, a quantitative measure of agreement between the compared maps (Eastman, 2016). This index takes values from -1 to 1; values close to 1 indicate strong agreement between maps, on the contrary, values close to -1 indicate low agreement (Martínez et al., 2000).
Results and Discussion
Criteria and sub-criteria identified for bark beetle infestation
Four relevant criteria were identified that influence infestation of D. mexicanus y D. frontalis. The order of importance is as follows: climate, fire, tree density and topography (Table 2). At the sub-criteria level, total annual precipitation, mean annual temperature, aspect, altitude and slope emerged as relevant.
Factor | Membership function | Inflection points used on x-axis |
---|---|---|
Total annual precipitation (mm) | Sigmoid symmetric | a: 500, b: 550, c: 700: d: 1 200 D. mexicanus |
a: 700, b: 750, c: 900: d: 1 100 D. frontalis | ||
Mean annual temperature (°C) | Sigmoid symmetric | a: 5, b: 18, c: 22: d: 26 D. mexicanus |
a: 12, b: 20, c: 26: d: 30 D. frontalis | ||
Basal area (m2∙ha-1) | Sigmoid symmetric | a: 12, b: 15, c:20, d: 28 D. mexicanus |
a: 12, b: 15, c: 20, d: 28 D. frontalis | ||
Altitude (m) | Sigmoid symmetric | a: 900, b: 2 300, c: 2 500, d: 3 300 D. mexicanus |
a: 800, b: 1 600, c: 2 000, d: 3 000 D. frontalis | ||
Slope (°) | Linear symmetric | a: 0, b: 0.5, c: 2, d: 4 D. mexicanus |
a: 0, b: 0.5, c: 2, d: 3 D. frontalis | ||
Aspect (°) | Linear symmetric | a: 90, b: 135, c: 225, d: 270 D. mexicanus |
a: 90, b: 135, c: 225, d: 270 D. frontalis |
Primary literature on bark beetles supports our finding on weighting climate as the main criterion affecting the outbreaks of these insects. According to Chapman et al. (2012) and Raffa et al. (2015), increased temperature and decreased precipitation favor the growth of insect populations. On the other hand, Raffa et al. (2008, 2015) mention that prolonged droughts lead to tree stress, causing them to weaken and become less resistant to bark beetle attack. Morales-Rangel et al. (2018) found that a mean annual temperature between 16 and 18 °C favors growth of D. frontalis and D. mexicanus, while Moser, Fitzgibbon, and Klepzig (2005) report that mean temperatures below 16 °C decreased the presence of both insects. González-Hernández, Morales-Villafaña, Romero-Sánchez, Islas-Trejo, and Pérez-Miranda (2018) indicate mean annual temperatures from 12 to 18 °C and precipitation in the range of 600 to 1 200 mm increase D. mexicanus growth.
Referring to fire analysis, according to Fonseca et al. (2009), trees subjected to fire action release volatile compounds that serve as attractants for some bark beetle insects. On the other hand, Billings et al. (2004) reported that fire weakens trees, causing them to produce less resin and have less defense against Dendroctonus attack.
Regarding the importance of a forest's structure, Raffa et al. (2008, 2015) note that tree density influences bark beetle outbreaks. Morris et al. (2018) discuss that stands with high tree density are more susceptible to infestation due to competition for resources. This agrees with that reported by Negrón et al. (2001), who found more attacks by Dendroctonus pseudotsugae Hopkins in stands of Pseudotsugae menziesii (Mirb.) Franco with basal area of more than 26.4 m2∙ha-1 and with greater intensity in those exceeding 39 m2∙ha-1.
Regarding topography, according to Krist et al. (2007), this variable plays an important role in tree mortality risk. According to Bennie, Huntley, Wiltshire, Hill, and Baxter (2008), slope and aspect influence solar radiation received by trees, affecting evapotranspiration and available soil moisture, mainly those on slopes with southern aspect. As a consequence, these trees suffer greater stress and are more exposed to bark beetle outbreak. Cuéllar-Rodríguez, Equihua-Martínez, Villa-Castillo, Estrada-Venegas, and Romero-Nápoles (2013) reported infestations of D. mexicanus in Pinus cembroides Zucc. forests established at altitudes ranging from 1 820 to 2 700 m and different aspects, including the southwest. Salinas-Moreno et al. (2004) mention a preferred altitudinal range between 2 100 and 2 500 m for this insect. Armendáriz-Toledano et al. (2018) and Salinas-Moreno et al. (2004) reported higher abundance of D. frontalis in altitudinal interval between 1 500 and 2 000 m.
Risk of infestation by Dendroctonus mexicanus and D. frontalis
Two risk maps were generated by insect species, one for the current scenario (1910-2009) and the other for the future scenario (2015-2039) (Figures 3 and 4). The maximum estimated infestation risk value for D. mexicanus was 0.78 for the current scenario and 0.83 for the future scenario. Maximum values of 0.84 and 0.85 were estimated for D. frontalis for the current and future scenarios, respectively. For both species, the areas of highest risk are located along the TVB, which coincides with that reported by Salinas-Moreno et al. (2010) and Armendáriz-Toledano et al. (2018), who indicated that the most favorable climatic conditions for bark beetles under analysis are located above the TVB. Mapped infestation risk is also consistent with reports of outbreak recorded in the period 2010 to 2019 by CONAFOR (Figures 3 and 4); 2 255 records were reported for D. mexicanus and 99 for D. frontalis.
According to Armendáriz-Toledano et al. (2018), in the SMS there are also ideal climatic conditions for both species of bark beetle to growth; however, the results show that the risk of an outbreak in this area is lower than that estimated for the TVB. In another study, Salinas-Moreno, Ager, Vargas, Hayes, and Zúñiga (2010) reported that the most common pine species in the TVB and SMS are very similar, so it may be the case that D. mexicanus and D. frontalis are found in the same area, because they have similar host preferences and can attack and coexist in the same tree (Armendáriz-Toledano et al., 2018; Moser et al., 2005). Both D. mexicanus and D. frontalis are aggressive species that attack various pine species in the northern hemisphere depending on latitude; although in appearance they are generalist species, the main hosts of D. frontalis in Mexico are Pinus oocarpa Schiede ex Schltdl., P. pringlei Shaw, P. teocote Schiede ex Schltdl., P. leiophylla Schiede ex Schltdl. et Chamisso, while for D. mexicanus they are Pinus leiophylla, P. teocote, P. pseudostrobus Lindl., P. devoniana Lindl., P. montezumae Lamb. and P. oocarpa (Salinas-Moreno et al., 2010).
Comparison of current and future infestation risk
Table 3 shows the cross-classification matrix for D. mexicanus; the analysis indicated that 767 013 ha (82.7 %) of the area classified for the period 1910-2009 remained unchanged in the future scenario. The rest of the area (160 194 ha; 17.3 %) underwent a change in the risk category initially assigned (period 1910-2009). The area classified as high risk of infestation increased by 10 135 ha (36 326 to 46 461) and the area classified as very high risk will remain constant in practical terms, as it increased by 508 ha (2 511 to 3 019). In contrast, the area classified in the low-risk category increased from 620 176 ha (66.9 %) in the 1910-2009 scenario to 671 204 ha (72.4 %) for the future scenario (Figure 5), meaning again of 51 028 ha, which results in a decrease in the medium risk category (loss of 61 672 ha). In general, and contrary to common perception, these results indicate that, under the climate change scenario used, the risk of D. mexicanus infestation will be lower in the future; however, gains should not go unnoticed, because, in absolute terms, the area classified as high and very high risk in the future scenario will increase by 10 136 and 508 ha, respectively.
Current period | Future risk (RCP 4.5 scenario) | Current total area (ha) | Losses (ha) | |||
---|---|---|---|---|---|---|
Low | Moderate | High | Very high | |||
Low | 582 992 | 36 726 | 458 | 0 | 620 176 | 37 184 |
Moderate | 88 195 | 158 569 | 21 369 | 61 | 268 195 | 109 625 |
High | 17 | 11 228 | 23 787 | 1 293 | 36 326 | 12 538 |
Very high | 0 | 0 | 846 | 1 665 | 2 511 | 846 |
Total area (ha) RCP 4.5 scenario | 671 204 | 206 523 | 46 461 | 3 019 | 927 207 | 160 194 |
Gains (ha) | 88 212 | 47 953 | 22 674 | 1 354 | 160 194 |
Based on Table 4, the comparison of current and future risk for D. frontalis indicated that 80.54 % of the area (746 784 ha) remained unchanged and that the area classified with the infestation risks of greatest interest to forest managers (high and very high risk) will be similar in the future scenario (Figure 6). In fact, the area with high-risk class decreased by 11 015 ha (100 353 to 89 338) and the very high-risk class decreased by 82 ha (3 901 to 3 819).
Current period | Future risk (RCP 4.5 scenario) | Total current area (ha) | Losses (ha) | |||
---|---|---|---|---|---|---|
Low | Moderate | High | Very high | |||
Low | 641 133 | 33 082 | 13 688 | 0 | 687 902 | 46 770 |
Moderate | 60 607 | 51 442 | 22 606 | 395 | 135 051 | 83 608 |
High | 13 751 | 33 925 | 51 731 | 946 | 100 353 | 48 622 |
Very high | 0 | 110 | 1 313 | 2 478 | 3 901 | 1 423 |
Total area RCP 4.5 scenario | 715 491 | 118 559 | 89 338 | 3 819 | 927 207 | 180 424 |
Gains (ha) | 74 358 | 67 117 | 37 607 | 1 341 | 180 424 |
The Kappa index was estimated at 0.91 for D. mexicanus and 0.89 for D. frontalis. These values suggest a very high concordance (0.81-1.00) (Landis & Koch, 1977) in the spatial distribution of the risk classes mapped in Figures 3 and 4. Kappa values indicate that a significant fraction (the largest part) of the infestation risk classes remain unchanged between the scenarios being compared, which is observed in the graphs in Figures 5 and 6.
Although in Mexico there is an “early warning and risk assessment system for forest pests” (SATERF), managed by the Comisión Nacional Forestal and the Sistema Integral de Vigilancia y Control Fitosanitario Forestal (CONAFOR & SIVICOFF, 2020), the methodological proposal and analysis presented here are new and could contribute to the development of this type of tools. Unlike SATERF, the methodology proposed here defines a hierarchy analyzed of relevant criteria for epidemic outbreaks. Additionally, it incorporates forest density (basal area) as an important factor for predicting the risk of bark beetle outbreaks (Morris et al., 2018; Negrón et al., 2001; Raffa et al., 2015). On the other hand, there is also an improvement in the scale on which the results are presented. SATERF publishes national and state maps per month; at the national level it issues statements which determines high risk in forest areas of Michoacán and by entities it reports high risk in most of the SMS and the southwest of the TVB (CONAFOR & SIVICOFF, 2020b), which makes it difficult to find references for the confirmation or contrast of the results.
The future climate scenario (RCP 4.5) predicts an increase in temperature and a decrease in precipitation. According to Sáenz-Romero, Rehfeldt, Crookston, Duval, and Beaulieu (2012), temperature will increase by 1.4 °C and precipitation will decrease by about 5.6 % in Michoacán by 2030. Considering that the most important criterion for the experts was climate (0.45), an increase in the level of risk was expected in the future scenario; however, precipitation and temperature maps analyzed in data processing indicated that some areas that make up the forests of Michoacán will have higher precipitation and lower temperatures in the future scenario. These forecasts are reflected in the spatial estimation of infestation risk, so that medium risk areas for the current scenario have a larger area than those of the future scenario.
Morales-Rangel et al. (2018) propose an approximate 3.0 °C increase in mean annual temperature by 2030, which would favor the incidence of D. mexicanus and D. frontalis outbreaks in places previously reported to be at low risk. The results of this research give reason to conclude that such an assertion is debatable, because it omits the future occurrence of temperature and precipitation conditions less favorable for the development of bark beetles, which will tend to balance the expected changes for infestation risk in forests of Michoacán.
Conclusions
Maps generated indicate that forests of Michoacán with higher risk of Dendroctonus mexicanus and D. frontalis infestation, for the future scenario (2015-2039), are those located in the Transversal Volcanic Belt; the highest risk of infestation by the two species is predicted in the northeastern part of the state. Continuous monitoring and silvicultural practices that favor adequate tree density are recommended activities to maintain the forest in a healthy condition and prevent outbreaks of pine bark beetles. It is important that density does not exceed the capacity of the site index. The methodology used represents a support tool for spatial analysis of the risk of bark beetle outbreaks and for decision making regarding monitoring, prevention and application of sanitary measures for the control of D. mexicanus and D. frontalis in Michoacán, but applicable to other areas if the required data is available.