Patterns and Drivers of Pine Processionary Moth Defoliation in Mediterranean Mountain Forests

Insect defoliations are one of the threats of forests health, and climate warming may enhance the occurrence of more severe outbreaks. Thus, understanding the spatial and temporal patterns of insect defoliation may provide valuable information to infer their responses to climate change. To contribute to this question, we studied a 42-year long record of defoliation by pine processionary moth (Thaumetopoea pityocampa, hereafter PPM), a major defoliating insect of Mediterranean conifers. We tested the hypothesis that climate warming is increasing PPM defoliation, which occurs in winter, and enhancing its upward shift. We analyzed data corresponding to PPM defoliation affecting 92 stands of four pine species distributed across an altitudinal (1,030–1,695 m) gradient (going downwards: Pinus sylvestris, Pinus nigra, Pinus pinaster, and Pinus halepensis) located near Mora de Rubielos (Teruel, Spain). We found a higher prevalence of PPM defoliation in P. nigra stands followed by P. sylvestris. PPM defoliation showed several peaks (1979–1983, 1991–1992, and 1995–1996) and it was positively related to winter minimum temperatures, with temperatures below −12°C reducing PPM defoliation. We found higher defoliation in stands dominated by P. nigra, located at intermediate altitude and with low tree density. We did not detect any upward shift of PPM defoliations. PPM defoliation dynamics depend on forest structure and on forest composition, but it might also be impacted by the occurrence of lethal temperatures below a threshold. These patterns should be considered when forecasting the range expansion of forest insects in response to climate warming.

Introduction

Even though a recent review of historical reports suggests that anthropogenic pressure threatens the diversity of entomofauna worldwide (Sánchez-Bayo and Wyckhuys, 2019), climate warming may increase the activity of tree-feeding insects, thus triggering more severe and frequent defoliation outbreaks (Battisti et al., 2005; La Porta et al., 2008; Sturrock et al., 2011; Pautasso et al., 2015). However, insects often exhibit cyclic outbreaks and cause periodic defoliation events as is the case of forest moths (Myers, 1993). While some hypotheses suggest that periodicity in moth defoliation is mainly driven by intrinsic factors such as changes in fecundity rates (Myers and Cory, 2013), others point to density-dependent trophic drivers including food quality (Hódar et al., 2004) or natural enemies as predators, pathogens, parasites and parasitoids (Myers, 1993; Liebhold et al., 2000, 2013; Bjørnstad et al., 2010; Klemola et al., 2010). Along these lines, weather oscillations are also considered as external modulators of cyclic outbreaks and defoliation incidence (Hódar et al., 2012; Allstadt et al., 2015). However, defoliation outbreaks also present spatial variability due to the occurrence of different tree species with diverse susceptibility against defoliators (Castagneyrol et al., 2014; Régolini et al., 2014), and the variation in regional climate or local environmental conditions (Li et al., 2015; Toïgo et al., 2017). These sources of variability lead to asynchronous defoliations across regions. Thus, analyzing how climate conditions influence long-term records of defoliation by insects across altitudinal or latitudinal gradients is crucial to understand the potential range expansion of forests months (Myers and Cory, 2013). This is the case of the pine processionary moth (Thaumetopoea pityocampa, hereafter abbreviated as PPM) a major Mediterranean defoliator moth whose expansion to higher latitudes and altitudes has been reported in some regions (see Battisti et al., 2015).

The PPM is one of the major forest pests across the Circum-Mediterranean region where it causes recurrent defoliation and forest damage affecting several pine and cedar species (Palacio et al., 2012; Jacquet et al., 2013; Roques et al., 2015). The PPM is regarded as an insect very sensitive to climate warming since the main feeding activity of larvae occurs in winter (Roques et al., 2015). This explains why links between PPM expansion to higher altitudes or latitudes and warmer winters have been explored by several researchers (Hódar et al., 2003; Hódar and Zamora, 2004; Battisti et al., 2005; Robinet et al., 2007, 2014), albeit human-mediated transport has also been involved in the PPM expansion (Robinet et al., 2012). It might be noted that anthropogenic activities, particularly agricultural practices, have been found to have more important effects than climate change on the decline of entomofauna worldwide (Sánchez-Bayo and Wyckhuys, 2019).

The PPM population dynamics follow a positive gradation phase until maximum PPM density is observed and defoliation peaks (Démolin, 1969). Then, a negative gradation phase occurs which is characterized by a decline in the PPM population due to reduced food availability until a latent phase is reached (Battisti et al., 2015). PPM larvae build silky nests in the crowns of trees which allow caterpillars being active in winter, resisting low air temperatures and feeding on needles (see the PPM cycle in Figure S1). Warmer winter conditions favor larvae performance and increase their survival (Démolin, 1969; Hoch et al., 2009; Battisti et al., 2015). In addition, larvae feeding activity increases with rising winter temperatures so warmer winters are expected to promote PPM expansion and accelerate population growth (Buffo et al., 2007; Robinet et al., 2007). However, very low winter temperatures lethal to PPM caterpillars (lower than −15°C; Démolin, 1969; Buffo et al., 2007) may be recorded in mountain and inland Mediterranean regions (Battisti et al., 2005; Robinet et al., 2014). Several studies suggest that night temperatures permissive for feeding activity (above 0°C) and day temperatures allowing for food digestion (9°C) are enough to account for the northward expansion (Battisti et al., 2005, 2015; Robinet et al., 2007, 2014; Hoch et al., 2009). However, the potential negative influence of the occurrence of lethal temperatures on PPM activity has been less documented (e.g., Buffo et al., 2007).

