Plants belonging to 12 aspen (P. tremula) genotypes were chosen from the SwAsp collection (genotypes number 4, 6, 7, 26, 41, 51, 69, 72, 79, 92, 98, and 100; Robinson et al., 2012). The selected clones represented the tremuloides-like chemotype, with four dominating SPGs (salicin, tremuloidin, salicortin, and tremulacin) (Robinson et al., 2012; Keefover-Ring et al., 2014). Aspen plants were produced from in vitro tissue cultures (Umeå Plant Science Centre, Umeå, Sweden) and planted into 5-L pots on May 23, 2014. The pots were placed in the greenhouse 20°C, 60% R.H., and a 16:8 h D/N cycle. Lateral branches were removed within the first 4 weeks after planting to promote apical growth.

Two individual plants from each of the 12 genotypes were randomly assigned to one of four treatments: control (C), nitrogen fertilization (N), herbivory (H), and their combination (NH) (n = 24 plants per treatment). The whole experiment thus comprised a total of 96 plants. For N treatment, a Weibulls Rika-S^® solution (containing 84 gr/liter nitrogen, NH₄NO₃) was applied weekly starting from June 17th for 4 weeks. This resulted in a final N input that corresponds to the level of industrial forest fertilization in Sweden, 150 kg N ha^–1 y^–1 (Nohrstedt, 2001; Decker et al., 2017). The first fully expanded leaf on each plant was marked on June 17th, and this leaf was later used to assess baseline leaf chemistry (see below).

After 4 weeks of fertilization treatment, a mousseline fabric net was placed over the first three fully expanded leaves that had developed above the marked leaf after the first fertilization event. Five adult aspen leaf beetles (Chrysomela tremula Fabricius) were placed into the nets after the last fertilization event and allowed to feed for 5 days (treatments H and NH). The beetles had been reared in the laboratory at Umeå Plant Science Centre since the beginning of May and belonged to the second generation of a culture captured close to Ekebo, in southern Sweden. Empty nets were placed in similar position also on plants in C and N treatments to control for the effect of the mousseline bag.

After 5 days of herbivory treatment, the experiment was ended. Leaves were collected from each plant for phytochemical analysis (the leaf below the netted leaves) and for fungal endophyte analysis (the netted leaves). Leaves for the chemical analyses were flash-frozen in liquid nitrogen, freeze-dried, and stored at −20°C until they were ground to fine powder using a bead mill (MM 301 Vibration Mill, Retsch GmbH and Co., KG, Haan, Germany) at 25 Hz for 3 min. The fine powder was stored in vials at −20°C until chemical analyses. Leaves designated for endophyte analysis were collected and stored at 4°C for a maximum of 3 days prior to surface sterilization and endophyte isolation. In case the petiole of a netted leaf was chewed by the beetle (causing wilting or the death of the leaf), those leaves were discarded from further analyses. A total of 90 plants were included in the endophyte analysis. The growth of plants and insects was recorded in the end of the experiment (Supplementary Table 1).

For the global GC/MS analysis, 6.00 (±1.00) mg of leaf powder was extracted in 1 ml methanol:chloroform:water (v:v:v) at 4°C. Deuterated salicylic acid [²H₆] (Isotec, Miamisburg, OH, United States) was used as an internal standard in all samples. Samples were centrifuged at 4°C, and 100 μl of the resulting supernatant was evaporated in vacuum. The residues were resolved in 30 μl of methoxyamine (15 μg/μl in pyridine) and 30 μl of MSTFA. Methyl stearate (30 μl of 15 ng/μl in heptane) was added before analysis, and 1 μl of each aliquot was injected by a CTC Combi Pal autosampler (CTC Analytics AG, Switzerland) into an Agilent 6890 gas chromatograph. The chromatograph was equipped with a 10 m × 0.18 mm fused silica capillary column with a chemically bonded 0.18-μm DB 5-MS UI stationary phase (J&W Scientific, Corston, Bath, United Kingdom). The compound detection was performed in a Pegasus III time-of-flight mass spectrometer, GC/TOFMS (Leco Corp., St Joseph, MI, United States). MATLAB™ R2011b (Mathworks, Natick, MA, United States) was used to quantify the mass by means of integrated peak areas. All pretreatment data procedures, such as baseline correction, chromatogram alignment, data compression, and hierarchical multivariate curve resolution (H-MCR), were performed using custom scripts according to Jonsson et al. (2005). The extracted mass spectra were identified by comparisons of their retention index and mass spectra with libraries of retention time indices and mass spectra (Schauer et al., 2005). Identification of compounds was based on comparison with mass spectra libraries (in-house database) as well as the Kovats retention index. The identified chemicals were quantified by peak area and assigned to metabolite classes (phenolics, amino acids, fatty acids).

