Teresa Kao

Ferdinand Kohle

Kai Ma

Tangi Aubert

Alexander Andrievsky

Ulrich Wiesner

Fluorescent Silica Nanoparticles with Well Separated Intensity Distributions from Batch Reactions

Keywords: Fluorescent Probes, Silica Chemistry, Organic Fluorescent Dye, Total Internal Reflectance Fluorescence Microscopy, Single Particle Imaging, Multiplexed Imaging

Silica chemistry provides pathways to uniquely tunable nanoparticle platforms for biological imaging. It has been a long-standing problem to synthesize fluorescent silica nanoparticles (SNPs) in batch reactions with high and low fluorescence intensity levels for reliable use as an intensity barcode, which would greatly increase the number of molecular species that could be tagged intracellularly and simultaneously observed in conventional fluorescence microscopy. Here, employing an amino-acid catalyzed growth, highly fluorescent SNP probes were synthesized with sizes < 40 nm and well separated intensity distributions, as mapped by single-particle imaging techniques. A seeded growth approach was used to minimize the rate of secondary particle formation. Organic fluorescent dye affinity for the SNP during shell growth was tuned using specifics of the organosilane linker chemistry. This work highlights design considerations in the development of fluorescent probes with well separated intensity distributions synthesized in batch reactions for single-particle imaging and sensing applications, where heterogeneities across the nanoparticle ensemble are critical factors in probe performance.

number of fluorescent tags that can be reliably distinguished using commercially available microscopes. This problem motivates a need for novel fluorescent tags that increase the number of molecular species that can be imaged in parallel, while maintaining the small sizes and target labeling capabilities of traditional fluorophores 7 .

Optical barcodes for highly multiplexed imaging have been realized in micron-sized beads exhibiting multiple spectral colors at distinguishable intensity levels [8][START_REF] Zhang | Microspheres with fluorescent spherical zones[END_REF][START_REF] Banerjee | Production of dyed polymer microparticles[END_REF] . These bead-based technologies typically rely on "bath dye" approaches for controlling fluorescence intensity, where solvent-swollen polymer particles are incubated in a solution of organic dye molecules, allowing dyes to absorb into the particles. Subsequent solvent removal traps dyes inside the particles [START_REF] Zhang | Microspheres with fluorescent spherical zones[END_REF][START_REF] Banerjee | Production of dyed polymer microparticles[END_REF] .

Using micron-sized particles, intensity barcodes with up to ten distinguishable intensity levels can be generated simply by changing the concentration of dyes in the bath [8][START_REF] Zhang | Microspheres with fluorescent spherical zones[END_REF][START_REF] Banerjee | Production of dyed polymer microparticles[END_REF] . Despite wide success in micron-sized particles, the development of nano-sized intensity barcodes remains a challenge 7 .

Highly-loaded microbeads may contain thousands of dyes per particle, allowing use of widely separated adjacent intensity levels 8 . For example, a factor of two difference in dye concentration may be used to generate microbeads that differ in dye-loading by an average of one thousand dyes. However, similar strategies cannot be applied in nano-sized systems, where small particle volumes necessitate use of substantially lower numbers of dyes per particle to prevent fluorescence quenching 11 . Under this constraint, a factor of two difference in dye concentration results in average differences of only a few dyes per particle and as illustrated in Figure 1, local fluctuations in dye concentration alone may cause overlap between intensity distributions [START_REF] Banerjee | Production of dyed polymer microparticles[END_REF] .

While previous reports have shown that this problem can be solved using programmable materials, including DNA and RNA 7,12 , only limited success in nanoparticle-based systems has been achieved 13,14 . Distinguishable intensity levels are difficult to achieve in nanoparticles using infiltration-based batch synthesis approaches even at increasing dye concentrations, due to low encapsulation efficiency 15 . Previous reports have shown that highly dye-doped nanoparticles with minimal fluorescence quenching may be synthesized by controlled covalent attachment of dyes during particle growth 16 . However, using these methods, the batch synthesis of nanoparticles with non-overlapping intensity distributions is complicated by heterogeneities across the ensemble, often exacerbated by secondary particle formation during growth 17 . Here, we have addressed these problems using a highly tunable synthesis method for covalently encapsulating organic fluorescent dyes inside silica nanoparticles (SNPs), by controlling secondary particle formation and varying the organosilane linker chemistry to the dye, demonstrating that high and low dye-loaded SNP intensity distributions are well separated as mapped by single particle fluorescence imaging.

Silica chemistry provides a uniquely tunable platform for nanoparticle synthesis, where particle size [18][19][20] , nanoscale morphology 21 , and surface properties 22,23 can be precisely controlled.

Recent advances in fabrication have extended the range of accessible particle sizes to below 10 nm, using water as a solvent 18 , and established efficient methods for surface functionalization with a variety of targeting moieties 22,23 , making such materials ideal candidates as probes for biological imaging. We use a previously reported method for amino acid catalyzed growth of SNPs in aqueous media 19,20 to investigate the synthetic parameters that govern single-particle intensity distributions in fluorescent dye-loaded SNPs. Using this basic framework, we demonstrate two different approaches for synthesizing bright SNPs with tunable average fluorescence intensity levels and sizes below 40 nm. Controlled distributions of dyes within individual SNPs are achieved by incremental shell growth, allowing dyes to remain unquenched even at high loadings. The resulting SNPs are up to 100 times brighter than the constituent dyes. Optimized intensity distributions, as mapped by single-particle imaging techniques, are achieved by using a seeded growth approach to minimize the rate of secondary particle formation and by tuning the affinity of the organosilane linker for the SNP during shell growth to improve the incorporation efficiency of the often-expensive fluorescent dye molecules, by means of functional groups able to hydrogen bond to silica surface silanol groups.

Fluorescent SNPs were prepared either by homogenous or seeded growth in the presence of a silane-conjugated, fluorescent cyanine dye, Cy3 (Figure 2a) (see Supporting Information for detailed methods). We utilized l-arginine as a base catalyst for SNP formation in aqueous media at elevated temperatures of 60 C 19,20,24 . As previously studied in detail 18,25 , under such conditions (i.e. aqueous medium at elevated temperatures) silica nanoparticle formation proceeds via the rapid formation of small (< 2 nm) silica clusters which subsequently aggregate into larger particles. In this system, however, a thin layer of tetraethoxysilane (TEOS) added to the top surface of the aqueous reaction volume serves as the silica source, with the interfacial area determining the rate of hydrolyzed TEOS diffusing into the water phase. For SNPs synthesized using what we call in the following the "homogenous" growth approach (Figures 2a andc), Cy3 dye was covalently coupled to an aminopropyltriethoxysilane (APTES) and added to the water phase before TEOS was layered on top. APTES was chosen for its relatively slow hydrolysis kinetics, as compared to aminopropyltrimethoxysilane (APTMS) 26 , to minimize self-aggregation and quenching of Cy3 dyes during particle nucleation. Due to the slow reaction rate of the sterically hindered Cy3/APTES conjugate, it is not fully consumed during the core synthesis. Three subsequent TEOS additions in 8 hr increments were therefore used to allow for continued attachment of unreacted Cy3/APTES to SNPs, resulting in highly dye-doped particles. After growth was completed, SNPs were stabilized by covalent attachment of polyethylene glycol (PEG) chains 22 . SNPs were subsequently purified using gel permeation chromatography (GPC) to remove any unreacted chemicals 18,23 (Figures S1c andd).