It has been shown that PPM defoliation outbreaks display contrasting temporal patterns in different regions (Hódar et al., 2012; Li et al., 2015). Hódar et al. (2012) reported oscillatory patterns in the percentage of pine stands with strong PPM defoliations across southern Spain and found a strong relationship between PPM outbreaks and winter climate conditions. Li et al. (2015) suggested that PPM dynamics across the French territory depended on several factors such as low host food quality following defoliation causing a strong negative density feedback. PPM outbreaks are also dependent on the host tree species, with some tree species being more susceptible to PPM defoliation than others. It has been shown that Pinus nigra and Pinus sylvestris are more susceptible of being defoliated than other species such as Pinus pinaster (Hódar et al., 2012; Sangüesa-Barreda et al., 2014). Recently, Castagneyrol et al. (2014) found that tree species mixing influence the degree of PPM defoliation severity and the occurrence of outbreaks. Therefore, interactions between spatial patterns and temporal trends of PPM defoliation must be investigated.

Since management actions to detect and evaluate PPM outbreaks are widely extended in several Circum-Mediterranean countries (e.g., Li et al., 2015), we capitalize on a regional assessment of PPM defoliation carried out in Mora de Rubielos (Teruel, eastern Spain) for a 42-year long period (1970–2011) and across a 665-m altitudinal gradient (Montoya, 1970). Given that stand defoliation, the percentage of attacked trees, and the density of nests observed in the canopy are more noticeable at the end of the winter, these are the variables often recorded in national and regional PPM monitoring programs (Cayuela et al., 2011). Some authors have cautioned against using counts of nests or measures of defoliation which may be difficult to estimate in tall trees, dense stands or at low PPM population densities during the early outbreak stages (Battisti et al., 2015), but these are the most available surrogates of PPM activity (Jactel et al., 2006). Moreover, defoliation incidence has been found to be a valid proxy to account for PPM population dynamics (Cayuela et al., 2014).

Here we: (i) analyze the spatiotemporal patterns of PPM defoliation, and (ii) evaluate whether PPM defoliation varies as a function of tree species, stand density, altitude and winter air temperatures. We hypothesize that warmer winter temperatures positively influence PPM dynamics by enhancing feeding activity, and, consequently, we expect an upward shift of PPM outbreaks because of climate warming unless lethal temperatures are reached and kill the populations.

Materials and Methods

Study Site and Variable Survey

The study region is located near Mora de Rubielos, Teruel, eastern Spain (40° 22′ 25′′ N; 0° 37′ 56′′ W; Figure 1) and it encompasses a 128-km² large area. Initially, the region was selected for studying PPM dynamics because: (i) it covers a wide altitudinal gradient (1030–1695 m a.s.l.); and (ii) this allows the presence of several pine species showing different susceptibility to PPM defoliation. Specifically, four pine species occur in different altitudinal belts along the studied gradient (from low to high elevations): Pinus halepensis Mill., Pinus pinaster Ait., Pinus nigra Arn. subpsp. salzmannii (Dunal) Franco (P. nigra hereafter) and Pinus sylvestris L. (see also Camarero et al., 2015). It should be remarked that P. nigra is the main host of PPM in this area (Sangüesa-Barreda et al., 2014). Expert forest technicians from the Aragón “Forest Health Laboratory” (Mora de Rubielos, Spain) have been assessing the degree of PPM defoliation annually from 1970 to 2011 in this area [see more details on this study area in Montoya (1970), Hernández Alonso et al. (2005) and Sangüesa-Barreda et al. (2014)]. PPM defoliation was recorded in a relative scale from 0 (absence of defoliation) to 5 (severe defoliation) during the late winter season (from December to February) following Montoya and Hernández (1991). PPM defoliation was estimated using a PPM damage visual scale at the stand level. First, a carefully inspection and delimitation of the stand was performed in the field. Second, the degree of defoliation was estimated according to the following scale: (0) absence of nests and defoliations; (1) undamaged trees with some PPM nests; (2) light defoliations affecting stand edges and isolated trees; (3) defoliations affecting stand edges and presence of nests in the center of the stand; (4) severe defoliations in the edges and abundance of nests in the center; (5) extreme and partial defoliations in the edge and center of the stands respectively. This monitoring was carried out in 92 forest stands distributed in the study region (Figure S2; Figure 1). In addition, for each forest stand they recorded information on the altitude, slope, and structure variables (density and basal area) (Table 1). The density and basal area of the stand were estimated at the beginning of the study period and they represent only one value per stand over the whole study period.

FIGURE 1

Figure 1. Views of (A) pine processionary moth (PPM) defoliation levels 3, 4, and 5 used in the monitoring of the Mora de Rubielos study area (Teruel, Spain), (B) illustrative winters with different PPM defoliation levels (0–5, color scales). In the plot (B) the inset shows the location of the study area (box) and the Circum-mediterranean countries where T. pityocampa has been reported are colored (orange areas).

