Relative impacts of gypsymoth outbreaks and insecticide treatments on forest resources and ecosystems: An experimental approach

1 Terrestrial Ecology Research Group, School of Life SciencesWeihenstephan, Technical University ofMunich, Freising, Germany 2 Department of Forest Protection, Bavarian State Institute of Forestry, Freising, Germany 3 Chair of Forest Growth and Yield Science, School of Life SciencesWeihenstephan, Technical University ofMunich, Freising, Germany 4 Field Station Fabrikschleichach, Department of Animal Ecology and Tropical Biology, Biocenter, University ofWürzburg, Rauhenebrach, Germany 5 EcosystemDynamics and Forest Management Research Group, School of Life SciencesWeihenstephan, Technical University ofMunich, Freising, Germany 6 Bavarian Forest National Park, Grafenau, Germany 7 Department of Biodiversity, Nature Conservation andHunting, Bavarian State Institute of Forestry, Freising, Germany


INTRODUCTION
The gypsy moth Lymantria dispar L. (Lepidoptera: Erebidae), a species of tussock moth indigenous to Europe and Asia, is widely acknowledged as one of the most critical defoliators in the Holarctic region. It exhibits cyclical population dynamics with outbreaks that are often periodic and spatially synchronous over large distances (Haynes, Liebhold, & Johnson, 2009) and believed to be driven by one or more of its trophic interactions (Johnson, Liebhold, Bjørnstad, & McManus, 2005). Gypsy moth larvae are voracious folivores that can feed for up to 10 weeks per year on more than 300 species of broadleaf and coniferous trees (Elkinton & Liebhold, 1990). This combination of outbreak dynamics, broad polyphagy and long larval duration granted the species an unmatched ability to rapidly defoliate large forested areas. Multiple large outbreaks have been reported during the last decades throughout its native range, including most of Europe (Alalouni, Schädler, & Brandl, 2013;Lentini et al., 2020;Villemant, 2010;Wulf & Graser, 1996), North Africa (Villemant, 2010), Central Asia (Orozumbekov, Liebhold, Ponomarev, & Tobin, 2009) and Japan (Liebhold, Higashiura, & Unno, 1998). However, the gypsy moth owes its infamous reputation mostly to the disastrous impacts of its eruptions in North American forests. Since its introduction in Massachusetts in 1869, it has colonized most of the Eastern United States and Southeastern Canada and keeps expanding its range (Liebhold, Halverson, & Elmes, 1992;Régnière, Nealis, & Porter, 2009;Tobin et al., 2004). Recent cost estimates including government expenditures and timber losses, appraised the overall financial toll of the gypsy moth invasion in North America to a staggering US$3.2 billion per year (Bradshaw et al., 2016). In Europe, the gypsy moth problem is generally considered less severe (McManus and Csóka, 2007) although a comprehensive synthesis of its economic impacts is lacking. Estimates of management costs can merely be extrapolated from published reports on the spatial extent of outbreak areas treated with insecticides (Alalouni et al., 2013), while data related to market costs remain scarce and scattered. Likewise, our current knowledge on the ecological impacts of gypsy moth outbreaks and their associated management practices is mostly drawn from studies conducted in the invasive range of the species. Both defoliation and control measures with synthetic and biological insecticides have adverse effects on forest ecosystems (Herms, 2003). Experimental approaches integrating economic and ecological aspects of both outbreaks and control practices would grant a powerful tool to improve the economic and environmental sustainability of gypsy moth management but are rare. In the present study, we first reviewed the different pathways through which a gypsy moth outbreak and insecticide use affect trees and animal communities in forest ecosystems. We then used the synthesized knowledge to produce biologically sound hypotheses ( Figure 1) as a baseline to develop a field experiment that aims to address the relative importance of these processes in the indigenous range of the gypsy moth.

