Direct and Indirect Analysis of the Elevational Shift of Larch Budmoth Outbreaks Along an Elevation Gradient

Larch budmoth (LBM) periodically defoliates alpine stands of European larch during vast outbreaks occurring generally at 8–10 year intervals. LBM outbreaks recently declined and the ongoing global warming has been pointed out as a possible cause of this decline. In this article, we reconstructed the recent history of LBM outbreaks at different elevations along a larch elevational gradient in the French Alps using direct and indirect observations based on tree-ring width and density analysis, and compared it with local climatic data. We found that LBM outbreaks time-series were better reconstructed with latewood density than with ring width. We also found that there was a recent but limited elevational shift of LBM outbreaks from medium toward higher elevations. We suggest that this elevational shift is a consequence of the variable effect of the global warming at the different elevations. Winter warming is expected to affect differently the timing of LBM egg hatch as well as that of larch bud flush, larvae being at present susceptible to emerge whereas no needles are available as food at the former optimal elevation. A better synchronization between larch and LBM may exist at higher elevations.


INTRODUCTION
Larvae of larch budmoth (Zeiraphera griseana; LBM) periodically defoliate alpine stands of European larch (Larix decidua) during massive outbreaks occurring most of the time at 8-10 year intervals, with an optimum at 1800 m elevation (Dormont et al., 2006;Baltensweiler et al., 2008). The defoliation decreases photosynthesis rate, lowers the production of carbohydrates and reduces annual ring growth (Dormont et al., 2006;Baltensweiler et al., 2008). Larch is able to refoliate 3-4 weeks after defoliation. The newly produced needles are shorter, as well as the needles produced in the following spring (Baltensweiler et al., 2008). Such radial growth reduction leaves a characteristic fingerprint in annual ring time-series used to retrospectively reconstruct LBM outbreaks (Weber, 1997;Rolland et al., 2001;Esper et al., 2007;Baltensweiler et al., 2008;Saulnier et al., 2017). For example, a long-term dendrochronology study of larch trees in the Alps have discovered 123 outbreaks between 781 A.D. and 1981 (Iyengar et al., 2016). Another one used and compared different methods to ultimately identify 15 outbreaks between 1774 and 1964 in the Central Italian Alps (Cerrato et al., 2019). Additionally, tree-rings analysis showed spatial heterogeneity of LBM events, with a 1-year delay between intermediate Alps, and Northern and Southern Alps outbreaks, and a decrease of the strength of the periodic signal from the epicenter to the outside of the distribution (Saulnier et al., 2017). Analysis of extensive records also showed elevational fluctuations of LBM outbreaks that coincide with temperature shifts (Johnson et al., 2010). LBM outbreaks recently weakened and the ongoing global warming has been pointed out as a possible cause of this weakening (Battipaglia et al., 2014;Saulnier et al., 2017). The warming is predicted to be maximum in mountain areas at medium to high elevations (Kotlarski et al., 2015). In some parts of the Alps, temperature already increased twice as fast as in the northern hemisphere during the previous century (Rebetez and Reinhard, 2008). This warming rate is superior at high elevation (Mountain Research Initiative Edw Working Group et al., 2015). In this article, we reconstructed the recent history of LBM outbreaks from 1960 to 2007 at different elevations between 1200 and 2400 m.a.s.l. along a larch elevational gradient in the French Alps, and compared it with local climatic data.
Annual ring reconstruction of LBM outbreaks show variable results according to the method of identification of the outbreak rings (Büntgen et al., 2009;Cerrato et al., 2019). In this article, we combined annual-ring time series with direct in situ observations in order to reconstruct consistently the history of LBM outbreaks for similar periods between 1960 and 2016 ( Table 1). In one forest stand located in the Briançon area of the Southeastern French Alps, we compared the temporal coincidence between the annual variations in LBM population density estimated by branch samplings, and annual-ring time series measured on the L. decidua host trees. We compared the results obtained with different annual-ring variables. We used these relationships to reconstruct the recent history of LBM outbreaks at different elevations along a larch elevational gradient located a few km part. We observed the variation with time of the LBM outbreaks at the different elevations and we discussed the results in the context of the local warming.