TABLE 1

Table 1. Topographic, structural, and defoliation data measured in the monitored pine stands affected by the pine processionary moth in Mora de Rubielos, Teruel, Spain.

In the study area the climate is Mediterranean continental with mean temperatures of 9.4°C and total precipitation of 728 mm at 1492 m a.s.l. (Figure S3). Climate data for the study area (mean monthly and seasonal minimum and maximum temperatures, total monthly and seasonal precipitations) were obtained by calculating regional climate series based on local stations situated at different altitudes (see Sangüesa-Barreda et al., 2014).

Data Analyses

Spatiotemporal Trends in PPM Defoliation and Their Determinants

To study spatiotemporal trends in PPM defoliation (PPMD) we used Generalized Additive Mixed Models (GAMM; Lin and Zhang, 1999; Wood, 2006). GAMM is a flexible semi parametric method that allows modeling complex non-linear relationships among a response variable and one or several covariates (Wood, 2004). GAMM allows handling the space-time trend and interaction in the response variable (Augustin et al., 2009). In our case, the temporal trend in PPMD might vary across the study region spatially due to variations in altitude as well as forest structure and composition since several pine species with different stand densities are spread all over the study area. Further, PPMD might vary temporally because of the variation in the climate conditions.

We considered as potential covariates for the GAMM: the year of defoliation, the altitude, tree density of each stand, the relative abundance of P. nigra in each stand and the average minimum winter temperature. We only included the abundance of P. nigra in the models because: (i) the distribution of tree species changed along the altitudinal gradient (Figure S4), with P. pinaster and P. halepensis dominating low altitudes and P. sylvestris dominant at high altitudes, whilst P. nigra was less correlated with elevation; (ii) defoliation in P. nigra stands was higher than defoliation in the other stands (Table 1). The distribution of P. nigra was positively related with the distribution of P. sylvestris (Spearman ρ = 0.27; P < 0.05) and negatively with that of P. pinaster (ρ = 0.48; P < 0.05). As climate variables we included the average minimum winter temperature because low temperatures can damage PPM larvae (Battisti et al., 2015). To test if monthly or climate data presented significant trends, we used the Kendall tau statistic (τ).

We proposed a GAMM with PPMD as response variable. Since the response variable represent ordered values in discrete classes (semi-quantitative variable), we used a Poisson response distribution with a log link (Wood, 2006). In addition, site (92 levels) was included as random factor to account for variations from site to site, and a first-order autocorrelation structure [AR(1)] was used to remove temporal autocorrelation in the model residuals due to the dependency of defoliation in year t for that in year t-1. Due to the use of the Poisson family, model comparisons based on log-likelihood parameters can be misleading (Wood, 2006). For that reason, we performed model selection in two steps. First, we created a single model including a cubic regression spline for calendar year as the only response variable (i.e., the temporal trend). After that, the rest of spatial covariates were included in the model using thin-plate regression splines. All covariates were retained because they had significant influences on PPMD (see Table 2). Second, the spatiotemporal trend in PPMD was modeled as a function of: calendar year, altitude, the relative abundance of P. nigra and their interaction. Specifically, the spatiotemporal trend is represented by a two-dimensional tensor product smooth including two-dimensional thin plate splines for space allowed to interact with a temporal trend represented by a cubic regression spline for calendar year (Augustin et al., 2009). Thus, the temporal trend in PPMD can vary along the altitudinal gradient and with the variation in the relative abundance of P. nigra. Further, tree density in each stand and the winter mean minimum temperature were retained in the model as thin-plate regression splines.

TABLE 2

Table 2. List of models selected to explain spatiotemporal trends in pine processionary moth defoliation at a regional scale.

Synchrony in PPMD Among Stands

A Mantel correlogram was used to assess if there was distance-decay in PPMD similarity (i.e., to test if sites showed more similar PPMD values as a function of distance than expected by random). We computed a defoliation distance matrix among sites using the Jaccard distance (Legendre and Legendre, 2012). We chose this index because it is suitable for variables recorded in ordinal scales. A matrix of geographical distance was computed using the geographical coordinates of each stand (n = 92). The Spearman rank correlation function was used to test for the change in PPMD distance with the change in geographical distance. The significance of each distance class was tested using 999 unrestricted permutations.

All analyses were performed in the R statistical environment (R Development Core Team, 2018). The package “mgcv” (Wood, 2006) was used to perform the GAMM, and the package “vegan” (Oksanen et al., 2019) was used to calculate the Mantel Correlogram.

Results

Spatiotemporal Trends in PPMD and Their Determinants

We observed strong spatiotemporal patterns in PPM defoliation in the study site (Figures 1B, 2). The temporal trend showed marked peaks in PPM defoliation (assessed in late winter), one from 1979 to 1983 and the others in two consecutive years in 1991–1992 and 1995–1996 (Figure 2). Low defoliation values were observed in 1984–1985 and 1994–1995. In the study area, warming trends have prevailed in spring and summer, particularly in March, whilst in winter only maximum temperatures showed a significant (P < 0.05) increase (Table S1). No significant trend was observed for precipitation data.