As a result of covalent dye encapsulation in the rigid silica matrix 4,27,28 , normalized absorption and emission spectra (Figures 3a) show that per dye brightness levels increase ~7 to 8 times over free dye in solution. Only slightly decreased per dye fluorescence enhancements are observed as the concentration of Cy3-silane is increased, suggesting that even at high loadings, dyes remain essentially unquenched. Cy3-labeled SNPs were further characterized using fluorescence correlation spectroscopy (FCS) 27 , where fluctuations in fluorescence intensity generated by particle diffusion into and out of a detection volume of a confocal setup are autocorrelated and fit using a model describing single-component diffusion 4,29,30 . Importantly, a single FCS measurement provides information over hydrodynamic particle size, particle concentration, as well as particle brightness 27 , making it a powerful tool in fluorescent nanoparticle characterization. A 543 nm HeNe laser was used as the excitation source. Figure 3b shows that the average hydrodynamic sizes obtained from FCS fits of Cy3-labeled SNPs synthesized by homogenous growth increase as Cy3/APTES concentration is increased. Using 1.3 M and 24.0 M Cy3/APTES resulted in SNPs with diameters of 27.8 nm and 37.4 nm, respectively. The difference in particle size measured by FCS as a function of Cy3/APTES concentration is corroborated in the corresponding GPC elugrams (Figures S1c andd) of the particles. In this system, a typical GPC elugram exhibits three peaks as a function of elution volume, corresponding to particle aggregates, PEGylated particle product, and free reactants, respectively 18 . Shorter elution times were observed for SNPs prepared using 24.0 M Cy3/APTES (Figure S1d), indicating increased particle size relative to SNPs synthesized using 1.3 M Cy3/APTES (Figure S1c). Transmission electron microscopy (TEM) images and corresponding size analysis (Figure S1a, b) corroborate the APTES-dependent size difference observed in FCS and GPC. Since the Cy3/APTES conjugate is prepared using a 10-fold excess of aminosilane and the reaction is performed at basic pH, these observations are consistent with early aggregation of negatively charged silica clusters (vide supra) facilitated by increasing concentrations of positively charged APTES, resulting in increased average particle sizes 25,31 . Interestingly, the GPC elugram in Figure S1d for particles from the higher dye concentration (24.0 M Cy3/APTES) exhibits a more skewed particle size peak. This suggests that higher concentrations of positively-charged APTES leads to continuous secondary particle formation, as more and more TEOS from the top layer hydrolyzes and diffuses into the water phase, thereby increasing particle size dispersity 32 . Despite increased polydispersity, Figure 3c clearly demonstrates that SNPs with increasing average fluorescence intensity levels could be obtained by adjusting the concentration of Cy3-silane used during particle synthesis, resulting in SNPs that are 15 and 47 times brighter than free Cy3 dye. Unfortunately, due to extreme particle size dispersity at even higher APTES levels, the concentration of Cy3silane could not be increased beyond 24.0 M without complete loss of particle uniformity. Accurate identification of "high" and "low" intensity tags requires minimal overlap between single particle intensity distributions. We therefore evaluated batches of fluorescent dyeloaded SNPs for use as an intensity barcode, assessing single particle analysis based intensity distributions using total internal reflection fluorescence (TIRF) microscopy 14 . For imaging, biotinylated SNPs were immobilized on a streptavidin-coated glass substrate. Small numbers of biotin molecules were introduced to SNP surfaces during PEGylation using ,heterobifunctional PEG molecules (Figures 2a andc). Streptavidin-coated substrates were prepared using a previously established method 33 with slight modifications (see Supporting Information). Briefly, a thiol-coated substrate (blue in Figure 2b) was prepared by silanization of a glass slide (gray) using mercaptopropyltrimethoxysilane (MPTMS). N--maleimidobutyryloxysuccinimide ester (GMBS) was attached to the thiol-coated substrate using its reactive maleimide group and the NHS ester group (green) at the opposite end of the GMBS was coupled to a primary amine on the streptavidin molecule (purple). Substrates were rinsed thoroughly after each step to remove unreacted chemicals. Imaging was performed over a 51.1 x 51.1 μm 2 area containing ~100 SNPs. Single images were acquired for 100 msec. After image acquisition, a spot detection algorithm was used to identify individual SNPs. The background level and full width at half maximum (FWHM) of spots in the x-and y-directions were derived from Gaussian fits. Spot intensity was defined as the background-corrected sum over all pixels within the FWHM 34 (see

Supporting Information and Figures S2 and S3 for details).

We compared spot intensities from Cy3-labeled SNPs synthesized using the homogenous growth approach and 1.3 μM or 24.0 μM Cy3/APTES. The resulting histograms (Figure 4a), while showing an overall intensity increase for high Cy3-loaded particles, consistent with the increased average brightness/particle values measured using FCS (Figure 3c), substantially overlap. To benchmark the performance of these SNPs, we quantified the overlap between the two intensity distributions shown in Figure 4a. First, a threshold value was identified that minimizes overlap between high and low intensity distributions (dashed black line in Figure 4a). Overlap for high Cy3-loaded SNPs was then defined as the fraction of particles with spot intensity below that threshold value, while overlap for low Cy3-loaded SNPs was calculated as the fraction of particles with spot intensity above the threshold. Using these criteria, 15.3% and 13.4% overlap were calculated for high and low intensity distributions, respectively. This result indicates that on a single particle level, nearly 30% of particles cannot be unambiguously identified based on intensity alone. Clearly, such low identification accuracy would preclude use of these SNPs as an intensity barcode.

To overcome these limitations of the homogenous growth approach (i.e. high concentrations of Cy3/APTES causing extreme polydispersity in SNPs and limited separation from low particle brightness distributions), we introduced a second synthesis approach which we will refer to as the "seeded" growth method. Here, undyed SNPs with sizes < 20 nm were used as seeds for further particle growth (Figures 2a andd). Furthermore, Cy3 dye was covalently coupled to N-(2-aminoethyl)-3-aminopropyltrimethoxysilane (AEAPTMS) instead of APTES. The comparatively fast hydrolysis kinetics of AEAPTMS 26 obviates the prolonged continued growth used to achieve high dye-loadings in the homogenous growth approach. Thus, to maintain relatively small particle sizes (< 40 nm), only one additional TEOS growth step was used. SNPs were again PEGylated prior to purification using GPC. Seeded growth of SNPs was monitored using dynamic light scattering (DLS), demonstrating that incremental shell growth results in steadily increasing particle sizes (see Figure S4a in Supporting Information). Figure 3e shows that as compared to SNPs synthesized by homogenous growth (Figure 3b), the hydrodynamic sizes of SNPs synthesized by seeded growth do not vary significantly with increasing Cy3/AEAPTMS concentration. Using 6.0 M, 24.0 M, and 30.0 M Cy3/AEAPTMS resulted in SNPs with diameters of 34.4 nm, 34.2 nm, and 38.2 nm, respectively. Using the same concentration of Cy3silane, SNPs prepared by seeded growth also exhibit more symmetric particle peaks in GPC (Figure S4e) than SNPs synthesized by homogenous growth (Figure S1d), suggesting that fewer secondary particles are formed. TEM images corroborated that these SNPs are narrowly dispersed (see Figure S4c and d in Supporting Information).