For the targeted LC/MS analysis, samples were prepared as above, but the residue was suspended in 25 μl of methanol and 25 μl of 0.1% v/v aqueous formic acid. A 2.0-μl aliquot was injected onto a C18 UPLC™ column (2.1 × 100 mm, 1.7 μm), and chromatographic separation was performed on a LCT Premier TOF/MS in negative mode (Waters, Milford, MA, United States), following the method by Abreu et al. (2011). The standard compounds used in this analysis were from our in-house library (the Swedish Metabolomics Centre: SMC, Umeå, Sweden) including salicinoid standards (salicin, tremulacin, salicortin, and tremuloidin). MassHunter™ Qualitative Analysis software package (version B06.00, Agilent Technologies Inc., Santa Clara, CA, United States) was used to acquire the mass feature extraction (MFE), and extracted features were aligned and matched between samples using Mass Profiler Professional™ 12.5 (Agilent Technologies Inc., Santa Clara, CA, United States). The concentrations of salicinoids were quantified according to the peak area of each compound using linear standard curves. Masses of either or both of the deprotonated ion ([M−H]−) and the formate adduct ([M−H + FA]−) were assessed based on molecular weights according to Abreu et al. (2011) and Keefover-Ring et al. (2014) and guided by retention times where available.

Condensed tannin concentrations were assessed based on a modified Porter’s assay (acid: butanol method) (Porter et al., 1985; Bandau et al., 2015): 10 ± 1 mg fine powder of freeze-dried leaf sample was added to 800 μl of an acetone/ascorbic acid solution (70% acetone, 30% Milli-Q water, with 10 mM ascorbic acid) then incubated with an iron and acid-butanol reagent in boiling water for 1 h. Absorbance was measured on a Spectra Max 190 microplate reader (Molecular Devices, Sunnyvale, CA, United States) at A 550 nm. Procyanidin B2 (C₃₀H₂₆O₁₂, Sigma- Aldrich^®, St. Louis, MO, United States) was used as the tannin standard.

Endophytic fungi were isolated following the method described by Albrectsen et al. (2010). From each leaf, 10 segments (about cm²) were cut using a sterilized scalpel. The segments were surface-sterilized using 95% ethanol and 2% sodium hypochlorite, and rinsed in sterile water before placing them individually on potato dextrose agar (PDA, Merck 70139-2.5KG) in petri dishes (12 × 12 cm) that had been autoclaved and received added streptomycin to avoid bacterial growth. The leaf segments on the plates were checked for emerging fungal growth every second day, and emerging colonies transferred to new PDA-dishes. Three leaves were harvested from every C and N treated aspen plant; and one to three leaves were harvested from the plants that had been exposed to herbivores. The total number of leaves (and petri dishes) was 228.