Life cycle of the gypsy moth
The gypsy moth has a univoltine life cycle throughout its range. Adult females lay egg masses indiscriminately on bare surfaces in the vicinity of host trees in midsummer, each cluster typically containing 100 to over 1200 eggs (Andresen et al., 2001;Doane & McManus, 1981). The larvae are already fully developed inside the egg 1 month after oviposition and enter diapause in the winter months, eventually hatching during the spring bud-burst. The young caterpillars hang from branch tips with long silken threads, allowing them to spread to other forests via airborne dispersal. The gypsy moth is described as a spring-early summer defoliator. In Central Europe, feeding generally starts late April and continues for up to 10 weeks, until pupation in late Juneearly July. Pupation takes place in various sheltered places, for example in the ground or between bark ridges, and lasts for approximately 2 weeks. Adults are in wings from July to early August and generally live for 6-10 days. Females produce pheromones to attract males and lay egg masses shortly after mating close to their site of emergence (Doane & McManus, 1981;Wellenstein & Schwenke, 1981).

Ecological impacts of gypsy moth outbreaks
Gypsy moth outbreaks primarily affect tree health through defoliation, which can in turn significantly reduce timber production. The primary effects of gypsy moth defoliation on host trees include reduced growth (Muzika & Liebhold, 1999;Naidoo & Lechowicz, 2001), reduced root biomass production (Thomas, Blank, & Hartmann, 2002) and reduced fruit production (Gottschalk, 1990). Defoliation can also induce substantial levels of tree mortality (Campbell & Sloan, 1977;Kegg, 1973;Lobinger, 1999), which vary in intensity depending on the frequency, magnitude and duration of defoliation, the vitality of host trees, as well as the additive action of biotic or abiotic stresses F I G U R E 1 Theoretical representation of an experimental block with the main processes connecting the focal elements of the study system. (1) tebufenozide treatment; (2) gypsy moth population; (3) foliage canopy; (4) tree growth; (5) regeneration / understory vegetation; (6) non-target Lepidoptera; (7) parasitoids community; (8) predator community. Each quadrant represents one treatment combination: high defoliation risk, unsprayed; low defoliation risk, sprayed with tebufenozide; low defoliation risk, unsprayed; and high defoliation risk, sprayed with tebufenozide. Treatment is randomly allocated to one plot of each defoliation risk class within each block. The arrows represent the bivariate relationship, coloured according to the hypothesized direction of the effect. Blue colour denotes a positive relationship; purple colour denotes a negative relationship; black colour denotes a relationship which direction is conditional, for example guild-or species-specific. The arrow thickness represents the expected strength of the effect (Davidson, Gottschalk, & Johnson, 1999;Lobinger, 1999). Defoliation and tree death can disrupt the community structure of regeneration and herbaceous vegetation through competitive effects among subcanopy species released by increased light availability in canopy gaps (Fajvan & Wood, 1996).
Damage caused to vegetative and reproductive parts of trees can, in turn, induce bottom-up effects on animal communities. Increased proportion of snags and deadwood may benefit species associated with deadwood and tree cavities (Koenig, Walters, & Liebhold, 2011). Strong reduction of green leaf biomass has been shown to negatively affect herbivores sharing gypsy moth host plants, especially specialist species (Work & McCullough, 2000), while damage-induced changes in leaf chemistry can impede growth and survival of competing herbivores (Nykanen & Koricheva, 2004;Redman & Scriber, 2000). Furthermore, secondary infestations by leaf pathogens can inhibit later-season herbivory (Csóka, Pödör, Nagy, & Hirka, 2015;Tack, Gripenberg, & Roslin, 2012) and mast failure may negatively affect species that rely on acorns as a food source (Clotfelter et al., 2007). Additionally, outbreaks can drive an increased activity of predators and parasitoids, intensifying top-down pressures on other arthropod species (Faeth, 1986;Redman & Scriber, 2000). However, gypsy moth impacts can be difficult to quan-tify at the community level, as the relative importance of bottom-up and top-down effects may vary dramatically among species (Timms & Smith, 2011). In fact, analyses of lepidopteran assemblages in gypsy moth-infested forests reported only minor alterations of community structure compared to stands with non-outbreak situations, despite strong negative responses by single species (Sample, Butler, Zivkovich, Whitmore, & Reardon, 1996;Timms & Smith, 2011;Work & McCullough, 2000).
Gypsy moth outbreaks generate resource pulses as erupting caterpillar densities offer a superabundance of live prey to insectivorous predators, influencing their behaviour and population dynamics. The activity of species well-adapted to feed on the gypsy moth, such as North American cuckoos Coccyzus sp., and the forest caterpillar hunter Calosoma sycophanta (Coleoptera: Carabidae), was shown to be strongly correlated with the density of caterpillars (Gale, DeCecco, Marshall, McClain, & Cooper, 2001;Weseloh, 1985). Moreover, nitrogen-rich faeces, corpses and prematurely abscised foliage cause a spike in nitrogen input in the forest soil (Lovett, Canham, Arthur, Weathers, & Fitzhugh, 2006), which can, in turn, lead to increased nitrogen leaching when associated with substantial tree mortality (Lovett et al., 2002). The stimulated development of subdominant vegetation in canopy gaps increases habitat availability to understory-nesting birds (Cooper, Dodge, Whitmore, & Smith, 1993;Gale et al., 2001), but increased canopy openness also promotes nest predation (Thurber, McClain, & Whitmore, 1994).