MATERIALS AND METHODS
We collected two types of time-series: time-series of direct observations of LBM outbreaks and time-series of indirect LBM outbreaks based on tree-ring analysis (LBM time-series and Treering index time-series, respectively, Table 1). Table 1 shows the location of the study site and of the direct and indirect LBM observations.

Direct Observations of the LBM Outbreaks
The "Les Combes" site was selected for surveying LBM population dynamics because is located on a sunny, south-facing slope, and the local population cycle therefore expected to be advanced by one year compared to other populations in the nearby valleys such as those of the north-oriented forests of Villard Saint Pancrace (Figure 1; Roques and Goussard, 1982). LBM density at the Les Combes site was measured since 1960 using two successive methods. From 1960 to 1978, 10 trees were randomly sampled in the site from early to mid-June, the sampling date being fixed each year with regard to the appearance of the first larval sheaths typical of 4th-instar larvae. At this moment, all LBM eggs were considered to have hatched (Auer, 1961). The LBM larvae of the different instars were recorded on the collected branches in order to measure the mean density of larvae per kilogram of branch per tree and per site (Auer et al., 1981). From 1977 on, we used a simpler method where larval density was estimated as the number of larvae per meter of branch (Roques and Goussard, 1982). When the 4th-instar LBM larvae were observed, five 40 cm-long branches were randomly cut per tree at mid-crown on 10 trees selected according to a zig-zag path covering the entire site. The mean density per meter of branch was then calculated at the site level. In order to compare the two periods, we used a conversion value of 44.2 m of branch for 1 kg (Roques and Goussard, 1982).

Tree-Ring Analysis
We collected the increment cores at breast-height with powerdriven increment borers, following a perpendicular-to-the slope orientation to avoid compression wood. The trees were all adult, dominant trees with no visible major defect. The freshly collected increment cores were stored in polycarbonate honeycomb boxes, brought back to the laboratory, where they were oven-dried, and sawn to a 2 mm uniform thickness. Then they were resinextracted with a 48 h pentane-bath. They were dried again and then X-rayed using the indirect procedure invented in 1966 by Polge (1978) and continuously improved since then (Polge, 1978;Mothe et al., 1998). Finally the X-ray films were scanned at 4000 dpi in order to acquire the microdensity profiles using the WINDENDRO software (Windendro 2008e Regent instruments Canada inc. 1 ).
The microdensity profiles were visually checked and compared to the sawn samples. Some profiles were discarded because ring boundaries were blurred by a very strong grain angle.

Annual Ring Variables
The raw microdensity profiles were imported into R (R Core Team, 2018) and the following conventional annual ring variables were calculated using specially designed R functions: ring width, earlywood width, latewood width, earlywood density, latewood density and standard deviation of ring microdensity (ring density standard deviation). Earlywood-latewood boundary was calculated using the conventional extreme average method, based on the mean of the extreme ring density values (Rozenberg et al., 2002).

Climatic Data
Climatic data are daily values of minimum and maximum temperature and precipitations coming from the closest Météo-France weather station (French national weather service). This weather station was first located in Briançon (N44.905856, W6.633046, 1306 m.a.s.l.) and started to record data in 1966. In 2003, a new station was created and started to record climatic data