The recovered fungal cultures (n = 1,642) were classified into 30 distinct morphotypes (MTs) based on their visual morphological traits and growth rate estimated on a four-step scale (fast, medium, slow, very slow). In order to obtain information about the taxonomy of the fungal community, a total of 59 cultures representing the 30 MTs were selected for sequencing. From each culture, plugs of PDA containing growing mycelium were transferred to malt extract broth. After an incubation period of 2 weeks (room temperature, darkness), mycelia were filtered, placed in 50-ml Falcon tubes, and lyophilized for 48 h. Freeze-dried samples were then pulverized in a FastPrep^®-24 homogenizer (MP biomedicals, Santa Ana, CA, United States). Total genomic DNA was extracted using the E.Z.N.A. SP Plant DNA Kit (Omega Bio-Tek, Inc., Norcross, GA, United States). The extracted nuclear DNA was measured using NanoDrop^® ND-1000 (Wilmington, United States). PCR amplification of the internal transcribed spacer (ITS) region of the rDNA was performed with the primers ITS1 and ITS4 (White et al., 1990). Each 50-μl PCR reaction mixture contained 5 μl of 10x PCR buffer, 0.4 μM of each primer, 0.2 mM dNTP, 1.5 mM MgCl₂, 1U Taq polymerase, and 10 ng fungal DNA. The PCR program consisted of 94°C for 5 min, 33 cycles of 94°C for 30 s, 55°C for 30 s, 72°C for 30 s, and 72°C for 10 min. PCR samples were kept at 4°C before analyzing via gel electrophoresis on 1.5% Agarose (BIO-RAD) gel and visualized under UV light with 0.1% GelRed (GelRed™ Nucleic Acid Gel Stain, 10,000X in DMSO). PCR products were purified with the HT ExoSAP-IT High-throughput PCR Product Cleanup (Affymetrix, Santa Clara, CA, United States) following the manufacturer’s instruction. After quantification of DNA concentrations with Qubit fluorometer 3.0 (Life Technologies, Carlsbad, CA, United States), samples were sent to the National Genomics Infrastructure (NGI) at Science for Life Laboratory (SciLifeLab, Uppsala, Sweden). Visualization of raw sequence data was done with the software Chromas (version 2.4.4, Technelysium, South Brisbane, QLD, Australia). All sequences were aligned and manually edited using BioEdit (Ibis Biosciences, Carlsbad, CA, United States). Sequencing data were blasted against best matches in the reference database at NCBI.¹ The ITS sequence homology for delimiting fungal taxa was set to >98.5% for presumed species (Supplementary Table 2).

Kruskal–Wallis non-parametric tests (Table 1) were used to assess the significance of effects on specific chemicals and the summary of the statistical analyses conducted for those relationships that are presented in Figure 2. PCA plots were performed for 114 targeted and non-targeted metabolites from LC/MS and GC/MS and for 11 targeted phenolics (mainly SPGs from LC/MS), respectively, using R version 3.6.1 (2019-07-05) packages: ggplot2,² FactorMineR³ (Lê et al., 2008), and factoextra.⁴

Fungal diversity (richness, S; Whittaker, 1972) was tracked as the number of distinct morphotypes (MTs) found in each plant, and the MT abundance (frequency) was recorded as the number of isolates of a given MT per sampled plant (n = 24 plants for C and N, and 21 plants for H and NH treatments). Venn diagram was constructed to illustrate the number of shared and unique MTs per treatment (Shade and Handelsman, 2012). Relative abundance was determined for each MT as the number of isolates of the given morphotype divided by the total number of isolates emerged from all samples (n = 1642). To examine the possible shifts in the community composition among the different treatments (C, N, H, and NH), a permutational multivariate analysis of variance (Permanova) of relative abundances was performed with the adonis function in vegan package⁵ using 999 permutations, followed by a pairwise comparison with pairwise.adonis (Martinez Arbizu, 2019).

Morphotype richness values were used to create sample-size-based rarefaction (interpolation) and extrapolation (prediction) curves (Colwell et al., 2012; Chao et al., 2014; Hsieh et al., 2016) with an endpoint of 42 individuals and 100 bootstrap repetitions. The curves were generated with the iNEXT (iNterpolation and EXTrapolation) R package⁶ and visualized with ggiNEXT, the ggplot2 extension for iNEXT. To compare the complexity of communities, Simpson’s (D) and Shannon (H′) diversity indices and Pielou’s index for evenness (J′) were calculated using the Vegan package (Oksanen et al., 2016; Supplementary Table 3). The effect of the treatments (C, N, H, NH) on the ecological diversity indexes and on the fungal richness (S) was assessed with an analysis of variance followed by a multiple comparison test (residuals were checked for normality, homoscedasticity, and linearity). When the data did not meet these requirements, a robust ANOVA based on trimmed means was performed with the WRS2 package.⁷ Before all the analyses were performed, data were balanced using a decision tree to generate new samples from the minority classes (H and NH treatments) with an oversampling algorithm (random walk oversampling) from the imbalance package in R⁸ (R Development Core Team, 2016).

The score plot (PCA) of global metabolite analysis (data included 114 metabolites) showed a weak separation between fertilized and non-fertilized treatments (Figure 1A). The plot of targeted metabolites (11 phenolics, including catechin but excluding condensed tannins) did not show any clear distinction between fertilized and non-fertilized plants (Figure 1B).

Nitrogen treatment reduced the concentration of the CT precursor, catechin, in aspen leaves, and a decreasing (non-significant) trend was observed also in the CT concentration (Table 1 and Figures 2A,B). The lowest within-treatment variation in CTs was found in N-treated plants (Table 1 and Figure 2E), despite the varied genetic background of the trees (12 clones). Overall, the effect of N on chemical characters appeared stronger compared to the effect of H (Table 1 and Figures 2A–F), and there were no significant N × H interactions.