FIGURE 2

Figure 2. Observed and predicted values of pine processionary moth (PPM) defoliation (0–5) at the regional scale (data from Mora de Rubielos study area, Teruel, Spain). The black line and dots represent the observed values of PPM defoliation at regional scale (mean across the 92 sites) and the vertical bars indicate the standard errors. The green circles and line indicate the predicted PPM defoliation according to Generalized Additive Mixed Models including the covariates presented in Table 2.

The results of the different GAMMs supported the existence of a strong temporal pattern in PPMD (Figure 2, Table 2). The results also indicated that part of the spatial variation in PPMD was attributable to variations in P. nigra relative abundance and tree density, with an increase in PPMD with the increase in these two factors (Table 2, Figure 3). The model also indicated that minimum winter temperatures influenced PPMD with values below −12°C lowering the incidence of PPMD (Figure S5; Table 2). The relationship between altitude and PPMD was more complex and suggested the existence of higher defoliations in P. nigra stands located at intermediate altitudes (from 1300 to 1500 m; Table 1, Figure 3; Figure S6). Finally, when a spatiotemporal interaction between P. nigra relative abundance, elevation and year was included the GAMM accounted for 44% of the variation in PPMD (Table 2). There was a significant three-way interaction between tree density, altitude and calendar year (i.e., spatiotemporal trend). Overall, PPMD increased in P. nigra stands located at intermediate altitudes and in response to warmer winter conditions but decreased with increasing tree density.

FIGURE 3

Figure 3. Temporal changes in pine processionary moth (PPMD) as a function of altitude (y axes, in m) in P. nigra stands and in stands dominated by other pine species according to the results of the GAMM. The PPMD varies from 0 to 5.

Synchrony in PPMD Among Stands

We found that the synchrony in defoliation across sites decreased with distance (Figure 4). The Mantel correlogram showed the existence of stands with similar PPMD up to a distance of 6 km. At distances from 6 to 9 km the Mantel correlation coefficients were negative (Figure 4).

FIGURE 4

Figure 4. Mantel correlogram of pine processionary moth defoliation data. The x axis shows distance between sites. Filled symbols indicate significant (P < 0.05) Mantel correlation coefficients.

Discussion

Spatiotemporal Trends in PPMD and Its Potential Expansion to Higher Elevations

The main aim of our study was to model the spatiotemporal dynamics of PPM defoliation and to assess its main drivers paying special attention to the potential expansion of PPM to higher elevations. In this sense, our spatiotemporal model accounted for more than 40% of the variation in PPM defoliation and suggested the influence of some factors already evaluated by other authors. However, we find no support for the increase in elevation of PPMD, because temporal patterns of defoliation where similar at different altitudes (Figure S7). The most feasible explanation of this lack of upward PPM expansion could be that winter temperatures, which are fundamental to allow larvae feeding (Battisti et al., 2005; Buffo et al., 2007; Robinet et al., 2014), have not increased significantly during the studied period (Table S1). Another explanation could be the change in tree species dominance along the altitudinal gradient if there is a trend toward a lower abundance of the preferred host, P. nigra, at high altitudes preventing PPM expansion (Sangüesa-Barreda et al., 2014). Finally, the potential occurrence of lethal temperatures during winter may reduce larvae survival thus preventing the altitudinal expansion of the species as suggested by the negative influence of winter minimum temperatures below −12°C on PPMD.

In agreement with previous studies we found a strong impact of the host pine species on PPMD and outbreaks (Hódar et al., 2012), being P. nigra the species with highest defoliation rates, particularly at intermediate altitudes (Table 1). Probably, intermediate altitudes present enough relative humidity and more favorable temperatures than high altitudes thus favoring PPM (Robinet et al., 2014). Interestingly, we found a negative relationship between stand density and PPMD which indicates that for a comparable number of insects the probability of a tree being affected by PPM is higher in open than in closed stands. A similar result was found in pure maritime pine plantations in southern France (Régolini et al., 2014). This negative relationship can be also due the higher occurrence of PPM outbreaks in edge forests than in the core of stands (Samalens and Rossi, 2011; Dulaurent et al., 2012). Thus, our results align with previous findings indicating that both stand density and composition are important factors determining PPMD.

We also found a significant influence of winter minimum temperatures on PPMD, with the occurrence of temperatures below −12°C strongly reducing the occurrence of severe PPMD (Battisti et al., 2005; Buffo et al., 2007). This agrees with empirical studies reporting lower lethal temperatures for the survival of PPM larvae of −16°C (Démolin, 1969; Hoch et al., 2009). In addition, winter temperatures below −7°C can also damage PPM nests (Battisti et al., 2015). Allstadt et al. (2015) suggested that weather oscillations can influence gypsy moth (Lymantria dispar) defoliation outbreaks directly through effects on mortality or survivorship, but they also indicated indirect trophic influences on this insect related to masting or predator (e.g., mice) abundance depending on oak production. It has been suggested that PPMD increases as winter temperatures warm with milder conditions favoring the occurrence of PPM outbreaks and very cold temperatures, below the thermal limit of tolerance, increasing PPM larvae mortality (Battisti et al., 2005). We found a significant positive effect of winter temperature on PPMD, which is in line with these previous findings, but our analyses also suggest other influencing factors as half of the variation in PPMD remains unexplained.