Data Analysis
Our reference is the longest-term time series of LBM outbreaks  obtained at one site, Les Combes, elevation 1825 m.a.s.l. To study the time variation of LBM outbreaks at different elevations from 1200 to 2300 m.a.s.l., we used fingerprints left by the periodic outbreaks on tree-ring data collected at different elevations. The tree-ring time-series were cross-dated using the software INTERDAT.exe (ver. 1.1) (Becker et al., 1994;Mérian et al., 2011) and pointer-years, then detrended. We used the residuals RCS (Regional Curve Standardization) method to eliminate the lowfrequency aging effect related with cambial age (Cook and Peters, 1997;Esper et al., 2003). To separate the growth and density reductions caused by LBM outbreaks from those caused by climatic or other biotic factors, we used the moving time-window approach proposed by Büntgen et al. (2009). This approach was tested and validated by comparing the direct observations of the LBM outbreaks with the tree-ring time-series of the 10 larch trees at the Les Combes site where we observed the LBM outbreaks time-series. For each time-window and each tree, the moving time-window approach detects and records the year with the minimum ring value. The identification is repeated for all successive time-windows with a 1-year step. The cumulated results point the years with minimum ring values over all trees in a given forest stand or plot. We tested different windows width, from 7 to 15 years and obtained the best agreements with the reference outbreaks time-series of Les Combes with a 7 and a 9-year window, according to the tested ring variables. Then we applied the same approach to the tree-ring time-series collected on the trees located at the other sites, at elevations ranging from approximately 1200-2350 m.a.s.l. (Table 1). We tested 2 year indexes: Index 1 is the cumulated sum of the number of trees over a given time-window for which the value of a ring variable is minimum. The second one is the sum of the values of the ring variable for which the value of a ring variable is minimum. For Index 1 the results are expressed as an annual unitless index while for the Index 2 the results are expressed in the same unit as the ring variable. Finally we also tested a combined index expressed as Index 1 × Index 2.
The elevational shift was tested with an ANOVA model where x ij is the LBM observation variable or the tree-ring index, µ is the mean, y j is the year variable, z j is the annual climate variable, and ε ij is the residual. We used the R aov function to perform the ANOVA and the gvlma package to test and validate the linear model assumptions.

Weather Variation at the Experimental Site During the Study Period
There are significant linear time trends for two climatic variables, minimum and maximum temperature, but not for rainfall and the De Martonne drought index (Figure 2). These relationships are quite strong for minimum temperature (Adj. R-square = 0.44) and for maximum temperature (Adj. R-square = 0.55). They correspond to a 0.03 and 0.06 • C annual increase for minimum and maximum temperature, respectively. Minimum and maximum temperature significantly increased between 1964 and 2016, while there was no significant linear variation of rainfall nor of the De Martonne aridity index during the period (Figure 2). We used the thresholds shown on Figure 2 to identify the coldest and driest years of the period. The direct measurements of the caterpillar abundance show six successive outbreaks and 2 years with missing data (Figure 3).

Relationships Between
According to Figure 3, there are three major wave-shaped outbreaks centered in 1962, 1971, and 1979, followed by three smaller ones centered in 1996, 2006, and 2015 at Les Combes, 1825 m.a.s.l. The major outbreaks peaks are preceded and followed by one, two, and, more rarely, three or more years with less important LBM density values. The seven outbreak episodes are separated by 2-5 years intervals with extremely low LBM density, very close to 0. The 1986 and 1987 data is not available. According to this graph, the time-lag between two successive outbreaks is between 8 and 10 years.
Three ring-variables index established on the Les Combes trees show a very significant agreement with the LBM outbreaks observed at Les Combes: ring width (r = 0.56), latewood density (r = 0.74), and ring density standard deviation (r = 0.75) (Figure 4).
The Figures 5-7 show the application of the ring index validated at Les Combes at the four elevations of an elevational gradient (1350, 1700, 2000, and 2300 m.a.s.l.).
The reconstructed LBM outbreaks index from 1967 to 2007 at the four elevations shows that the agreement with the Les Combes outbreaks is better for the Latewood density index than for the Ring density standard deviation index and overall than for the Ring width index (r = 0.76 at the maximum for latewood density, while it is respectively 0.63 and 0.34 for Ring density standard deviation and Ring width, Figures 5-7).
The Les Combes plot is located on a south facing slope at 1826 m.a.s.l., intermediate between that of the gradient 1700 and 2000 plots ( Table 1). The tree-ring index of these two plots show generally the best agreement with the Les Combes LBM outbreaks, compared with the two other plots, located at 1355 and 2347 m.a.s.l. The four plots of the altitudinal gradient are located on a north-facing slope right in front of Les Combes (number 5, 6, 7, and 8 on Figure 1).
We used the same method for calculating LBM outbreaks index at two extreme elevations for larch in the Briançon region, 1200 and 2350 m.a.s.l. (Figure 8). In this case, the agreement with the Les Combes outbreaks is very poor (never significant). The outbreak intensity is low, compared to the intensity measured at Les Combes and estimated along the elevational gradient.