Ecological impacts of insecticide treatments used against the gypsy moth
In the United States, insecticides have been used to combat gypsy moth populations since the inception of the invasion in the late 19th century (Liebhold & McManus, 1999). Even though hopes of eradication were definitively abandoned during the 1960s, insecticides are still used as the primary control measure to contain gypsy moth expansion and protect trees from defoliation in the current distribution area (Natural Resources Canada, 2020;Tobin et al., 2004;USDA Forest Service, 2020). Insecticides have also been frequently employed throughout the species' native range, although their use differs among countries (Alalouni et al., 2013;Lentini et al., 2020;Villemant, 2010;Wulf & Graser, 1996).
The moulting disruptor diflubenzuron was widely used throughout the 1980s in the United States (Liebhold & McManus, 1999) and was still applied in some European forests until recently (Schönfeld, 2009).
Although diflubenzuron is not toxic to vertebrates, concerns were raised with regard to its adverse effects on invertebrates (Durkin, 2004a) and the potential effects of its metabolite 4-chloroalinine on human health (European Food Safety Authority, 2012). The microbial insecticide Bacillus thuringiensis var. kurstaki (BTK) and the moulting hormone agonist tebufenozide were later approved as biorational alternatives due to their Lepidopteran-specific activity (Liebhold & McManus, 1999) and are nowadays the principal insecticides used in gypsy moth management. While BTK is applied over most of the treated area in North America (USDA Forest Service, 2020), tebufenozide is the preferred option some European countries such as Germany, due to its greater reliability in effectively suppressing gypsy moth populations (Lemme, Lobinger, & Müller-Kroehling, 2019).
Toxic side-effects of BTK and tebufenozide appear to be restricted to non-target Lepidoptera (Durkin, 2004b and references therein; Durkin & Klotzbach, 2004 and references therein), though environmental risk assessment studies involving tebufenozide in forests are comparatively scarce (but see Butler, Kondo, &Blue, 1997 andLeroy et al., 2019).
Forest spraying is usually carried out in mid-spring soon after gypsy moth larvae start feeding on leaves. At this time, insect populations in tree canopies are large and strongly dominated by Lepidoptera larvae (Martinat, Coffman, Dodge, Cooper, & Whitmore, 1988;Southwood, Wint, Kennedy, & Greenwood, 2013). Therefore, concerns have been raised about the impacts of depressed caterpillar abundances on insectivorous predators and parasitoids, mainly on forest birds that strongly rely on caterpillars to fulfil their energy requirements during the breeding season (Cooper, Dodge, Martinat, Donahoe, & Whitmore, 1990).
While most studies showed no substantial effect on bird communities, alterations in foraging behaviour, reproductive success and growth were observed in individual species in treated forests (Awkerman et al., 2011;Durkin, 2004a, 2004b andreferences therein;Holmes, 1998;Lih, Stephen, Smith, Nagy, & Wallis, 1994). Diet shifts and relocation of foraging territories were also documented in rodents (Bellocq, Bendell, & Cadogan, 1992;Seidel & Whitmore, 1995). Indirect effects on predatory arthropods are comparatively poorly known. Klenner (1996) observed reduced activity density of Carabid beetles in stands treated with diflubenzuron but was unable to determine whether this pattern was primarily driven by direct toxicity or reduced prey availability. While host scarcity is expected to negatively affect parasitoid populations, parasitism of sub-lethally intoxicated hosts appears to differ among species and could be influenced by insecticide-specific physiological effects on the gypsy moth hosts (Erb, Bourchier, van Frankenhuyzen, & Smith, 2001;McCravy, Dalusky, & Berisford, 2001;Weseloh et al., 1983).