Cyclicity and Synchrony in PPM Defoliation

Previous studies have reported either oscillatory patterns or cyclic behavior of PPM defoliations and attributed them to climate forcing and density-dependent population dynamics (Buffo et al., 2007; Hódar et al., 2012; Li et al., 2015). Contrastingly, other authors have reported lack of PPM cyclicity (Tamburini et al., 2013). First, differences appear in the studied regions ranging from a entire country (e.g., Li et al., 2015) to local regions (e.g., Hódar et al., 2012). Second, discrepancies between studies can also arise because of the different methods used to assess PPMD and the diverse statistical methods applied, but also due to the short length of the studied periods. In fact, we used a 42-year long record of PPMD which represents one of the longest and most complete time series studied for this species to the best of our knowledge. Our results align with the second group of studies as we found no evidence for the existence of cyclic PPM defoliations. However, these findings should be interpreted with caution due to: (i) the method applied to assess PPM incidence based on a relative scale of defoliation recorded at the stand level, and (ii) local factors modifying the regional-scale influence of climate and other large-scale forcings. These local factors include the structure and composition of affected forests since mixed or natural pine forests may be less susceptible to PPM defoliation than pure forests and pine plantations (Cayuela et al., 2011). In addition, the phase of PPM population dynamics may also affect cyclicity because PPM populations in the positive gradation phase, before defoliation incidence peaks, may be more prone to reflect large-scale climatic influences, and thus cyclicity, than PPM populations in the negative gradation phase whose density is probably controlled by food availability or natural enemies (Klemola et al., 2010). Future studies should consider the possibility of reconstructing PPM defoliations using tree-ring data in order to decipher if synchronous defoliations have occurred in the past at higher or lower altitudes (Myers and Cory, 2013).

PPMD and related outbreaks may also depend on global forcings such as sunspot numbers (Galkin, 1992; Selås, 2014) and winter atmospheric patterns as those characterized by the North Atlantic Oscillation (Hódar et al., 2012), but PPM activity is also strongly influenced by factors operating at local scales such as climate conditions, food availability or the density of predators and parasites (Battisti et al., 2005; Li et al., 2015). The occurrence of severe frost damages in one region may strongly influence larvae survival and thus disrupt the cyclicity in PPM outbreaks (Robinet et al., 2014). The influence of local factors on regional PPMD patterns may reduce the similarity in defoliation with the increase in distance and over large areas of its distribution (Myers and Cory, 2013; Li et al., 2015). Thus, synchrony depends on different factors such as the distance between the insect populations and the relative importance of external factors (Myers and Cory, 2013; Allstadt et al., 2015). It is expected that synchrony will decrease with distance if the main cause of the synchrony is the dispersal between populations (Myers and Cory, 2013). However, if the synchrony is due to external forcing then the relative influence of those forcing factors and their spatial variation will determine the synchrony (Selås, 2014; Allstadt et al., 2015). For example, Li et al. (2015) found that synchrony between PPM defoliation across France was clustered in regions and attributed this clustering to regional climate, specifically to the influence of winter temperatures. Thus, some PPM populations could be more sensitive than others to climate due to the relative influence of local factors including tree density, stand composition or the PPM gradation phase. In our study, we found that synchrony between stands in PPMD was significant until distances up to 6 km (Figure 4), which is close to the mean radius of the area dominated by P. nigra forests. These results can be due to the change in the tree species dominating the stand at distances higher than 6 km or due to limited PPM dispersal between stands at larger distances. If PPM moths can disperse large distances (Régolini et al., 2014; Battisti et al., 2015), it is more plausible to think that changes in stand density and composition can create the observed distance-decay in PPMD similarity.

Conclusions

Our results do not support the main hypothesis that PPM defoliation intensity have increased upwards because of climate warming. This lack of an upward shift of PPM defoliation could be explained because winter minimum temperatures have not significantly increased in the study area, whereas maximum temperatures did. Along this, our results indicate that PPM defoliation dynamics are dependent on minimum winter temperature, which points to the importance of lethal temperatures below a threshold (in our case −12°C) which may compromise PPM survival. However, our results also suggest that spatial covariates related with the distribution of different tree species and tree density play an important role on explaining PPMD spatiotemporal patterns. Thus, our results align with previous studies which have identified that forest composition and structure might be important factors determining the potential expansion of PPM (Samalens and Rossi, 2011; Castagneyrol et al., 2014; Régolini et al., 2014). All in all, the main findings of this study indicate that considering both, the spatial pattern and the temporal trend of defoliation, as well as their interactions, allow a better understanding of the factors determining PPM dynamics.

Data Availability Statement

The datasets generated for this study are available on request to the corresponding author.

Author Contributions

RH-A organized the data collection. AG analyzed the data. AG and JC lead the writing of the manuscript. All authors contributed to discussing and revising the manuscript.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We acknowledge funding by the Spanish Ministry of Economy Fundiver project (CGL2015-69186-C2-1-R). We thank the Spanish Environment Ministry and Aragón Government who supported PPM monitoring programs at Mora de Rubielos (Teruel). We acknowledge support of the publication fee by the CSIC Open Access Publication Support Initiative through its Unit of Information Resources for Research (URICI).

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo.2019.00458/full#supplementary-material