Figure 3d shows that Cy3-labeled SNPs synthesized by seeded growth and using increasing amounts of Cy3/AEAPTMS exhibit per dye fluorescence enhancements of 6 to 7 times. As in SNPs synthesized via homogenous growth, increasing average brightness/particle was achieved by increasing the concentration of Cy3/AEAPTMS conjugate. Figure 3f shows that Cy3-labeled SNPs using 6.0 M, 24.0 M, and 30.0 M Cy3/AEAPTMS are 28, 91, and 106 times brighter, respectively, than free Cy3 dye. Comparison between Figures 3c andf shows that from 24.0 M Cy3-silane concentration upwards, the seeded growth approach is a more efficient method of generating bright fluorescent tags. Cy3-labeled SNPs prepared by seeded growth using 24.0 μM Cy3/AEAPTMS are nearly 2 times brighter than those synthesized by homogenous growth using Cy3/APTES at the same concentration, even though the per dye brightness is slightly decreased in these particles. This increased dye incorporation efficiency is attributed to reduced rates of secondary particle nucleation/formation achieved using the seeded growth approach and to the relatively fast reaction rates of Cy3/AEAPTMS as compared to Cy3/APTES. As this result could not be achieved with Cy3/APTMS (data not shown), we hypothesize that the secondary amines present in AEAPTMS may also contribute to the enhanced association with the growing SNPs by hydrogen bonding to silica surface silanol groups [START_REF] Kroger | Proc. Natl. Acad. Sci[END_REF] , further increasing the incorporation efficiency of Cy3/AEAPTMS relative to Cy3/APTES. The advantages of the seeded growth of fluorescent SNPs using Cy3/AEAPTMS (i.e. increased brightness/particle in FCS) particularly manifested in TIRF microscopy. Single particle intensity distributions for SNP batches synthesized using 6.0 M, 24.0 M, and 30.0 M Cy3/AEAPTMS show that an overall intensity increase is reproducibly observed as the concentration of Cy3/AEAPTMS is increased (see Figure S5 in the Supporting Information).

Importantly, single particle intensity distributions for batches generated using the seeded growth approach and 24.0 μM Cy3/AEAPTMS (Figure S5b) are also shifted to higher spot intensities as compared to particle batches derived from the homogeneous growth approach at the same concentration of Cy3/APTES (Figure 4a). Designating SNP batches synthesized by homogenous (1.3 M Cy3/APTES) and seeded growth, (30.0 M Cy3/AEAPTMS) as low and high intensity level tags, we approach near perfect identification accuracy. Figure 4d shows that single particle intensity distributions exhibit 0.2 to 0.3% overlap for low (purple) and high (green) Cy3-loaded SNPs, indicating that >99% of SNPs are unambiguously identifiable based on intensity alone. We applied our intensity barcoding scheme to fluorescence images of mixed high and low intensity probes (Figure 4c). Spots above and below the intensity threshold defined in Figure 4d are colored green and purple, respectively. Comparison to images of only high (Figure 4e) or only low intensity (Figure 4b) probes illustrates that highly accurate intensity-based identification can be achieved. While for these proof-of-principle demonstrations, batches from homogeneous and seeded growth were compared, we expect that after careful parameter optimization other combinations, for example of two batches synthesized using seeded growth, will show similar results.

We have described two batch approaches for the synthesis of fluorescent silica nanoparticles using either homogenous or seeded growth. Controlled fluorescence intensity levels can be achieved in these systems simply by varying the concentration of a Cy3-silane conjugate in the batch. We have shown that the precise spatial distributions of dyes within particles can be controlled by silica shell growth, minimizing self-quenching of dyes even in highly dye-doped particles. Comparison of the two approaches highlights that the nucleation of secondary particles during shell growth presents a major challenge in the synthesis of highly dye-doped SNPs with narrowly distributed single particle fluorescence intensity distributions. Using the seeded growth the rate of Cy3 incorporation by choice of silane linker (e.g. AEAPTMS versus APTES), minimizing the effects of secondary nucleation in this system. We find that single particle analysis based fluorescence intensity distributions of SNPs synthesized in these ways are well separated and can therefore be used as an intensity-based barcode. These SNPs are small (< 40 nm), highly biocompatible, and readily functionalizable, making them ideal fluorescent tags for imaging in living biological systems. We have used a single fluorescent dye (Cy3) as proof-of-concept, but expect that our results can be extended to combinations of dyes with multiple spectral colors, greatly increasing the number of SNP tags that can be imaged in parallel using conventional fluorescence microscopy. Additionally, since this strategy relies on the improved incorporation efficiency of a silane linker, we expect that similar results can be achieved using a range of organic fluorescent dyes and linker chemistries, providing an efficient route to encapsulation of these often-expensive molecules. single particle imaging, a thiol-coated substrate (blue) was prepared by incubating a glass slide (gray) in a solution of MPTMS in ethanol (EtOH). Heterobifunctional GMBS was attached to the thiol-coated substrate using its reactive maleimide group and the NHS ester group (green) at the opposite end was coupled to a primary amine on the streptavidin molecule (purple). Substrates were rinsed thoroughly after each step to remove unreacted chemicals. Biotinylated SNPs were attached to streptavidin-coated substrates just before imaging.

Figure 1 .

 1 Figure 1. Schematic illustration of a bath dye approach for synthesizing fluorescent polymer

Figure 2 .

 2 Figure 2. Schematic representations of the preparation of fluorescent SNPs and substrates for

Figure 3 .

 3 Figure 3. Characterization of fluorescent SNPs synthesized by (a-c) homogenous and (d-f) seeded

Figure 4 .

 4 Figure 4. Comparison of fluorescence intensity levels for low and high dye-loaded SNPs as

Supporting Information

The Supporting Information is available free of charge on the ACS Publications website. Materials, synthetic methods, characterization, and image analysis are described. Additional characterization of SNPs prepared by homogenous and seeded growth is shown. Datasets showing reproducibility of the synthesis and imaging method are also included.


Toigo Maude
email: maude.toigo@cefe.cnrs.fr

Nicolas Manuel

Jonard Mathieu

Luc Croisé

Nageleisen Louis-Michel

Jactel Hervé

Temporal trends in tree defoliation and response to multiple biotic and abiotic stresses

Keywords: drought, pest insect, foliar nutrients, die-off, interactions, temperate forest, crown defoliation, ICP Forests, Fagus sylvatica, Quercus robur, Quercus petraea

The vitality of forests depends on multiple biotic and abiotic stresses that can act in an interactive way. Tree defoliation levels are therefore subject to long-term changes and shorter-term hazards that need to be monitored in the context of ongoing climate change. In this study, we analysed the combined effects of drought, insect attack and tree nutritional status on the level of defoliation in three major hardwood species in Europe: European beech (Fagus sylvatica), sessile oak (Quercus petraea) and pedunculate oak (Quercus robur). The forest condition data set was collected in the long-term intensive monitoring network in France (RENECOFOR). The average level of defoliation observed was 20% but varied considerably between plots and trees. From 1997 to 2015, beech defoliation worsened by 10% while oak defoliation stabilized over time. For the three tree species, crown defoliation was generally higher in the event of insect attack (increased by 2.6% to 7.5% depending on the tree species) and under drought (increased by 5.9%), while the effect of the tree nutritional statue was less consistent (from -5.6% to +10%) these factors acting synergistically. This study highlights the fact that the different stresses a tree species undergoes act in a complex and interactive way, with species-specific responses. In a context of increasing abiotic and biotic stresses, their combined analysis appears to be a necessity at a time when forest owners need to adapt their management to cope with climate change, particularly through the choice of tree species to be favoured in the future.

1.INTRODUCTION

