Forest Fire Influence on Tomicus piniperda-Associated Fungal Communities and Phloem Nutrient Availability of Colonized Pinus sylvestris

Sampling Locations

For this study, four 20- to 30-year-old pine forest stands in Sweden were selected to explore the effect of forest fire on phloem nutrient quality and T. piniperda-associated fungal communities, Hälleskogsbrännan (HSB: 59°54′35.5″N, 16°08′42.2″E), Ecopark Öjesjön (OJ: 59°51′04.2″N, 16°14′57.8″E), Knutby (KB: 59°56′46.1″N, 18°18′00.4″E), and Skyttorp (SK: 60°06′32.4″N, 17°49′34.0″E) (Fig. 1). Two of the sites (HSB and OJ) were affected by a single forest fire event during summer 2014, and the other two were not (Fig. 1b). Following forest fire, substrate availability is high, and beetles do not fly far from these abundant food sources [38, 39]. Thus, we expect that the distance between HSB and OJ (approximately 9 km) is enough to assume that the bark beetle populations and associated fungal communities sampled from the two burnt sites are independent of each other.

Sample Collection

The study was initiated in February 2016, at which time we expected T. piniperda to have established in the area and have completed one generation post-fire. Five trees at each site, spaced 30 m apart, were felled in mid-February, before the first T. piniperda flight. The trees were placed on plastic cylinders approximately 20 cm above the forest floor and left without additional treatment (Fig. 1c).

In early March before warming temperatures would initiate beetle flight, three free-hanging window traps were mounted above each of the stems, and two funnel traps were mounted by the felled trees at each site. Traps were baited with 70% ethanol and ( −)-α-pinene (Sigma-Aldrich, Saint Louis, MO, USA) to attract adult T. piniperda to the traps during flight (Fig. 1c). To avoid cross-contamination between beetles, an artificial pinecone made by carbon paper (Fig. S1) was placed in each collection cup so that the beetles could hide from cannibalism and other enemies while in the cup. To monitor beetles flight, traps were continuously checked from March 11–April 6, 2016. When beetles were observed in the cups, the artificial cone was transferred to a watertight plastic bag, and a new artificial cone was placed in the collection cup. Cones were handled separately and transported on ice to the lab where they were stored at 4 °C. Within 24 h, T. piniperda collected from each artificial cone were sorted by sex and dissected under a stereomicroscope. Samples were prepared by placing two individuals of the same sex into a 1.5-mL microcentrifuge tube. To collect spores and fungal propagules from their bodies and gut contents, the beetle samples were washed by vortexing for 30 s in 50 µL RNAlater™ Stabilization Solution (Invitrogen, Carlsbad, CA, USA), and the resulting body wash solution was collected. Gut samples of both individuals were then collected by dissecting the beetles by splitting the back body of the insect using a forceps. The gut was put in RNAlater™ Stabilization Solution (Invitrogen, Carlsbad, CA, USA). Samples were collected to represent all of the traps with beetles from each sampling occasion. Between two and four sampling occasions were needed to reach the target of 24 samples of each sex from each site. In the end, four different beetle sample types were thus collected: female body wash, male body wash, female gut, and male gut, all with 24 replicates per site for a total of 384 beetle samples. All samples were stored frozen until DNA extraction.

After felling, a phloem sample (pre-colonization time point) was taken from each of the stems using a leather punch (28 mm diameter) and a rubber hammer (Fig. 1d). Bark was removed prior to phloem sampling. Each sample was a composite of five phloem plugs placed in an air-tight plastic bag and transported on ice back to the laboratory where they were frozen at − 20 °C and later freeze-dried (Heto LyoLab 3000 Freeze-Dryer, Thermo Fisher Scientific, Waltham, MA, USA) for at least 24 h.