References

Allstadt, A. J., Liebhold, A. M., Johnson, D. M., Davis, R. E., and Haynes, K. L. (2015). Temporal variation in the synchrony of weather and its consequences for spatiotemporal population dynamics. Ecology 96, 2935–2946. doi: 10.1890/14-1497.1

PubMed Abstract | CrossRef Full Text | Google Scholar

Augustin, N. H., Musio, M., von Wilpert, K., Kublin, E., Wood, S. N., and Schumacher, M. (2009). Modelling spatiotemporal forest health monitoring data. J. Am. Statist. Ass. 104, 899–911. doi: 10.1198/jasa.2009.ap07058

CrossRef Full Text | Google Scholar

Battisti, A., Avci, M., Avtzis, D. N., Jamaa, M. L. B., Berardi, L., Berretima, W., et al. (2015). “Natural history of the processionary moths (Thaumetopoea spp.): new insights in relation to climate change,” in Processionary Moths and Climate Change: An Update, ed A. Roques (Dordrecht: Springer), 15–79. doi: 10.1007/978-94-017-9340-7_2

CrossRef Full Text | Google Scholar

Battisti, A., Stastny, M., Netherer, S., Robinet, C., Schopf, A., Roques, A., et al. (2005). Expansion of geographic range in the pine processionary moth caused by increased winter temperatures. Ecol. Appl. 15, 2084–2096. doi: 10.1890/04-1903

CrossRef Full Text | Google Scholar

Bjørnstad, O. N., Robinet, C., and Liebhold, A. M. (2010). Geographic variation in North American gypsy moth cycles: subharmonics, generalist predators, and spatial coupling. Ecology 91, 106–118. doi: 10.1890/08-1246.1

PubMed Abstract | CrossRef Full Text | Google Scholar

Buffo, E., Battisti, A., Stastny, M., and Larsson, S. (2007). Temperature as a predictor of survival of the pine processionary moth in the Italian Alps. Agric. For. Entomol. 9, 65–72. doi: 10.1111/j.1461-9563.2006.00321.x

CrossRef Full Text | Google Scholar

Camarero, J. J., Gazol, A., Tardif, J. C., and Conciatori, F. (2015). Attributing forest responses to global-change drivers: limited evidence of a CO₂-fertilization effect in Iberian pine growth. J. Biogeogr. 42, 2220–2233. doi: 10.1111/jbi.12590

CrossRef Full Text | Google Scholar

Castagneyrol, B., Régolini, M., and Jactel, H. (2014). Tree species composition rather than diversity triggers associational resistance to the pine processionary moth. Basic Appl. Ecol. 15, 516–523. doi: 10.1016/j.baae.2014.06.008

CrossRef Full Text | Google Scholar

Cayuela, L., Hernández, R., Hódar, J. A., Sanchez, G., and Zamora, R. (2014). Tree damage and population density relationships for the pine processionary moth: prospects for ecological research and pest management. For. Ecol. Manage. 328, 319–325. doi: 10.1016/j.foreco.2014.05.051

CrossRef Full Text | Google Scholar

Cayuela, L., Hódar, J. A., and Zamora, R. (2011). Is insecticide spraying a viable and cost-efficient management practice to control pine processionary moth in mediterranean woodlands? For. Ecol. Manage. 261, 1732–1737. doi: 10.1016/j.foreco.2011.01.022

CrossRef Full Text | Google Scholar

Démolin, G. (1969). Bioecología de la “procesionaria del pino” Thaumetopoea pityocampa Schiff. Incidencia de los factores climáticos. Bol. Serv. Plagas For. 23, 9–24.

Google Scholar

Dulaurent, A.-M., Porté, A. J., van Halder, I., Vétillard, F., Menassieu, P., and Jactel, H. (2012). Hide and seek in forests: colonization by the pine processionary moth is impeded by the presence of nonhost trees. Agric. For. Entomol. 14, 19–27. doi: 10.1111/j.1461-9563.2011.00549.x

CrossRef Full Text | Google Scholar

Galkin, G. I. (1992). [Periodicity of mass multiplication of pests.] Zashch. Rast. 12, 37–39.

Google Scholar

Hernández Alonso, R., Pérez Fortea, V., Camarero, J. J., Montoya, R., and Sánchez Peña, G. (2005). “Efectos de la defoliación inducida por la procesionaria del pino (Thaumetopoea pityocampa) sobre el crecimiento y la supervivencia de Pinus nigra e interacciones con el clima durante el periodo 1992–2004,” in Actas del IV Congreso Forestal Español (Zaragoza), 366–374.

Google Scholar

Hoch, G., Toffolo, E. P., Netherer, S., Battisti, A., and Schopf, A. (2009). Survival at low temperature of larvae of the pine processionary moth Thaumetopoea pityocampa from an area of range expansion. Agric. For. Entomol. 11, 313–320. doi: 10.1111/j.1461-9563.2009.00431.x