Coping with climate change is a major challenge for forest managers because the longevity of trees limits the capacity for rapid renewal of the species or varieties best adapted to the new abiotic constraints. To adapt existing forests or plan future plantations, they need an integrative vision of the effects of the multiple risks to which forest stands are exposed. Climate change has multiple important and interacting impacts on forests globally [START_REF] Mcdowell | Pervasive shifts in forest dynamics in a changing world[END_REF]. Forests are the terrestrial ecosystem that contributes the most to carbon sequestration worldwide [START_REF] Bonan | Forests and Climate Change: Forcings, Feedbacks, and the Climate Benefits of Forests[END_REF][START_REF] Janssens | Europe's Terrestrial Biosphere Absorbs 7 to 12% of European Anthropogenic CO2 Emissions[END_REF][START_REF] Pacala | Consistent Land-and Atmosphere-Based U.S. Carbon Sink Estimates[END_REF][START_REF] Piao | The carbon balance of terrestrial ecosystems in China[END_REF] but global warming could alter drastically functioning, reducing primary productivity and thus carbon sequestration [START_REF] Lindner | Climate change impacts, adaptive capacity, and vulnerability of European forest ecosystems[END_REF]. Forest growth is favoured by an increase in CO 2 concentration [START_REF] Norby | Forest response to elevated CO2 is conserved across a broad range of productivity[END_REF][START_REF] Walker | Decadal biomass increment in early secondary succession woody ecosystems is increased by CO2 enrichment[END_REF], moderate increase in temperature, and in nitrogen deposition [START_REF] Myneni | Increased plant growth in the northern high latitudes from 1981 to 1991[END_REF] as far as the other needs in water and nutrients are satisfied. Conversely, above a critical value, an increase in temperature and atmospheric compounds such as ozone, nitrogen and sulfur is expected to be detrimental to forest primary productivity through altered photosynthesis, reduced leaf area and leaf longevity [START_REF] Ainsworth | The Effects of Tropospheric Ozone on Net Primary Productivity and Implications for Climate Change[END_REF][START_REF] Juknys | Dynamic response of tree growth to changing environmental pollution[END_REF][START_REF] Magill | Ecosystem response to 15 years of chronic nitrogen additions at the Harvard Forest LTER, Massachusetts, USA[END_REF][START_REF] Wittig | Quantifying the impact of current and future tropospheric ozone on tree biomass, growth, physiology and biochemistry: a quantitative meta-analysis[END_REF].

However, these long-term trends can be affected by catastrophic events occurring at shorter time scale. Climate change is indeed characterized by an intensification of abiotic (fire, windstorms, drought) and biotic (insects and pathogens) hazards, leading to severe alteration of tree vitality [START_REF] Seidl | Forest disturbances under climate change[END_REF]. In particular, drought episodes have intensified in severity and frequency [START_REF] Diffenbaugh | Quantifying the influence of global warming on unprecedented extreme climate events[END_REF] and are predicted to increase further in the future [START_REF] Lehner | Projected drought risk in 1.5°C and 2°C warmer climates: Drought in 1.5°C and 2°C Warmer Climates[END_REF][START_REF] Ruosteenoja | Seasonal soil moisture and drought occurrence in Europe in CMIP5 projections for the 21st century[END_REF]. This highlights the need to analyse past longterm changes in forest conditions along with short-term variations due to biotic and abiotic perturbations in order to better predict the complex effects of climate change on forest ecosystem productivity and vitality [START_REF] Matyssek | Forests under climate change and air pollution: gaps in understanding and future directions for research[END_REF].