The concentrations of amino acids increased in plants receiving N treatment (Table 1 and Figures 2A,B). Most sugars were somewhat reduced by H and increased by N treatment (data for glucose are shown as an example; Table 1 and Figure 2C). A deviating pattern was found only for maltose, which was reduced by treatments, especially by N, alone and in combination with H (Table 1 and Figure 2D).

Nitrogen treatment had a strong effect on both the above-ground biomass and height of the trees [Supplementary Table 1; ANOVA: F_{(df = 3)} = 44.07; p < 0.0001, F_{(df = 3)} = 15.96; p < 0.0001, respectively], and there was no effect of herbivory (p > 0.9).

Among the 30 MTs, 18 were identified to class (2 MTs), genus (10 MTs), or species (6 MTs) level (Table 1 and Supplementary Table 2). Of the identified isolates, 17 belonged to Ascomycota (class Dothideomycetes, Eurotiomycetes, Hypocreales, and Sordariomycetes) and one to Basidiomycota (MT25: Rhodotorula sp.). The three most common morphotypes across the treatments, Ramularia sp. (MT1), and the unidentified morphotypes MT2 and MT3 had 576 isolates (35% of all isolates), 336 (20% of all isolates), and 329 (20% of all isolates), respectively. A total of 18 MTs were captured as a maximum of three isolates (<0.2% relative abundance) and were thus considered as rare MTs.

The sample size-based rarefaction and extrapolation curves (Figure 3A) are predicted to accurately estimate the number of MTs per sampling unit (tree) within the double number of species when compared to the reference sample size (Chao et al., 2014). At least two fungal isolates were recovered from each plant, and the highest increase in diversity was obtained when increasing the number of sampled trees from 1 to 5 within all the treatments. The non-asymptotic rarefaction curves suggest that additional sampling efforts would yield more MTs, especially in H treatment. However, the coverage-based rarefaction and extrapolation curves (Figure 3B) suggested that the sample size (24 trees for C and N each, and 21 for H and NH) allowed high completeness with the method that was used: the sample coverage values were 0.95, 0.97, 0.91, and 0.93 for C, N, H, and NH, respectively.

The treatments caused significant differences in endophyte richness (robust ANOVA, F = 3.70, p = 0.02) and relative abundance (Permanova, F = 1.67, p = 0.04). These differences were located between N and H treatments (p = 0.01 for richness; and F = 2.92, p = 0.02, for abundance), while the rest of combinations did not result in significant differences. Together, the isolates belonging to the three dominating MTs made up 81% of all isolates in C plants, but their proportion tended to decrease by treatments: in N and H treatments, they made up 73% of all isolates, and in NH 71%. At the same time, the relative abundance of the other (less frequent, and rare) morphotypes increased correspondingly (Table 1).

The Shannon index (H’) did not differ among the treatments (ANOVA, F = 1.54, p = 0.21), with mean (±se) and SE values 1.35 (±0.07) for C; 1.50 (±0.06) for N, 1.24 (±0.10) for H, and 1.36 (±0.08) for NH. Likewise, the Simpson index (D) values did not vary among treatments (0.68 ± 0.02 for C; 0.72 ± 0.02 for N; 0.62 ± 0.04 for H; and 0.68 ± 0.03 for NH; robust ANOVA, F = 0.41, p = 0.74) nor did the values of Pielou’s measure of evenness (J) (0.85 ± 0.02 for C; 0.84 ± 0.02 for N, 0.82 ± 0.03 for H and 0.86 ± 0.02 for NH; robust ANOVA, F = 1.03, p = 0.39).

In total, 11 MTs were shared among all four treatments (Figure 4), of which six morphotypes were identified at least to the genus level (Ramularia sp., Cladosporium sp., Penicillium olssonii, Arthrinium sp. and 2, Penicillium sp.) and five remained unknown (Table 2). Three rare MTs were recovered only from C plants, including a Dothideomycetes and a Coleosporium species, and three others from N-treated plants (a Dothideomycetes, Physalospora scripi, and an unknown species). Three rare MTs, including the identified Aureobasidium species (A. microstictum and A. pullulans), were found only in H treatments (Table 2).