CrossRef Full Text | Google Scholar

Hódar, J. A., Castro, J., and Zamora, R. (2003). Pine processionary caterpillar Thaumetopoea pityocampa as a new threat for relict Mediterranean Scots pine forests under climatic warming. Biol. Conserv. 110, 123–129. doi: 10.1016/S0006-3207(02)00183-0

CrossRef Full Text | Google Scholar

Hódar, J. A., Cayuela, L., and Zamora, R. (2012). Climate change and the incidence of a forest pest in Mediterranean ecosystems, can the North Atlantic Oscillation be used as a predictor? Clim. Change 113, 699–711. doi: 10.1007/s10584-011-0371-7

CrossRef Full Text | Google Scholar

Hódar, J. A., and Zamora, R. (2004). Herbivory and climatic warming, a mediterranean outbreaking caterpillar attacks a relict, boreal pine species. Biodivers. Conserv. 13, 493–500. doi: 10.1023/B:BIOC.0000009495.95589.a7

CrossRef Full Text | Google Scholar

Hódar, J. A., Zamora, R., Castro, J., and Baraza, E. (2004). Feast and famine: previous defoliation limiting survival of pine processionary caterpillar Thaumetopoea pityocampa in Scots pine Pinus sylvestris. Acta Oecol. 26, 203–210. doi: 10.1016/j.actao.2004.05.004

CrossRef Full Text | Google Scholar

Jacquet, J. S., Bosc, A., O'Grady, A. P., and Jactel, H. (2013). Pine growth response to processionary moth defoliation across a 40-year chronosequence. For. Ecol. Manage. 293, 29–38. doi: 10.1016/j.foreco.2012.12.003

CrossRef Full Text | Google Scholar

Jactel, H., Menassieu, P., Vétillard, F., Barthélémy, B., Piou, D., Frérot, B., et al. (2006). Population monitoring of the pine processionary moth (Lepidoptera: Thaumetopoeidae) with pheromone-baited traps. For. Ecol. Manage. 235, 96–106. doi: 10.1016/j.foreco.2006.08.002

CrossRef Full Text | Google Scholar

Klemola, N., Andersson, T., Ruohomäki, K., and Klemola, T. (2010). Experimental test of parasitism hypothesis for population cycles of a forest lepidopteran. Ecology 91, 2506–2513. doi: 10.1890/09-2076.1

PubMed Abstract | CrossRef Full Text | Google Scholar

La Porta, N., Capretti, P., Thomsen, I. M., Kasanen, R., Hietala, A. M., and Von Weissenberg, K. (2008). Forest pathogens with higher damage potential due to climate change in Europe. Can. J. Plant. Pathol. 30, 177–195. doi: 10.1080/07060661.2008.10540534

CrossRef Full Text | Google Scholar

Legendre, P., and Legendre, L. (2012). Numerical Ecology, 3rd Edn. Amsterdam: Elsevier.

Google Scholar

Li, S., Daudin, J. J., Piou, D., Robinet, C., and Jactel, H. (2015). Periodicity and synchrony of pine processionary moth outbreaks in France. For. Ecol. Manage. 354, 309–317. doi: 10.1016/j.foreco.2015.05.023

CrossRef Full Text | Google Scholar

Liebhold, A., Elkinton, J., Williams, D., and Muzika, R.-M. (2000). What causes outbreaks of the gypsy moth in North America? Popul. Ecol. 42, 257–266. doi: 10.1007/PL00012004

CrossRef Full Text | Google Scholar

Liebhold, A. M., Plymale, R., Elkinton, J. S., and Hajek, A. E. (2013). Emergent fungal entomopathogen does not alter density dependence in a viral competitor. Ecology 94, 1217–1222. doi: 10.1890/12-1329.1

CrossRef Full Text | Google Scholar

Lin, X., and Zhang, D. (1999). Inference in generalized additive mixed models by using smoothing splines. J. R. Stat. Soc. Series B Stat. Methodol. 61, 381–400. doi: 10.1111/1467-9868.00183

CrossRef Full Text | Google Scholar

Montoya, R. (1970). Zona de mora de rubielos. descripción, programas y trabajos efectuados en 1970. Bol. Serv. Plag. Forest 26, 119–129.

Google Scholar