Relative importance of defoliation and insecticides on forest ecosystems: Consequences for an experimental approach
Based on a comprehensive review of the respective effects of defoliation and BTK on non-target Lepidoptera, Scriber (2004) suggested that the decision of non-intervention could be as damaging for Lepidoptera communities as insecticides. This hypothesis has yet to be tested in an integrative experimental approach to address the relative importance of the numerous processes at play (Figure 1). Such an approach should consider the fact that outbreak intensity is spatially variable. Due to differences in abiotic factors and trophic interactions, the development of gypsy moth populations may widely differ among forest stands, potentially leading to a considerable variation in the magnitude of their economic and ecological impacts (e.g. Sample et al., 1996). Because insecticides must be applied shortly after larvae start feeding in order to be effective, decisions for spraying are mostly based on counts of egg masses on tree trunks conducted during the previous winter. These surveys cannot systematically predict severe defoliation, as most of the damage is typically caused by late instar larvae, up to 10 weeks after the initiation of feeding (Doane & McManus, 1981). Therefore, an experimental approach tailored to address the relative impacts of gypsy moth outbreaks and insecticide treatment should fulfil the following conditions: (1) selection of forest stands with broadly differing egg mass densities, to allow contrasting outbreak from non-outbreak conditions; (2) selection of sufficiently large experimental areas within stands, to be able to study processes over a biologically realistic scale; (3) exclusion of sites with a recent spraying history so that results are not confounded by potential carry-over effects of past treatments; (4) sufficient replication to account for unforeseeable variation in gypsy moth population dynamics among sites; (5) use of a randomized block design with closely located and structurally comparable stands within experimental blocks, to control for spatial heterogeneity among sites; (6) selection of blocks with different stand and climatic conditions, to allow the generalization of the results.
Here we describe a large-scale field experiment in Germany that fulfils these conditions. We established a full factorial block design including 12 blocks each composed of two plots with predicted outbreak densities of gypsy moth, and two plots with low predicted densities. An insecticide treatment with tebufenozide was randomly attributed to one plot in each density class. The experiment has three main objectives: (1) fostering our knowledge of gypsy moth population dynamics and impacts in its native range; (2) providing comprehensive information on the direct and indirect impacts of tebufenozide on forest ecosystems; (3) evaluating the relative impact of gypsy moth outbreaks and tebufenozide-based treatments to help to improve the sustainable management of outbreaks. In the present article, we first introduce the site selection process, the experimental design and the variables that will be measured in the experiment. We then discuss the viability of our design using gypsy moth population and defoliation data from the first experimental year.

Experimental area and gypsy moth population surveys
The experiment was set up in the region of Franconia, in north-western

F I G U R E 2
Map of the study design. Ellipses with capital letter labels represent the experimental blocks. Plots at high defoliation risk are displayed in red, plots at low defoliation risk in blue. Lighter colours indicate aerial treatment with tebufenozide; darker colours indicate unsprayed controls. Note: Three plots (one in block F, two in block D) are not displayed on the map as local landowners only agreed to take part in the experiment on the condition that the exact location of their stands is not published Local district foresters carried out surveys of gypsy moth egg masses in four administrative regions -Upper Franconia, Middle Franconia, Lower Franconia and Swabia -during fall 2018, following a standardized protocol. The number of gypsy moth egg masses were counted on the lowest 2 m of tree trunks along a transect comprising, in most cases, 10 trees of the dominant social class. The abundance of egg masses on the underside of lower canopy branches of each tree, stand vitality, stand age and history of previous outbreaks were also reported. These data were used to calculate a 'defoliation risk index' (DRI) to identify areas at high or low risk of defoliation in the summer 2019 (Supplementary Information, file S1). In total, 26823 single trees were surveyed along 2802 transects.  (Table S4-1 and Figure S5-1 in the Supporting Information). Application proceeded in dry weather and low wind conditions (i.e. wind speed below 2.5 m s −1 ) and block-wise when applicable. block T could not be included due to time constraints). In addition to the central tree, its six nearest neighbours with a diameter at breast height higher than 7 cm were included to account for competitive influences (Prodan, 1968). Starting from the centre of each plot, 20 of these 'six-tree samples' were taken along transects in the four cardinal points (Figure 3(a)). Sample circles were positioned at 25, 50, 75, 100 and 125 m from the centre of the infested area (origin of the coordinate system in Figure 3(b)). These sample circles served to record the stand characteristics (e.g. basal area, standing stock). Within each sample circle, the tree with medium diameter was selected as the sample tree for the detailed sampling. Species, diameter, position and crown transparency of the neighbour trees were scored.