DISCUSSION
The global warming trend at Briançon is mostly linear, especially for maximum temperature (Figure 1, R-square = 0.55). The increase rate of annual temperature is nearly 3 • C for the last 50 years and is consistent with other estimations showing that the current warming is greater at upper elevation (Mountain Research Initiative Edw Working Group et al., 2015). Rainfall seems stable and while the de Martonne index does not vary linearly with time, the number of dry years is much higher after 1990 than before (two before, six after according to the selected threshold in Figure 1).
We observe in Figures 4-8 that the LBM outbreaks correspond to a massive reduction of ring width, latewood density and ring density standard deviation. The ring width reduction is systematically delayed by one year, while the  latewood density and the ring density standard deviation are synchronized with the LBM outbreaks. The 8-10 years frequency of Les Combes outbreaks is consistent with previous observations (Dormont et al., 2006). The outbreak intensity is high in 1962, 1971, and 1979 and much lower in 1996, 2006, and 2015. Such weakening of the LBM outbreaks since 1981 has been already noted by several authors (Baltensweiler, 1993;Battipaglia et al., 2014;Saulnier et al., 2017).
The quality of the reconstructed outbreaks is variable according to the ring variable used: among the six variables that we tested (ring width, earlywood width, latewood width, ring density, earlywood density, latewood density, and standard deviation of ring microdensity) we find that three provide robust and usable reconstructions of LBM outbreaks: ring width, latewood density, and standard deviation of ring microdensity. Latewood density and standard deviation of ring microdensity provide the best reconstructions, with little difference between them (r = 0.74 and 0.75 p < 0.001, respectively, Figure 4). The quality of the ring width reconstruction is slightly lower (r = 0.56, p < 0.001). Arbellay et al. (2018) found that latewood width and blue intensity provided more precise reconstructions that ring width, but they did not test the ring density variables. The Les Combes reconstruction demonstrate that the 1986-1987 missing data correspond to a LBM outbreak of high or quite high intensity.
If we compare the LBM reconstructions along the elevational gradient with the LBM outbreaks at Les Combes, we find again a good agreement for latewood density (from r = 0.60 to 0.76 according to the elevation, p < 0.001), a lesser one but still quite good for standard deviation (from r = 0.52 to 0.63 according to the elevation, p < 0.001) and a poor one for ring width (from non-significant at 1350 and 2300 m.a.s.l. to 0.34 and 0.33, p < 0.05, respectively at 1700 and 2000 m.a.s.l.): latewood density seems to be a promising alternative to ring width for LBM reconstructions (Figures 5-7). To our knowledge, it is the first time that latewood density and ring density standard deviation were shown to be interesting alternatives to ring width for LBM outbreaks reconstruction. Latewood width was sometimes found to show relevant patterns (Weber, 1997;Baltensweiler et al., 2008;Arbellay et al., 2018), but in our case this ring variable was not better associated to the LBM outbreaks than ring width.
There is a systematic 1-year shift between the LBM outbreaks and the ring width reconstruction, while there is a precise synchronization between the LBM outbreaks and the latewood density reconstruction: the LBM defoliation decreases latewood density of the current annual ring, while it decreases ring width of the following year ring. The LBM defoliation decreases the photosynthetic activity at a time of the growing season where assimilates are mainly used for latewood formation. The decrease of latewood density can be the consequence of the formation of wider-lumen latewood cells, or of a lower proportion of narrow-lumen latewood cells (Bjorklund et al., 2017), associated for example with an earlier cessation of latewood formation at the end of the growing season. Such modification of latewood anatomy may have consequences on ring hydraulic properties and on larch resistance to drought. Ring-width indices of resistance to drought have been developed in larch (George et al., 2017). The key-role of latewood characteristics on hydraulic properties has been shown in other conifer species like Douglasfir (Dalla-Salda et al., 2014). But to our knowledge no literature is available concerning the consequence of a wood density decrease (either in early-or latewood) on larch hydraulic properties. At the same time, the diminution of the photosynthetic activity and/or the earlier growing season end probably decreases reserve stock. Reduced reserve may slow down earlywood formation when growth resumes during the next spring and explains the overall narrower next-year ring (Weber, 1997).
We observe that the cold 1984 year is associated to apparent LBM outbreaks in four reconstructed time-series, confirming that climatic anomalies can introduce confusion into the reconstructed time-series (Baltensweiler et al., 2008). We grouped the three low-elevation plots and the two high-elevation plots with the longer (up to 2017) reconstructed LBM time-series into two elevation categories (low and high). There is a good association between the Les Combes LBM outbreaks and the LBM reconstruction along the elevational gradient. On the other hand, the LBM reconstructions in the five extreme-elevation plots (grouped into low and high elevation plots) shows very little, mostly non-significant, agreement with the Les Combes LBM outbreaks. These plots are located at the low and high elevation margins of the local distribution of L. decidua: most of the time LBM outbreaks are observed at elevations somewhere between 1600 and 2200 m.a.s.l. (Weber, 1997;Rolland et al., 2001;Baltensweiler et al., 2008;Saulnier et al., 2017). The five-extreme elevation plots are located at 1200 (number 2, 3, and 4, Table 1 and Figure 1) and close to 2400 m.a.s.l. (number 9 and 10, Table 1 and Figure 1), in climatic conditions that could be too extreme for Z. griseana: too warm at the bottom and too cold at the top.
Our results confirm the recent perturbations of the LBM cycle observed in different regions of the Alps, with a weakening of the outbreaks since the 1980s (Baltensweiler, 1993;Battipaglia et al., 2014;Saulnier et al., 2017). The last outbreak we observed occurred mainly in 2015 and 2016 at 2000 and 2200 m.a.s.l. elevation, that is 10 years after the previous one. The direct observations at Les Combes show a significant weakening of the LBM outbreaks between 1960-1990 and 1990-2015. But this trend is not observed when the LBM outbreaks time series are based on the tree-ring index. The number of LBM outbreaks effectively captured by the tree-ring index is low (five) thus probably too low to highlight any time trend. Furthermore, some relatively cold (1996) and very dry (2004, 2005, and 2007) years may have affected the outbreak reconstruction and decreased their statistical power.
Larch budmoth outbreaks are often described as traveling waves originating from regions with a high concentration of favorable conditions (Bjørnstad et al., 2002;Johnson et al., 2006;Price et al., 2006;Saulnier et al., 2017). A very low speed of traveling wave could explain some of the spatial and temporal discrepancies observed between our plots, which are located only a few km apart the ones from the others. Other local factors like terrain, aspect and continuity of the forest cover may play a more important role. Several authors suggested that the ongoing climate change was responsible of these perturbations. Iyengar et al. (2016) proposed a model showing that the periodicity of the LBM outbreaks could increase to 40 years. So far our results show only a small increase of the periodicity, from 8-9 years to 10-11 years at the maximum. The available Météo-France climate data demonstrate that there is a conspicuous warming from 1964 to 2017, with no significant change in the precipitations. This warming already affected the larch growing conditions along its elevational distribution, degrading its growth conditions through an increase of the drought stress during the growing season at the bottom of its elevational distribution (Saderi et al., 2019). Several authors hypothesized that the global warming will improve the growing conditions at the top of the distribution of the mountain forest species like larch (Primicia et al., 2015;Sidor et al., 2015). The spatial and temporal changes in the tendencies of the LBM outbreaks could reflect the changes in the hosttree and in the insect physiology and the modifications of their biological interactions.

CONCLUSION
The study of the annual-ring time series of direct observations at Les Combes and along the elevational gradient showed a recent but limited elevational shift of LBM outbreaks from medium toward higher elevations, over 2000 m. We suggest that this elevational shift is a consequence of the variable effect of the global warming at the different elevations. Longer and warmer summer improves larch life conditions at the top of the gradient while increased summer drought may have a harmful effect at the bottom. Winter warming is also expected to affect differently the timing of LBM egg hatch and that of larch bud flush, larvae being at present susceptible to emerge whereas no needles are available as food at the former optimal elevation. A better synchronization between larch and LBM may exist at higher elevations. Whether the larch ecosystem will be able to adapt at all points of its elevational distribution is still an open question.

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

AUTHOR CONTRIBUTIONS
AR collected and analyzed the insect direct observations. PR and LP collected and analyzed the tree-ring data. FH provided and analyzed the climatic data. PR analyzed the data and wrote the first version of the manuscript. AR and LP reviewed the manuscript. All authors contributed to the article and approved the submitted version.