Tree growth and defoliation are important indicators of tree vitality [START_REF] Dobbertin | Tree growth as indicator of tree vitality and of tree reaction to environmental stress: a review[END_REF] since long interruption of tree growth or intense defoliation represent relevant indices of severe decline or imminent death of trees. Numerous studies have documented long-term changes in forest growth as well as biotic and abiotic determinants of tree growth through dendrochronological approaches [START_REF] Charney | Observed forest sensitivity to climate implies large changes in 21st century North American forest growth[END_REF][START_REF] Graumlich | Subalpine Tree Growth, Climate, and Increasing CO_2: An Assessment of Recent Growth Trends[END_REF][START_REF] Piovesan | Drought-driven growth reduction in old beech (Fagus sylvatica L.) forests of the central Apennines, Italy[END_REF], but similar studies on tree defoliation are less common as they require long-term monitoring of permanent plots. In Europe, no clear pattern emerged from observation of long-term trend in tree defoliation, the overall mean defoliation remains around 20% and changes in defoliation over time are not significant or show idiosyncratic responses of tree species (Michel, Seidling, Prescher, & editors, 2018) suggesting that local scale might be more relevant to detect anomalies (de la [START_REF] De La Cruz | Defoliation triggered by climate induced effects in Spanish ICP Forests monitoring plots[END_REF][START_REF] Eichhorn | Temporal development of crown condition of beech and oak as a response variable for integrated evaluations[END_REF][START_REF] Ferretti | Plot-scale modelling to detect size, extent, and correlates of changes in tree defoliation in French high forests[END_REF][START_REF] Popa | Influence of climate on tree health evaluated by defoliation in the ICP level I network (Romania)[END_REF][START_REF] Sousa-Silva | Tree diversity mitigates defoliation after a drought-induced tipping point[END_REF]. Moreover tree defoliation has long been considered a relevant marker of environmental changes due to its sensitivity to abiotic factors such as temperature, precipitation, wind, ozone concentration, soil pH and fertility [START_REF] Bussotti | Visible injury, crown condition, and growth responses of selected Italian forests in relation to ozone exposure[END_REF][START_REF] Bussotti | Leaf morphology and chemistry in Fagus sylvatica (beech) trees as affected by site factors and ozone: results from CONECOFOR permanent monitoring plots in Italy[END_REF][START_REF] Drobyshev | Crown condition dynamics of oak in southern Sweden 1988'999[END_REF][START_REF] Ferretti | Plot-scale modelling to detect size, extent, and correlates of changes in tree defoliation in French high forests[END_REF]. Biotic factors such as forest pests (e.g. defoliators), pathogens (e.g. leaf rusts), wild game but also tree competition (e.g. stand density and composition) are also known to cause tree defoliation [START_REF] Ferretti | Plot-scale modelling to detect size, extent, and correlates of changes in tree defoliation in French high forests[END_REF][START_REF] Pollastrini | Physiological significance of forest tree defoliation: Results from a survey in a mixed forest in Tuscany (central Italy)[END_REF]. Tree defoliation can thus be used as an integrative marker of multiple interactive stresses [START_REF] Vries | Effects of environmental stress on forest crown condition in Europe. Part I: Hypotheses and approach to the study[END_REF]. It is thus surprising that only few studies investigated the combined and interactive effects of biotic and abiotic stresses on tree defoliation.

Many studies have already reported tree mortalities caused by severe droughts around the world [START_REF] Allen | On underestimation of global vulnerability to tree mortality and forest die-off from hotter drought in the Anthropocene[END_REF][START_REF] Mcdowell | Drivers and mechanisms of tree mortality in moist tropical forests[END_REF][START_REF] Meir | Short-term effects of drought on tropical forest do not fully predict impacts of repeated or long-term drought: gas exchange versus growth[END_REF] and this is worsening as climate becomes hotter and drier [START_REF] Van Mantgem | Widespread Increase of Tree Mortality Rates in the Western United States[END_REF]. Water stressed trees are also known to be more prone to many pest insects and disease damages [START_REF] Bréda | Temperate forest trees and stands under severe drought: a review of ecophysiological responses, adaptation processes and long-term consequences[END_REF][START_REF] Jactel | Drought effects on damage by forest insects and pathogens: a meta-analysis[END_REF]. Not only drought but also higher temperatures can enhance insect population dynamics [START_REF] Jactel | Responses of forest insect pests to climate change: not so simple[END_REF], as shown by the largest bark beetle (Dendroctonus ponderosa) outbreak ever observed in North America [START_REF] Carroll | Effect of climate change on range expansion by the mountain pine beetle in British Columbia[END_REF] after a series of mild winters. Several hundred million m 3 of trees were killed, which transformed the forest ecosystem into a net source of atmospheric carbon [START_REF] Kurz | Mountain pine beetle and forest carbon feedback to climate change[END_REF]. Those examples illustrate the fact that multiple biotic and abiotic stresses can cascade or synergize to alter forest vitality.

Tree health is also closely linked to soil fertility, disequilibrium in soil nutrients and atmospheric depositions. Deficit or excess in a particular element may affect tree growth or crown condition. For example, high crown transparency of Picea abies has been related to deficiency in foliar potassium (K), nitrogen (N), phosphorus (P) and iron (Fe) nutrients [START_REF] Ewald | Ecological background of crown condition, growth and nutritional status of Picea abies (L.) Karst. in the Bavarian Alps[END_REF]. Fertilization trials showed that adding potassium could increase foliar K concentration and decrease defoliation in young balsam fir [START_REF] Ouimet | Effects of fertilization and liming on tree growth, vitality and nutrient status in boreal balsam fir stands[END_REF]. Moderate nitrogen atmospheric deposition can exert a fertilising effect, enhancing tree growth and decreasing crown defoliation [START_REF] Thimonier | Does exceeding the critical loads for nitrogen alter nitrate leaching, the nutrient status of trees and their crown condition at Swiss Long-term Forest Ecosystem Research (LWF) sites?[END_REF].

Conversely, high nitrogen deposition may cause nitrate leaching [START_REF] Thimonier | Does exceeding the critical loads for nitrogen alter nitrate leaching, the nutrient status of trees and their crown condition at Swiss Long-term Forest Ecosystem Research (LWF) sites?[END_REF] and soil acidification [START_REF] Flückiger | Nitrogen and Its Effect on Growth, Nutrient Status and Parasite Attacks in Beech and Norway Spruce[END_REF] affecting tree health. High level of foliar nitrogen concentration has been associated with increased tree sensitivity to frost, drought, pest and diseases (W De [START_REF] Vries | Impacts of acid deposition, ozone exposure and weather conditions on forest ecosystems in Europe: an overview[END_REF][START_REF] Flückiger | Nitrogen and Its Effect on Growth, Nutrient Status and Parasite Attacks in Beech and Norway Spruce[END_REF]. Moreover, the fertilising effect of nitrogen supply can be sustained only if other critical nutrients are not limiting. For example, nitrogen addition along with deficit in phosphorus results in unbalanced N:P ratio [START_REF] Flückiger | Nitrogen and Its Effect on Growth, Nutrient Status and Parasite Attacks in Beech and Norway Spruce[END_REF], which can reduce tree growth [START_REF] Braun | Does nitrogen deposition increase forest production? The role of phosphorus[END_REF] and increase crown defoliation [START_REF] Veresoglou | Exploring continental-scale stand health-N: P ratio relationships for European forests[END_REF]. In European forests, high nitrogen deposition are currently associated with unbalanced foliar N:P ratio (W De [START_REF] Vries | Impacts of acid deposition, ozone exposure and weather conditions on forest ecosystems in Europe: an overview[END_REF], which is widely used as marker of nutrition limitation in trees [START_REF] Güsewell | N: P ratios in terrestrial plants: variation and functional significance[END_REF][START_REF] Tessier | Use of nitrogen to phosphorus ratios in plant tissue as an indicator of nutrient limitation and nitrogen saturation[END_REF]. Deficit or excess in foliar nutrients might thus predispose or amplify climate change effects and biotic damages affecting tree vitality.

Interactions between those multiple driving stress-factors of tree health were rarely tested (but see [START_REF] Seidl | Forest disturbances under climate change[END_REF], yet studies are accumulating that show tree declines under more stressful conditions. To make a step forward in the understanding of these interactions, we conducted a thorough analysis of tree defoliation data collected during a time series of 18 years in a network of forest monitoring plots (belonging to ICP Forests Level II) in France. Focusing on three major forest species in Europe, European beech (Fagus sylvatica), pedunculate oak (Quercus robur) and sessile oak (Quercus petraea), our main objectives were:

1. to characterize the long-term change in tree defoliation; 2. to investigate the interactive effects of local climate, foliar N:P ratio and insect damage on tree crown defoliation.

MATERIAL AND METHODS

Study sites

Tree health data and forest stand characteristics were retrieved from the database of the RENECOFOR network, which is dedicated to the long-term monitoring of forest vitality in France [START_REF] Ulrich | Le réseau renecofor : objectifs et réalisation[END_REF]. The network was set up in 1992 to cover major European tree species growing in temperate forests of France. In this study, we focused on three broadleaved species: sessile oak (Quercus petraea), pedunculate oak (Quercus robur) and European beech (Fagus sylvatica), because they exhibited more defoliation than coniferous species and they were present in a sufficient number of plots to allow reliable statistical analyses (plots characteristics are summarized in Table 1). The monitoring plots were initially established in healthy and preferably even-aged stands dominated by one main tree species [START_REF] Ulrich | RENECOFOR -Manuel de référence n1 pour la sélection et l'implantation des placettes permanentes[END_REF]. In order to have full time series and plots in comparable ecological conditions, we removed from our dataset the plots damaged by the 1999 storms and the beech plots in the few highland sites.

Defoliation

The defoliation assessments started in 1994 on 52 dominants trees of the main tree species in each plot. Then, every summer (from July 1 st to August 31 st ), trained observers recorded defoliation using classes of 5%, in comparison to a healthy referent tree nearby in the same plot, following ICP Forests manual [START_REF] Eichhorn | Criteria for Harmonized Sampling, Assessment, Monitoring and Analysis of the Effects of Air Pollution on Forests[END_REF]. In particular, observers focused on defoliation in the upper and lateral parts of the crown, which were not shaded. Preliminary data analysis revealed for some plots a sharp increase in defoliation from 1994 to 1997. This was attributed to a methodological bias due to a progressive improvement of observers' skills in the early years of the monitoring [START_REF] Ferretti | Plot-scale modelling to detect size, extent, and correlates of changes in tree defoliation in French high forests[END_REF]. We thus decided to analyze tree defoliation from 1997 to 2015. Records from 2003 were missing due to budget shortage.

In addition to defoliation, observers recorded the plausible causes and amount of leaf damage, according to the French Forest Health Department's reference database, in the spring and summer of the same year. Damage assessment protocols have changed over time. For example, the number of possible causes and levels of damages increased with time. To ensure data comparability across time, we decided to base our analyses on the presence/absence of insect damage instead on damage intensity and to focus on the three main causes of damages (see Table A1 for more details on the types of insects causing damage by tree species studied and by year).

To avoid any uncertainty about damage identification (e.g. confusion between insect and fungal pathogen damage on leaves) we defined two categories of trees: those with the "presence of insect damage" when at least one insect damage was reported in the absence of any other type of damage and those with "absence of pest damage" when tree did not present any type of damage or presented only traces of damage (Table 1).

Foliar nutrient analyses

We characterized the nutritional status of trees by considering two major limiting nutrients for tree growth in French forests: nitrogen (N) and phosphorus (P) and their ratio N:P (Table 1). Reference thresholds for tree foliar nutrients in Europe [START_REF] Mellert | Comparison of new foliar nutrient thresholds derived from van den Burg's literature compilation with established central European references[END_REF] were used to define the balanced or unbalanced foliar nutrient status of a tree (Table 1). The lower threshold for foliar N:P ratio was calculated using the lowest threshold of foliar N content and the highest threshold of foliar P content, the upper threshold for foliar N:P ratio was derived from the highest threshold of foliar N content and the lowest threshold of foliar P content. The nutritional status of trees was assessed every second year in all plots. Mature leaves in full sun were sampled in summer on a minimum of 5 dominant or sub-dominant trees in the immediate vicinity of the plot [START_REF] Rautio | Part XII: Sampling and Analysis of Needles and Leaves[END_REF]. Leaves were then dried at 70°C for 24h and ground for chemical analyses. We used a linear interpolation to approximate nitrogen and phosphorus values in the alternate years (when leaves were not sampled).

Table 1 Overview of stand and tree characteristics and range of ecological conditions in the selected plots for pedunculate oak (Quercus robur), sessile oak (Quercus petraea) and European beech (Fagus sylvatica) from 1997 to 2015: average, minimum and maximum values (in brackets) are specified. (-4-5.5) 0.71 (-4.9-5.9) 0.77 (-3.9-5.4)

Drought index

We estimated drought intensity by using the Standardized Precipitation-Evaporation Index (SPEI) [START_REF] Vicente-Serrano | A Multiscalar Drought Index Sensitive to Global Warming: The Standardized Precipitation Evapotranspiration Index[END_REF]. It is extensively used because it can be compared across locations and aggregated at different time scales. Time series of SPEI were calculated on a monthly basis using the difference between the sum of precipitations and the potential evapotranspiration. Precipitations and potential evapotranspiration were obtained from the SAFRAN model of METEO FRANCE at a daily time step for every sampled plot from 1997 to 2015.

Potential evapotranspiration was estimated using the Penman-Monteith method according to FAO recommendations [START_REF] Smith | Revision of FAO methodologies for crop water requirements[END_REF].

We estimated drought intensity for two distinct periods corresponding to the previous year and the year of defoliation assessment, respectively. We calculated SPEI of the previous year on a timescale of 6 months from May to October, to match with the period of broadleaved tree growth. The length of tree growth period was estimated using leaf phenology of the 36 dominants trees, every year on every plot. We used the date at which budburst occurred in at least 90% of dominants trees as the start of the growth period and the date at which 90% of dominants trees were yellowing as the end of the growth period (Table 1). We calculated SPEI of the year of defoliation assessment on a timescale of 3 months (May to July). July is the month at the end of which most of the plots have been visited for defoliation assessment. Finally, to account for the drought experienced by trees during both the previous year and the year of defoliation assessment, we averaged the two values of SPEI. To facilitate the interpretation of analyses outcomes, we used minus SPEI ('-SPEI') instead of 'SPEI' as drought index, so that a high and positive value of SPEI indicates an intense drought (Table 1).

Plot Basal Area

A complete dendrometric inventory of all trees with a diameter at breast height (DBH) above 5 cm within the inner part of the plots (0.5 ha) was made every 5 years as well as before and after logging. Species identity, dominance status and circumference at 1.30 m were recorded for each tree. We then used the basal area (m 2 /ha) of all dominant trees as an index of stand density (Table 1). Missing values (during the 4 years period in between two inventories) were interpolated assuming a linear increase in basal area with time.

Statistical analyses

Statistical analyses were performed with R 3.4.3. We used the log-transformed median of the 5% class of defoliation to normalize the model residuals. The biotic and abiotic determinants of crown defoliation were investigated with a Linear Mixed-Effects Models (function lme package nlme). We considered the nested sampling design and the repeated measurements of trees over time by using a random effect of trees nested in plots. Residual variance was modelled with a power function of fitted values and an autoregressive moving average process estimating the residual temporal autocorrelation [START_REF] Pinheiro | Mixed-effects models in S and S-PLUS[END_REF]. We investigated the long-term temporal trend in tree defoliation by introducing a quadratic effect of the calendar year. Moreover, we tested the interactions between stress drivers by including into the models two-way interactions between drought index, foliar N:P ratio and presence/absence of insect damage. We measured model quality by calculating the root mean square error (RMSE), which is based on the difference between predicted and observed values [START_REF] Janssen | Calibration of process-oriented models[END_REF]. We thus fitted the following model for each tree species: ln(DEFi,j) = β0 + b0,i + b0,i,j + β 1 x YEAR + β 2 x YEAR 2 + β 3 x G + β 4 x INSECT + β 5 x DROUGHT + β 6 x N:P + β 7 x INSECT x DROUGHT + β 8 x INSECT x N:P + β 9 x DROUGHT x N:P + εi,j b0,i ~ Ɲ(0,σ1 2 ) ; b0,i,j ~ Ɲ(0,σ2 2 ) ; εi,j ~ Ɲ(0,σ 2 )

Where DEF is the median level of defoliation classes (in %) of tree i in plot j; β 0: 9 are model parameters for the effects of time (YEAR), stand density (G in m 2 /ha), presence/absence of insect damage (INSECT), drought intensity (DROUGHT) and foliar N:P ratio and their interactions two by two; b0,i is the plot random effect and b0,i,j is the tree random effect. We eventually applied a correction factor when we back-transformed the values fitted by the model in the logarithmic scale to avoid the systematic underestimation of crown defoliation inherent to the log transformation [START_REF] Chave | Improved allometric models to estimate the aboveground biomass of tropical trees[END_REF][START_REF] Finney | On the distribution of a variate whose logarithm is normally distributed[END_REF][START_REF] Shen | Efficient mean estimation in log-normal linear models[END_REF]:

= exp (log ( ) + )
Where is the average crown defoliation that is predicted by the model; log ( ) is the average log-transformed value of crown defoliation that is predicted by the model and is the correction factor:

= 2

Where RSE is the residual square error.

RESULTS

The mean level of defoliation observed was 20%, but the standard deviations associated with random effects in the models indicated important inter-plot and inter-individual variability (Table 2).

Effect of stand density on tree defoliation

Crown defoliation of European beech and pedunculate oak increased significantly (Table 2) with plot basal area (tree density) respectively by 0.4% and 0.3% per m 2 /ha, whereas it decreased significantly by 0.4% per m 2 /ha in sessile oak.

Long term changes in tree defoliation

The three tree species exhibited contrasting temporal trends of tree defoliation over the period 1997-2015 (Table 2, Fig. 1). European beech defoliation increased significantly all along the period (Fig. 1) at a rate of 0.5 % per year. Long-term changes in sessile oak defoliation exhibited a humped shape with an increase from 1997 to 2008 followed by a slight decrease in defoliation from 2008 to 2015. Defoliation of the pedunculate oak showed an opposite trend with a gradual decrease until 2011.

Tree defoliation in presence/absence of insect damage

Trees damaged by pest insects were significantly more defoliated than trees not damaged by insects in the three species (Table 2), at an average value of drought index and foliar N:P nutrient content.

The increase of crown defoliation due to insect damage was of 7.5 % in European beech, 4.7% in sessile oak and 2.6% in pedunculate oak, as estimated by models.

Interaction between insect damage and drought

For European beech, crown defoliation was significantly aggravated by drought, but only in trees damaged by pest insects (Table 2, Fig. 2), i.e. significant drought x insect interaction (estimate=-0.08 ± 0.03; P=0.01). Tree defoliation increased by 5.9% from very wet (drought index =-1.5) to very dry (drought index = 1.5) conditions in trees damaged by insects, as estimated by the model.

Regarding sessile and pedunculate oak, the aggravation of crown defoliation was similar in the presence and in the absence of insect damage (P > 0.05, Table 2, Fig. 2) and was in average of 2.9%

for sessile oak and 3.4% for pedunculate oak between very wet (drought index =-1.5) and very dry (drought index = 1.5) conditions.

Interaction between insect damage and foliar N: P ratio

Crown defoliation of European beech increased with higher foliar N:P ratio in trees that were not damaged by insects (Table 2). The variation in tree defoliation estimated by the model was of +10% between low (16) to high (33) level of foliar N:P ratio (Fig. 3). Conversely, the slope of the relationship between defoliation and foliar N:P ratio was not significantly different from zero for European beech trees damaged by insects, i.e. significant interaction N:P x insect (estimate=-0.045±0.006; P<0.001).

In sessile oak, the effect of foliar N:P ratio was not significant in trees that were not damaged by insects but was negative in trees with an insect damage (-3.5% between low (16) and high (33) level of foliar N:P ratio Table 2), i.e. significant N:P x insect interaction (estimate=-0.007±0.004; P=0.04).

The trends were more pronounced in pedunculate oak: crown defoliation decreased significantly with higher values of foliar N:P ratio with a steeper slope in trees damaged by insects than in trees not damage by insects (-5.6 % for tree damage by insect and -2.6% between low (16) and high (33) level of foliar N:P ratio), i.e. significant N:P x Insect interaction (estimate=-0.008±0.004; P=0.03).

Interaction between drought and foliar N: P ratio

In pedunculate oak, there was a significant interaction between drought and foliar N:P ratio (estimate=-0.008±0.002; P<0.001, Table 2). Trees were more defoliated under dryer conditions and at lower levels of foliar N:P ratio (Fig. 4). The interaction between drought and foliar N:P ratio was not significant in European beech and sessile oak (P>0.05).

Table 2 Summary of the effects of year, stand density (Basal area), presence/absence of insect damage, drought index and foliar N:P ratio on crown defoliation (log transformed) at the tree level in European beech (Fagus sylvatica), sessile oak (Quercus petraea) and pedunculate oak (Quercus robur). Model estimates (est.) are given with associated standard deviation (sd) and p-values (pval). Significant effects (pval<0.05) are in bold. Random effects (sd) correspond to the standard deviation of the random effect at the plot level and the random effect of trees nested within plots, RMSE corresponds to the root mean square error of the linear mixed models. Long-term changes in defoliation that have been observed at the European scale are low or nonexistent for European beech and temperate hardwood oaks (Michel et al., 2018). However, in France, until the end of the 2000s, oaks and beeches seemed to show an increase in defoliation [START_REF] Ferretti | Plot-scale modelling to detect size, extent, and correlates of changes in tree defoliation in French high forests[END_REF]. Despite a high variability between plots and between trees, our study updated the trends at the French level and confirmed the increase in beech defoliation. We also showed a stabilization of defoliation for both oaks more in line with what was observed at the European level. Furthermore, we found that tree defoliation was in general exacerbated by the presence of insect damage and drought events. The effect of the tree nutritional status was more species specific. Interestingly, we also detected interactive effects between these factors, which lead to additive or synergistic effects of stresses for most tree species.

Long term change in tree defoliation

We found a long-term trend of declining crown condition in European beech forests, as tree defoliation was consistently increasing with time. This is consistent with a previous study conducted in the same forest monitoring network in France [START_REF] Ferretti | Plot-scale modelling to detect size, extent, and correlates of changes in tree defoliation in French high forests[END_REF], which reported a similar increase during the period 1997-2009. Recently, Sousa-Silva et al. ( 2018) observed a similar pattern in Belgian forests from 1997 to 2015 and suggested that the decline in European beech condition was related to a progressive increase in drought intensity over time. European beech is a species well known for its high susceptibility to dry conditions [START_REF] Rubio-Cuadrado | Long-term impacts of drought on growth and forest dynamics in a temperate beech-oak-birch forest[END_REF] with strong negative correlations between drought intensity and growth loss [START_REF] Farahat | Growth-climate relationship of European beech at its northern distribution limit[END_REF][START_REF] Lebourgeois | Climate-tree-growth relationships of European beech (Fagus sylvatica L.) in the French Permanent Plot Network (RENECOFOR)[END_REF][START_REF] Seynave | Is the spatial distribution of European beech ( Fagus sylvatica L.) limited by its potential height growth?[END_REF]) and most species distribution models have predicted European beech decline in the future under dryer conditions (e.g. [START_REF] Saltré | How climate, migration ability and habitat fragmentation affect the projected future distribution of European beech[END_REF].

The temporal pattern was quite different in oak species. Sessile oak exhibited a deterioration trend approximatively until 2007, which was followed by a period of slight recovery, while pedunculate oak showed a general trend for a slow crown condition improvement. Oak species are more drought tolerant than European beech [START_REF] Rubio-Cuadrado | Long-term impacts of drought on growth and forest dynamics in a temperate beech-oak-birch forest[END_REF]. It is thus possible that in the recent years, from 2007 onward, precipitations were below the tolerance threshold of European beech but above the tolerance thresholds of oaks. Higher CO2 concentration causes an increase in carbohydrates and a reduction in N, thus increasing C:N ratio which is unfavorable to leaf feeders. It can also enhance C-based tree defenses against defoliating herbivores [START_REF] Robinson | A meta-analytical review of the effects of elevated CO2 on plant-arthropod interactions highlights the importance of interacting environmental and biological variables: Tansley review[END_REF] and particularly against gypsy moth, a common oak defoliator [START_REF] Foss | Effects of Elevated Co 2 Leaf Diets on Gypsy Moth (Lepidoptera: Lymantriidae) Respiration Rates[END_REF].

The lack of degradation, or even the slight improvement, in oak crown conditions might thus be due to positive effects of CO2 enrichment on tree growth and defense against leaf feeders counterbalancing negative effects of droughts on tree health.

4.2 Effects of Intraspecific competition, insect damage, drought and foliar N:P High tree density emerged as a cause of defoliation in all three species, probably due to increased competition for canopy space [START_REF] Getzin | Scale-dependent competition at the stand level assessed from crown areas[END_REF][START_REF] Hajek | Mechanical abrasion, and not competition for light, is the dominant canopy interaction in a temperate mixed forest[END_REF]. Competing trees may also be more prone to pest attacks, due to trade-offs between growth and defenses [START_REF] Huang | Eyes on the future-evidence for trade-offs between growth, storage and defense in Norway spruce[END_REF].

This calls for new studies targeting the effect of forest thinning on resistance to multiple stresses [START_REF] Jactel | The influences of forest stand management on biotic and abiotic risks of damage[END_REF].

Tree defoliation usually increases with summer drought of the current year [START_REF] Seidling | Signals of summer drought in crown condition data from the German Level I network[END_REF] because leaf abscission allows avoiding disruption of sap flow flux and xylem cavitation [START_REF] Bréda | Temperate forest trees and stands under severe drought: a review of ecophysiological responses, adaptation processes and long-term consequences[END_REF].

Defoliation can also arise from lagged effect of previous summer drought [START_REF] Seidling | Signals of summer drought in crown condition data from the German Level I network[END_REF], which can affect carbohydrates reserves and bud viability [START_REF] Bréda | Temperate forest trees and stands under severe drought: a review of ecophysiological responses, adaptation processes and long-term consequences[END_REF].

The negative effect of foliar N:P ratio on tree defoliation did not exhibit a consistent pattern among the studied tree species contrary to what was observed at the European scale. [START_REF] Veresoglou | Exploring continental-scale stand health-N: P ratio relationships for European forests[END_REF] found a segmented response of tree defoliation to foliar N:P, with a positive relationship between nutrient content and defoliation beyond the breaking point of 14.8 in broadleaved species (corresponding to the minimum level of the foliar N:P ratio observed in our study). We found such a positive correlation (higher defoliation at higher level of N:P ratio, i.e. high level of N and low level of P) only for European beech not damaged by insect. European beech covered a larger range of N:P ratio than oak species for which approximatively 80% of the observations belong to the intermediate group. The increase in defoliation observed in European beech may thus arise from the fact that European beech occurs in sites with a stronger nutritional imbalance of N and P than oaks.

In all three species, insect damaged trees exhibited significantly higher defoliation than non-attacked trees, which confirms the significant impact of insect feeding on crown condition [START_REF] Kozlov | Global patterns in background losses of woody plant foliage to insects: Latitudinal patterns in insect herbivory[END_REF]. We found a higher sensitivity of sessile oak than pedunculate oak to insect attack, which is consistent with recent findings on the sensitivity of sessile oak to the oak processionary moth [START_REF] Damestoy | Tree species identity and forest composition affect the number of oak processionary moth captured in pheromone traps and the intensity of larval defoliation[END_REF]. From a management point of view, the sessile oak appears be more resistant to drought events than the pedunculate oak, suggesting that sessile oak should be favored to adapt oak forests to climate change. However, this choice could be made at the cost of making oak stands more vulnerable to attack by defoliators. This trade-off calls for caution when choosing to favor one tree species over another in a substitution logic and suggests that mixing tree species would be a more effective long-term strategy for dealing with multiple risks [START_REF] Jactel | Tree Diversity Drives Forest Stand Resistance to Natural Disturbances[END_REF][START_REF] Sousa-Silva | Tree diversity mitigates defoliation after a drought-induced tipping point[END_REF]. Both beech and oak species can be infested by many insect defoliators, with for example 41 species of Lepidoptera associated with European beech and 189 with pedunculate oak and sessile oak in Britain [START_REF] Kennedy | The Number of Species of Insects Associated with British Trees: A Re-Analysis[END_REF]. In addition, it has been demonstrated that insect damage can accelerate leaf abscission, thus explaining higher defoliation recorded in early summer on insect damaged trees [START_REF] Zvereva | Effects of herbivory on leaf life span in woody plants: a metaanalysis[END_REF].

Interactions between pest damages, drought and foliar N: P ratio

In all tree species, we could observe significantly higher damage (defoliation) in trees that were both attacked by insects and under drought conditions. Interestingly, in our study, we found that the combined effects of drought and insect attack were additive for oaks but synergistic for European beech (Fig. 2). This suggests that, in oak trees, insect damages and drought cumulated their effect on tree defoliation whereas, in beech trees, insect damages increased with water stress intensity. This is consistent with the general pattern of higher forest insect damage following droughts [START_REF] Ferretti | Plot-scale modelling to detect size, extent, and correlates of changes in tree defoliation in French high forests[END_REF][START_REF] Logan | Assessing the impacts of global warming on forest pest dynamics[END_REF][START_REF] Rouault | Effects of drought and heat on forest insect populations in relation to the 2003 drought in Western Europe[END_REF][START_REF] Toïgo | Geographical variation in climatic drivers of the pine processionary moth population dynamics[END_REF] and more particularly the higher susceptibility of water stressed trees to defoliating insects [START_REF] Jactel | Drought effects on damage by forest insects and pathogens: a meta-analysis[END_REF]. The direct effect of drought on forest insect herbivores, e.g. on survival and fecundity, lacks experimental evidence [START_REF] Jactel | Responses of forest insect pests to climate change: not so simple[END_REF]. By contrast, it is well known that drought can indirectly benefit forest defoliators through improved host plant quality (e.g. higher nitrogen content [START_REF] Jactel | Drought effects on damage by forest insects and pathogens: a meta-analysis[END_REF]) or reduced plant resistance (Castagneyrol, Jactel, et al., 2018;[START_REF] Holopainen | Climate Change Effects on Secondary Compounds of Forest Trees in the Northern Hemisphere[END_REF]. It would now be necessary to better identify the species of defoliating insects attacking oaks on the one hand and beeches on the other hand to understand why the synergy between drought and insect attacks was only observed in the case of European beech.

We also showed cumulative, negative effects of pest damage and low N:P ratio on all tree species. This confirms that nutrient imbalance has to be considered when searching the causes of tree crown decline. One possible explanation for the combined effect of pest feeding and low N:P is that insect herbivory is higher on leaf with low N content as more leaf tissue have to be consumed to meet the nutritional requirement of insect herbivores (i.e. compensatory feeding, (Castagneyrol, Moreira, et al., 2018)). In one species, beech, we detected a synergistic effect of insect damage and low N:P leaf ratio, with no effect of nutrient imbalance in undamaged trees and large defoliation increase with increasing N:P in damage trees. This might be due to specific nutrient requirements of Orchestes fagi, the main insect defoliator of European beech in our sample.

PERSPECTIVES AND LIMITATIONS

Our study has several limitations. Tree crown conditions, and in particular level of defoliation, are difficult to estimate accurately. However, observers involved in ICP forests network were trained each year with inter-calibration workshops. In addition, a selection of plots (16 out of 102) between 2001-2005 and 2008-2009 was controlled to assess the quality of the observations of crown defoliation [START_REF] Ferretti | Plot-scale modelling to detect size, extent, and correlates of changes in tree defoliation in French high forests[END_REF]. The observation bias (the difference in score between the field team and the control team) was between ± 5% for 76% of the defoliation scores and ± 10% for 92% of them. This confirms the reliability of crown defoliation score attributed on the field. In addition, we focused on the three tree species with largest sample size and only retained trees for which all necessary data were available. We therefore believe that our results on long-term trend and explanatory factors are conservative. More problematic is the lack of precise information on the insect species that made the observed leaf damage. This would have required collecting branches and assessing a sample of leaves, which is technically and economically difficult. Nevertheless, a more focused long-term study on a few well-known and easily monitored insect defoliators would provide more insights into physiological processed involved in drought, nutrient imbalance and pest damage. Likewise, other potential interactive drivers of crown conditions like storm, air pollution or presence of mycorrhizae would be worth considering.

CONCLUSIONS

Our findings confirm the gradual worsening of crown condition in European beech and the stabilization in oaks, which has been observed at the European scale, suggesting species-specific response to ongoing climate changes. They also demonstrate that events of crown defoliation result from the complex interplay between biotic and abiotic stresses, which are more or less relevant to particular tree species. In particular, insect attacks deteriorated crown condition, especially after dry periods for the three broadleaved species while the effect of imbalanced nutritional status was less consistent among tree species. More importantly, we showed that these stress factors generally act synergistically. It therefore seems important for the adaptive management of forests to address issues related to forest conditions in an integrative way, especially in a context where biotic and abiotic hazards are likely to intensify. Long-term monitoring of forest condition remains more useful than ever, but it should be accompanied by a more accurate analysis of the causes of decline, in order to better predict future trends and better adapt forests to future threats.

Fig. 1 .Fig. 2 .Fig. 3 .Fig. 4 .

 1234 Fig. 1. Temporal change in percentage of defoliation from 1997 to 2015 in European beech, sessile oak and pedunculate oak. Boxplots represent the distribution of raw data. Solid lines represent the crown defoliation of a tree predicted by models after a back-transformation of the logarithmic values at an average level of stand density, N:P ratio and drought index. Shaded areas represent the confidence intervals at 95%

ACKNOWLEDGMENTS

This work and MT's grant were financed by the French National Forest Office (ONF). We are grateful to Frédéric Delport, Guy Landmann and Thierry Belouard for their advice on preliminary results and to Sebastien Macé and Fabien Caroulle for they in-depth methodological advice.