Data collection
In order to investigate long-term effects on the growth pattern of the oaks due to defoliation, the central oaks were additionally equipped with a permanent girth tape. During the 2019 growing season, the tapes were read and checked five times (April, June, July, September and November). The assessment was repeated in 2020 and should be extended for up to three additional years. In the future, the reading and checking will be carried out annually in autumn after completion of annual ring formation and in spring immediately before the start of the growth, in order to eliminate artefacts caused by winter swelling and shrinkage, and defects caused by manipulation or overstretching of the tension springs.
In January 2020, four fences ( caterpillars, including gypsy moth, to identify each specimen and its associated parasitoids to species. The by-catch will be sorted to order for further analyses. The second and third sampling rounds, that is the acute toxicity phase and the peak defoliation sampling were repeated in 2020 in the same fogging areas that were swapped between both time points. As the post-treatment recovery of oak-dwelling Lepidoptera is expected to take more than 2 years, canopy arthropod surveys should be extended for at least one additional year. was conducted at the same intensity in 2020. As no outbreak population remained in any plot after 2020, lighter surveys (e.g. on a smaller number of trees per plot) will be performed during the next years as routine monitoring of population density.

Statistical procedures
Statistical analyses will be performed in R 4.0.2 and upcoming versions (R Core Team, 2020).

Missing data
In cases of non-random missing data (e.g. an entire block could not be sampled) in independent variables, missing data will be dropped, and the analysis only performed on complete data. Covariates with non-random missing data could be dropped depending on their importance in the analysis. For data missing at random in essential variables, individual data points will be dropped when the proportion of missing values is below 5%. For randomly missing data ranging from 5 to 30% of the total, multiple imputations (n = 500) will be performed with the R package mice to impute NAs with realistic values computed based on information from relevant variables in the dataset (Buuren & Groothuis-Oudshoorn, 2011). Statistical models will then be applied to each of the imputed datasets and the results pooled by Rubin's rules (Rubin, 1987).

Statistical models
In order to analyse the effects of gypsy moth outbreaks and tebufenozide application on the various focal variables, we will use the following core linear mixed model based on the characteristics of the experimental design: For variables with non-normal error distribution, generalized linear mixed models will be used with the family distribution that best fits the response variable. Zero-inflation models will be used to investigate the post-spray response of groups for which strong insecticide effects are expected, that is gypsy moth and non-target Lepidoptera. Adjustments to the core model, such as inclusions of relevant covariates or nested random effects (e.g. 1|block/plot) will be made on a model-to-model basis, considering the existence of specific hypotheses justifying the inclusion of additional variables and the hierarchical level at which the effects are investigated (e.g. plot or individual tree). Because our experiment has a full factorial design, (generalized) linear models may be performed instead of mixed effect models in the absence of F I G U R E 4 Structural equation model to test the relative impact of direct and indirect effects of gypsy moth outbreaks and tebufenozide on the response variables. The 'treatment' predictor depicts the combination of defoliation risk (DR; high/low) and tebufenozide application (TBF; control/sprayed) characterizing the experimental plots. Two different outcomes (A and B) may be integrated into the model to study cascading effects (e.g. trophic interactions). Each arrow represents a (generalized) linear or generalized additive (mixed) model with a random effect structure relevant to the scale at which the analysis is conducted, when applicable. Additional covariates are not shown but may be included when relevant nested structure (i.e. plot-level analysis) and strong block effect. Lastly, generalized additive (mixed) models may be used to fit nonlinear relationships between the dependent and independent variables.
To separate direct and indirect effects of tebufenozide on the various non-target groups under study, we will perform structural equation models with gypsy moth density or foliation rate as mediators of the relationship between treatment and the independent variable(s).
The choice of the appropriate mediator(s) will be motivated by the pathways through which indirect effects on the independent variable are expected to occur, namely whether they are related to change in resources (foliage or caterpillars) or habitat degradation (Figure 4).

Covariates
A covariate will be included in a model only if it complies with the following criteria: (1) expected biological relevance, that is there is a sound hypothesis motivating its consideration in the analysis; (2) independence, that is the effect of the focal variable is not confounded with that of treatment and other covariates. The presence of confounding effects will be tested for each covariate in regression models with treatment as a predictor. Failure to accept the null hypothesis (i.e. no significant correlation between the covariate and the treatment) will lead to the exclusion of the focal covariate. To test for multicollinearity of covariates, we will perform a principal component analysis (PCA) on all potential covariates using the subset of plots considered for the analysis. In cases when two or more covariates are highly correlated, we will favour the variable which best describes the expected relationship with the independent variable. The covariates potentially included in further analyses are listed in the overview of all measured variables in  Hartig, 2020) and mgcv (additive models; Wood, 2017) Influence measures. To assess the influence of outliers on the results of the regression models, we will compute Cook's distance (D) for each observation and examine the observations with values of D that are substantially outstanding from the rest. These observations will be dropped from the data, and the model refitted. Outliers will be considered influential if they substantially bias the estimates, in which case they will be dropped from the analysis.
Inference. To test our hypotheses, we will use t-, F-, Wald or likelihood ratio tests depending on the family error distribution and the structure of the focal model. For (generalized) linear mixed models, we will use the options currently available, following recommendations from Bolker (2020)