Montoya, R., and Hernández, R. (1991). “La procesionaria del pino,” in Plagas de Insectos en las Masas Forestales Españolas, eds N. Romanyk and D. Cadahia (Madrid: Ministerio de Agricultura, Pesca y Alimentación, 59–73.

Google Scholar

Myers, J. H. (1993). Population outbreaks in forest Lepidoptera. Am. Sci. 81, 240–251.

Google Scholar

Myers, J. H., and Cory, J. S. (2013). Population cycles in forest Lepidoptera revisited. Annu. Rev. Ecol. Evol. 44, 565–592. doi: 10.1146/annurev-ecolsys-110512-135858

CrossRef Full Text | Google Scholar

Oksanen, J., Blanchet, F. G., Friendly, M., Kindt, R., Legendre, P., McGlinn, D., et al. (2019). Vegan: Community Ecology Package. R package version 2.5-4. Available online at: https://CRAN.R-project.org/package=vegan

Google Scholar

Palacio, S., Hernández, R., Maestro-Martínez, M., and Camarero, J. J. (2012). Fast replenishment of initial carbon stores after defoliation by the pine processionary moth and its relationship to the re-growth ability of trees. Trees 26, 1627–1640. doi: 10.1007/s00468-012-0739-y

CrossRef Full Text | Google Scholar

Pautasso, M., Schlegel, M., and Holdenrieder, O. (2015). Forest health in a changing world. Microb. Ecol. 69, 826–842. doi: 10.1007/s00248-014-0545-8

PubMed Abstract | CrossRef Full Text | Google Scholar

R Development Core Team (2018). R: A Language and Environment for Statistical Computing. Vienna: R Foundation for Statistical Computing.

PubMed Abstract | Google Scholar

Régolini, M., Castagneyrol, B., Dulaurent-Mercadal, A.-M., Piou, D., Samalens, J.-C., and Jactel, H. (2014). Effect of host tree density and apparency on the probability of attack by the pine processionary moth. For. Ecol. Manage. 334, 185–192. doi: 10.1016/j.foreco.2014.08.038

CrossRef Full Text | Google Scholar

Robinet, C., Baier, P., Pennerstorfer, J., Schopf, J., and Roques, A. (2007). Modelling the effects of climate change on the pine processionary moth (Thaumetopoea pityocampa L.) expansion in France. Glob. Ecol. Biogeogr. 16, 460–471. doi: 10.1111/j.1466-8238.2006.00302.x

CrossRef Full Text | Google Scholar

Robinet, C., Imbert, C. E., Rousselet, J., Sauvard, D., Garcia, J., Goussard, F., et al. (2012). Human-mediated long-distance jumps of the pine processionary moth in Europe. Biol. Invasions 14, 1557–1569. doi: 10.1007/s10530-011-9979-9

CrossRef Full Text | Google Scholar

Robinet, C., Rousselet, J., and Roques, A. (2014). Potential spread of the pine processionary moth in France: preliminary results from a simulation model and future challenges. Ann. For. Sci. 71, 149–160. doi: 10.1007/s13595-013-0287-7

CrossRef Full Text | Google Scholar

Roques, A., Rousselet, J., Avc, M., Avtzis, D., Basso, A., Battisti, A., et al. (2015). “Climate warming and past and present distribution of the processionary moths (Thaumetopoea) in Europe, Asia Minor and North Africa,” in Processionary Moths and Climate Change: An Update, ed A. Roques (Dordrecht: Springer, 81–161.

Google Scholar

Samalens, J.-C., and Rossi, J.-P. (2011). Does landscape composition alter the spatiotemporal distribution of the pine processionary moth in a pine plantation forest? Popul. Ecol. 53, 287–296. doi: 10.1007/s10144-010-0227-4

CrossRef Full Text | Google Scholar

Sánchez-Bayo, F., and Wyckhuys, K. A. (2019). Worldwide decline of the entomofauna: a review of its drivers. Biol. Conserv. 232, 8–27. doi: 10.1016/j.biocon.2019.01.020

CrossRef Full Text | Google Scholar

Sangüesa-Barreda, G., Camarero, J. J., García-Martín, A., Hernández, R., and de la Riva, J. (2014). Remote-sensing and tree-ring based characterization of forest defoliation and growth loss due to the Mediterranean pine processionary moth. For. Ecol. Manage. 320, 171–181. doi: 10.1016/j.foreco.2014.03.008

CrossRef Full Text | Google Scholar

Selås, V. (2014). Linking ‘10-year' herbivore cycles to the lunisolar oscillation: the cosmic ray hypothesis. Oikos 123, 194–202. doi: 10.1111/j.1600-0706.2013.00716.x

CrossRef Full Text | Google Scholar

Sturrock, R. N., Frankel, S. J., Brown, A. V., Hennon, P. E., Kliejunas, J. T., Lewis, K. J., et al. (2011). Climate change and forest diseases. Plant Pathol. 60, 133–149. doi: 10.1111/j.1365-3059.2010.02406.x

CrossRef Full Text | Google Scholar

Tamburini, G., Marini, L., Hellrigl, K., Salvadori, C., and Battisti, A. (2013). Effects of climate and density-dependent factorson population dynamics of the pine processionary moth in the Southern Alps. Clim. Change 121, 701–712. doi: 10.1007/s10584-013-0966-2

CrossRef Full Text | Google Scholar

Toïgo, M., Barraquand, F., Barnagaud, J.-Y., Piou, D., and Jactel, H. (2017). Geographical variation in climatic drivers of the pine processionary moth population dynamics. For. Ecol. Manage. 404, 141–155. doi: 10.1016/j.foreco.2017.08.024

CrossRef Full Text | Google Scholar

Wood, S. (2004). Stable and efficient multiple smoothing parameter estimation for generalized additive models. J. Am. Statist. Ass. 99, 673–686. doi: 10.1198/016214504000000980

CrossRef Full Text | Google Scholar

Wood, S. (2006). Generalized Additive Models: An Introduction With R. Boca Raton, FL: Chapman and Hall/CRC Press. doi: 10.1201/9781420010404

CrossRef Full Text | Google Scholar