At each site, three of the five felled trees were chosen in April to be included in the study based on successful beetle colonization during beetle flight. This resulted in a total of 12 felled trees across the four sites. Following T. piniperda colonization, phloem samples were collected from each tree at two time points corresponding to two larval developmental stages. Two control samples were taken from areas ≥ 30 cm from the closest larval gallery (Fig. 1e) at both times, and eight phloem samples in galleries were taken in the first instar (mid-April; Fig. 1f), as well as in the third instar (mid-May; Fig. 1g). At both time points, a total of 10 phloem samples were collected from each felled tree (as described previously for pre-colonization samples). Each larval gallery represents a separate brood resulting from one infestation event. This resulted in a total of 252 phloem samples: ten samples (eight colonized and two uncolonized phloem samples) at two post-colonization time points from each of the 12 trees, as well as the sample taken pre-colonization from each tree. Just as with the pre-colonization phloem samples, each sample was a composite of five phloem plugs. All samples were transported on ice to the lab and frozen at − 20 °C within 12 h. Samples were later freeze-dried (Heto LyoLab 3000 Freeze-Dryer, Thermo Fisher Scientific, Waltham, MA, USA) for at least 24 h and stored for DNA extraction.

DNA Isolation, Amplification, and Sequencing

Prior to DNA extraction, T. piniperda samples (body and gut) were homogenized in 2-mL screw cap tubes with three stainless steel hexagon screw M3 nuts (BAHAG AB, Mannheim, Germany) and approximately 40 borosilicate glass balls (1 mm diameter) (Sigma-Aldrich, Saint Louis, MO, USA) at a frequency of 30 Hz for 3 min in a bead beater (TissueLyser II, Qiagen, Hilden, Germany). Phloem samples were homogenized by grinding using a kitchen mixer at 16,000 rpm for 5–10 s (Bamix, Mettlen, Germany). DNA was extracted from the beetle samples using the NucleoMag Plant Kit (Macherey–Nagel, Düren, Germany). DNA extraction was conducted using the same kit for homogenized phloem samples, but with double the amount of lysis buffer due to the high absorption capacity of the phloem samples. The internal transcribed spacer 2 (ITS2) region of nuclear ribosomal DNA was amplified by PCR using the forward primer fITS9 (5′–GAACGCAGCRAAIIGYGA–3′) [40] and reverse primer ITS4 (5′–TCCTCCGCTTATTGATATGC–3′) [41], each with unique barcoded primers. Three parallel PCRs were conducted in 20 µL volumes containing 1.5 U DreamTaq DNA polymerase (Thermo Fisher Scientific, Waltham, MA, USA) with 10 × DreamTaq™ buffer containing 20 mM MgCl₂, 0.2 mM dNTP, 1 mM fITS9 primer and 0.3 mM ITS4 primer, and 10–20 ng genomic DNA template. MgCl₂ concentration was optimized for each sample type to a final concentration of 22.75–29 mM. Reactions were carried out using an Applied Biosystems 2720 Thermal Cycler (Applied Biosystems, Foster City, CA, USA) with the following cycling protocol: an initial denaturation step at 95 °C for 10 min, 30–35 cycles of denaturation at 95 °C for 30 s, annealing at 58 °C for 30 s, and extension at 72 °C for 50 s, followed by a final extension step 72 °C for 3 min. Negative and positive controls were included in each PCR (DNA from an Agaricus bisporus fruitbody was used as the positive control). Products from the three parallel PCRs were pooled and amplification was confirmed via gel electrophoresis using a 1.5% Agarose Basic (PanReacc AppliChem ITW Reagents, Chicago, IL, USA) and GelRed Nucleic Acid Gel Stain (Biotium, Inc., Freemont, CA, USA). DNA from 423 of the 636 samples was successfully amplified. Amplified DNA was cleaned using the AMPure XP kit (Beckman Coulter, Brea, CA, USA) and then quantified using a Quant-iT PicoGreen dsDNA Assay Kit (Life Technologies, CA, USA). Equal amounts of PCR product from each sample, approximately 50 ng, were pooled into one library for body and gut samples and a separate library for phloem samples. The two libraries where vacuum centrifuged at 1300 rpm (Scan Speed 32, Labogene ApS, Lynge, Denmark) until the volume reached 60 µL and sequenced at SciLifeLab/NGI (Uppsala, Sweden) on a PacBio RS II system (Pacific Biosciences, Menlo Park, CA, USA) using two SMRT cells (one for beetle samples, one for phloem samples). Sequence data were delivered to us as error-corrected FASTQ files (containing circular consensus sequencing reads).