Evaluation of the experimental design
In order to assess the suitability of the study design for addressing our research questions, we measured gypsy moth population density and the intensity of defoliation in the study sites during the treatment year (2019). We gathered data on gypsy moth population density at various stages of its development: egg masses, early-instar larvae (canopy; pre-spray and post-spray), late-instar larvae (burlap bands; live and dead), pupae (burlap bands) and imagines (see section 2.3 'Data collection' for a description of the sampling methods). The data were analysed following the statistical procedures described in the previous section. Characteristics of the individual models are shown in Table S7-1 in the Supporting Information.

DISCUSSION
Gypsy moth outbreaks and insecticide treatments can strongly impact tree health and forest animal communities and are thus an important subject in forest ecology and management. Our understanding of these effects, however, is still somewhat limited due to a number of methodological shortcomings. Designing a suitable experiment to study these effects is a challenging undertaking that requires a careful selection of forest plots with sharp differences in gypsy moth densities while controlling for spatial heterogeneity and confounding factors such as past All response variables were fitted to generalized linear mixed models with negative binomial error distribution, except the foliation ratio that was fitted to a beta regression model (Table S7-1 in the Supporting Information). Fixed effects include defoliation risk, accounting for pre-spray gypsy moth density and stand vulnerability to defoliation (Supplementary Information, file S1), and insecticide treatments. The interaction between defoliation risk and insecticide treatment was added when treatment effects were hypothesized to vary with the pre-spray density of the gypsy moth. Significant effects are highlighted in bold. All analyses were conducted on data aggregated at the plot-level when applicable. Wald tests were used to test the significance of fixed effects. Further information on the models are given in the Supplementary Information (Table S7-1) spraying, within a sufficiently large region to allow the generalization of results. We successfully implemented an experimental approach that overcomes these hurdles to address the relative impacts of gypsy moth outbreaks and the insecticide-based management of these outbreaks. While a factorial design crossing high and low gypsy moth densities with insecticide treatment has already been implemented in a previous study (Sample et al., 1996), our design should considerably increase the robustness of this approach, notably with the inclusion of twice as many replicates and the use of blocking to account for spatial heterogeneity. Additionally, our high-resolution monitoring of gypsy moth population across the included stands offers a unique perspective on the different trajectories followed by different spatially synchronized gypsy moth populations. This data should help us investigating critical factors driving gypsy moth densities and closely follow the ecosystem response as a function of the outbreak magnitude.
Multiple obstacles face insecticide trials in forests: aerial pesticide applications are highly regulated practices that require considerable administrative effort to obtain clearance and involves complex organization and logistics. High heterogeneity among stands and unfavourable issues in negotiations with landowners often prevent the inclusion of stands within an experimental design, such that a considerable number of stands must be surveyed in order to reach reasonable objectives with regard to replication. These problems are usually overcome with pseudo-replication (Schönfeld et al., ), reduced numbers of true replicates (Cadogan & Scharbach, 2003) or small-scale approaches . In a research field where the compliance with the principles of replication, randomness and blocking are more often the exception than the rule the inclusion of 12 replicates in a large-scale aerial forestry trial constitutes a great effort.
One year of gypsy moth monitoring showed that our site selection met the conditions required for the success of the approach.
Tebufenozide applications worked as intended by successfully suppressing gypsy moth populations, which remained low until adult emergence in all sprayed plots ( Figure 5). In line with our objectives, the selected high-and low-defoliation risk classes sharply differed in terms of gypsy moth densities, with egg masses and early instar larvae in average 12-and 10-fold more abundant in high-than in lowdefoliation risk plots, respectively (Table 1 and Figure 5). Outbreak plots experienced the highest defoliation, with foliation rates in high-control 71%, 67% and 64% that of low-control, high-treatment and low-treatment plots, respectively ( Figure 5). However, defoliation was characterized by considerable variation among outbreak populations. In eight out of the 11 scanned high control stands, the gypsy moth population rose to very high densities, causing near-complete defoliation of the oaks but also showed clear signs of collapse, with high larval mortality and a drop in pupal abundance in seven plots ( Supplementary Information, file S8). On average, nearly two-thirds of the larvae collected under burlap bands in high-control plots were dead, while mortality was comparatively low in plots at low defoliation risk ( Figure 5). Such high mortality suggest an infestation of the caterpillar populations by pathogens such as the nuclear polyhedrosis virus, which often plays a crucial role in the collapse of erupting populations (Elkinton & Liebhold, 1990). These findings emphasize the importance of a large sample size to account for the highly variable and poorly predictable course of multiple outbreak populations within the same region.
Unlike the other life stages, the behaviour of the adult populations significantly deviated from our initial expectations based on defoliation risk and insecticide application. The number of adult gypsy moths caught in high-risk plots was twice as high as in low-risk plots, while tebufenozide only reduced adult catches by 25%, such that the difference between sprayed and control plots was not statistically significant within defoliation risk class (Table 1 and Figure 5). This pattern was driven by male moths, while females did not differ among treatments ( Figure S9-1 in the Supporting Information). This result suggests spillover between sites, likely facilitated by the pairing of adjacent plots at similar defoliation risk as control and treatment areas within eight out of the 12 blocks ( Figure 2). Such short distances between treatment and control plots may be considered as a lack of independence. However, this is consistent with the practice of aerial spraying in Bavaria, F I G U R E 5 Plot-level abundance of different life stages of the gypsy moth and foliation ratio in the treatment year as a function of predicted defoliation risk (H = high; L = low) and insecticide treatment (C = unsprayed control, T = tebufenozide). From left to right, top to bottom: gypsy moth life cycle with time points corresponding to the sampling of response data; egg masses; canopy larvae (pre-spray); canopy larvae (post-spray); burlap larvae (live); burlap larvae (dead); foliation ratio; pupae; imagines (males + females). Boxplots show the raw data aggregated to the plot-level with the type and number of sampling units per plot indicated in the y-axis title of each graph. Points and error bars correspond to estimated marginal means and 95% confidence intervals for each treatment group. Different letters indicate significant differences among groups (multivariate-t-adjusted comparisons of estimated marginal means, α = 0.05) where treated areas are usually small (e.g. 7.5 ha median plot size in the operational application conducted in parallel to our project in 2019; Insecticide applications in forests are controversially discussed in the political sphere, on the one hand because of the potential longerterm effects of insecticides on biodiversity and on the other hand because of the potential income loss of forest owners when outbreaks are nleft unmanaged. Insect outbreaks in forests can have long-lasting consequences for ecosystems (Carson, Cronin, & Long, 2008), such as a decrease in tree growth following defoliation, or higher susceptibility of defoliated trees to secondary biotic and abiotic stresses. Similarly, defoliation and insecticide application may have long-lasting effects on the structure of animal communities. To better understand the implications of these processes for forest management, it is essential to monitor the impacts of gypsy moth outbreaks and insecticide treatments over several years. Our experiment provides the opportunity to study both the effects of gypsy moth outbreaks and insecticide application on short and medium timescales. Data produced as part of this effort during the coming years should help decision-makers to develop well-informed management strategies considering both the economic and ecological impact of gypsy moth outbreaks.