Quantification of C, N, and P in Phloem Samples

From each homogenized phloem sample, approximately 12 mg was used to determine total C and N by combustion (Costech Elemental Combustion System 4010, Costech Analytical Technologies, Inc. Italy). Total P was quantified using an ICP AVIO200 (Perkin Elmer, Waltham, MA, USA) in nitric acid extracts following the national standard protocol SS 28,311:2017 (SIS, Stockholm, Sweden). Extraction and quantification of P were conducted in the certified laboratory at the Department of Soil and Environment, Swedish University of Agricultural Sciences, Uppsala, Sweden. Quantification of C and N was successful for 249 of the 252 phloem samples, and P concentration was successfully quantified for all samples. Phloem N, C, and P concentration (mg/kg) and C:N, C:P, and N:P ratio (molar) were calculated for each sample.

Read Quality Filtering and OTU Generation

Sequence reads were quality filtered, demultiplexed, and clustered into Operational Taxonomic Units (OTUs) using the Sequence Clustering and Analysis of Tagged Amplicons (SCATA) pipeline (available at https://scata.mykopat.slu.se/ and first described in [40]. For this we used the following settings: clustering distance of 1.5% sequence similarity, minimum sequencing length threshold of 150 bp, and minimum allowed base quality of 2. After read clustering, all global singletons were removed. For simplicity, we will refer to cluster names “scata3997_number” as “OTU_number” throughout the text. Samples with invalid tagged primer combinations were removed from the resulting OTU/sample matrix, as well as any singletons created by their removal. The resulting OTU/sample matrix consisted of 15,854 reads (412 samples, 362 OTUs). Any OTUs with zero reads in the dataset were also removed from the FASTA file containing the representative sequence for each OTU cluster.

The ITS2 region was then extracted using ITSx (version 1.1.2) [42], and taxonomy prediction was accomplished using the SINTAX classifier (USEARCH v11.0.667; [43] with a bootstrap confidence cut-off of 0.8. The UNITE USEARCH/UTAX dataset (version 8.0, all eukaryotes, https://doi.org/10.15156/BIO/786346, UNITE [44] was used as the reference dataset and was modified following [45] so that it would be compatible with the SINTAX classifier [43]. Sequences for which the ITS2 region was not successfully extracted by ITSx or were not predicted to represent fungi were filtered out of the dataset for downstream analysis. In addition, sequence reads corresponding to A. bisporus (OTU_1) were removed from the dataset, since DNA from an A. bisporus fruitbody was used as the positive control for PCRs and this organism is unlikely to occur naturally in this environment. This resulted in a dataset consisting of 12,313 reads (401 samples, 289 OTUs). Samples with less than 10 reads were also excluded from downstream analyses, as well as any singleton OTUs in the dataset created by the removal of these samples, resulting in 11,962 reads (331 samples, 286 OTUs) used for analysis of fungal communities.

Data Visualization and Multivariate Statistical Analysis

Fungal Communities

Non-metric multidimensional scaling (nMDS) was used to visually explore the relative (dis)similarity between samples based on the fungal community detected in each. Prior to ordination, samples (or subsets of samples) were normalized to relative read abundance, and ordination based on Bray–Curtis dissimilarity was conducted using the vegan R package (version 2.5–7; [46] with R version 3.6.2 [47] using a maximum of 200 random starts and condensed to two dimensions. Ordination was conducted on beetle samples and phloem samples prior to conducting ordination with samples combined. There were 7155 sequence reads (167 samples, 258 OTUs) and 4807 sequence reads (164 samples, 96 OTUs) included in ordination for the beetle and phloem subsets, respectively. (For an overview of the number of samples representing each sample type, see Table S1).

Permutational multivariate analysis of variance (PERMANOVA) tests were used to test for the significance of observed patterns via the “adonis2” function of the vegan R package, and in some cases homogeneity of dispersion was tested for groups using the “betadisper” function coupled with the base R “anova” function. For both beetle and phloem sample datasets, standard PERMANOVA tests were used to test for a significant difference in groups based on fire history, and both marginal and standard PERMANOVA tests were carried out on the full dataset to investigate the relative significance of differences in fungal community composition based on fire history and beetle vs. phloem samples. We also used PERMANOVA to test for significance of groupings based on sample type (male and female, body and gut) for beetle samples and relating to sampling time and bark beetle colonization status for phloem samples. Pairwise PERMANOVAs with a Bonferroni correction was then used to test for significant differences between uncolonized phloem and first instar colonized phloem, as well as first instar vs. third instar colonized phloem samples. Pairwise PERMANOVA with Bonferroni correction was also used to test for effect of fire history within uncolonized samples, within first instar colonized samples, and within third instar colonized samples.

Based on the results of ordination and PERMANOVA tests, we identified eight groups of samples (four categories of samples from both burnt and unburnt sites): beetle (all beetle samples), uncolonized pine phloem (no larvae; pre-colonization samples, and uncolonized control samples from first and third instar sampling time points combined), pine phloem colonized by first instar larvae, and pine phloem colonized by third instar larvae. For an in-depth distribution analysis of the most prominent OTUs across these eight categories, we identified a core community including OTUs with a minimum of 1% average relative read abundance and detection in a minimum of 30% of the samples in any of the eight groups. A heatmap was used to display the shifts in fungal community across the eight sample groups. Color intensity represents the scaled (but not centered) average relative read abundance within each OTU, and the average relative read abundance is reported within each cell.

Taxonomic assignment of the 15 OTUs identified as the core community was manually curated in October 2021 by running BLAST queries to the online UNITE database [48]. A Unite Species Hypothesis (USH) threshold of 1.5% dissimilarity was selected to reflect the clustering used in this study. Species identification was often not possible since the ITS2 region is highly conserved across multiple closely related species for most of the lineages represented in the core community [49]. The curated taxonomy was used for discussing possible functions and lifestyles represented by the core OTUs in this study system (Table S2).

Phloem Chemistry

To analyze phloem chemistry, a partially Bayesian mixed-effects model was fitted for each of the three nutrient ratios (C:N, C:P, N:P) using the blme R package (version 1.0–5; [50]. In each model, fire history (burnt or unburnt site), bark beetle colonization (colonized or uncolonized phloem), and sampling time (first or third instar) were treated as predictor variables and the nutrient ratio was treated as the response variable. Interaction terms for each pair of the three explanatory variables and a three-way interaction were included in the model. The nutrient ratio measured pre-colonization for each tree was treated as a nuisance co-variate, and the interaction between tree and site was treated as a random effects variable. The assumption of equal slopes was tested by comparing the model described here to a model that also includes an interaction between the pre-colonization nutrient ratio and each of fire history, bark beetle colonization, and sampling time. In all cases, the pre-colonization nutrient ratio values were centered (but not scaled), and sum-to-zero contrast coefficients were applied to each of the three explanatory factors. A model including interactions was significantly better for C:N (χ²₍₃₎ = 9.7412, Pr(> χ²) = 0.0209) and for N:P (χ²₍₃₎ = 17.974, Pr(> χ²) = 0.0004) but not for C:P (χ²₍₃₎ = 6.2103, Pr(> χ²) = 0.1018), indicating that the ANCOVA assumption of equal slopes for a co-variate is violated for the C:N and N:P ratio models. The inclusion of interaction terms has little impact in the outcome of the analysis for the C:P ratio model, so the interaction terms were retained for all three models. However, it seems that the interaction between pre-colonization ratio value and sampling time may be driving this, so a model fitted with only the interaction between pre-colonization nutrient ratio and sampling time was fitted to see if the more complex model with all three co-variate–explanatory variable interactions was still significantly better. In all three comparisons, the more complex model was not significantly better (C:N: χ²₍₂₎ = 0.3075, Pr(> χ² = 0.8575; C:P: χ²₍₂₎ = 5.8374, Pr(> χ²) = 0.054; N:P: χ²₍₂₎ = 2.4065, Pr(> χ²) = 0.3002). Therefore, the final models included the three explanatory variables (bark beetle colonization, fire history, and sampling time) and their interactions, pre-colonization nutrient ratio, the interaction between the pre-colonization nutrient ratio and sampling time, and the interaction between site and tree as a random effects variable.

Assumptions of mixed-effects models were inspected both visually and statistically. We tested for homogeneity of variance by plotting the standardized residuals against predicted values. The data for all three nutrient ratios were mildly heteroscedastic, however reasonably so (data not shown). Levene’s tests indicate that the assumption of equal variance of residuals is violated for several grouping factors across the three models (Tables S3–S5). Normality of errors for both the residuals and random effects were inspected visually by histogram for residuals, Q-Q plots for random effects and statistically via Shapiro–Wilk normality tests. Shapiro–Wilk tests indicate that the distributions of residuals deviate from normality for all three nutrient ratios (Tables S3–S5). Histograms show that the distribution looks close to normal for all three nutrient ratios, with the number of values more than two standard deviations from the median within reason given the number of samples (data not shown). Random effects did not have errors that deviated from a normal distribution for any of the three models, as determined via Q-Q plots (data not shown) and Shapiro–Wilk normality tests (Tables S3–S5). Although these datasets may violate some of these model assumptions, a recent study showed that mixed-effects models are quite robust to such violations [51]. We therefore proceeded with analysis.

Type II Wald χ² tests via the car R package (version 3.0–10; [52] were used to determine the marginal significance of each term (Tables 1, 2, and 3. Significance of each term was also evaluated by 95% confidence intervals (CI; Tables 1, 2, and 3). The performance R package (version 0.7.2; [53] was used to determine the marginal and conditional R² calculated using the Nakagawa method, and the marginal and conditional intra-class correlation coefficient (ICC) values for each model. Conditional R² values of 0.731, 0.774, and 0.777 for C:N, C:P, and N:P, respectively, show that the model in each case accounts for between 73 and 78% of the variation in the dataset. Marginal R² values of 0.469, 0.352, and 0.320 for C:N, C:P, and N:P, respectively, indicate that of the variance in the model, 47%, 35%, and 32% is attributable to the fixed effects. All in all, the model is a good fit for the data in each case, although variation between tree and site (random effects) explains a large proportion of the variance. ICC values confirm this observation of large proportions of variance explainable by random effects for C:N (ICC_adj = 0.494, ICC_cond = 0.262), C:P (ICC_adj = 0.651, ICC_cond = 0.422), and N:P (ICC_adj = 0.672, ICC_cond = 0.457).

The results show that for all three nutrient ratios, bark beetle colonization had a significant main effect, but it also had significant interactions with other explanatory variables, indicating that the significance of bark beetle colonization is dependent on the value of these explanatory variables. Fire history had no significant main effect on any of the three nutrient ratios, but it did have significant interactions with other explanatory variables. This indicates that while fire history alone cannot explain much variation on nutrient ratio, fire history can play a role in determining the effects (either in significance or direction) of other explanatory variables. Each nutrient ratio dataset was therefore subset by fire history, and the same model (excluding the fire history term) was used to tease apart these interactions. For each, type II Wald χ² tests and 95% CIs were used to evaluate the significance of the main and interaction effects for each remaining explanatory variable (sampling time and bark beetle colonization). The six subset models also had similar R² and ICC values as the full models (Table S6).