Microbiome Variation Across Two Hemlock Species With Hemlock Woolly Adelgid Infestation.

The hemlock woolly adelgid (Adelges tsugae, HWA), an invasive insect, is devastating native hemlock populations in eastern North America, and management outcomes have so far had limited success. While many plant microbiomes influence and even support plant immune responses to insect herbivory, relatively little is known about the hemlock microbiome and its interactions with pathogens or herbivores such as HWA. Using 16S rRNA and ITS gene amplicon sequencing, we characterized the needle, branch, root, and rhizosphere microbiome of two hemlock species, Tsuga canadensis and T. sieboldii, that displayed low and high levels of HWA populations. We found that both archaeal/bacterial and fungal needle communities, as well as the archaeal/bacterial branch and root communities, varied in composition in both hemlock species relative to HWA population levels. While host species and plant-associated habitats explained a greater proportion of the variance in the microbiome than did HWA population level, high HWA populations were associated with enrichment of 100 likely fungal pathogen sequence variants across the four plant-associated habitats (e.g., needle, branch, root, rhizosphere) compared to trees with lower HWA populations. This work contributes to a growing body of literature linking plant pathogens and pests with the changes in the associated plant microbiome and host health. Furthermore, this work demonstrates the need to further investigate plant microbiome effects across multiple plant tissues to understand their influences on host health.

Introduction

A growing body of literature recognizes that microorganisms living inside or in close association with plant tissues are integral to plant health and survival (Compant et al., 2005; Santoyo et al., 2016). In some cases, microorganisms can increase their hosts’ resistance to insect herbivory (Pineda et al., 2017) by affecting plant secondary metabolism (Badri et al., 2013; Hubbard et al., 2019). Plant inoculation with foliar fungal isolates has been shown to reduce herbivory by virtue of fungal metabolites toxic to insects (Tibbets and Faeth, 1999), by fungi acting directly on herbivores as insect pathogens (Marcelino et al., 2008), or by “priming” production of salicylic and jasmonic acids used in plant resistance to pests and pathogens (Thaler et al., 2010). However, the extent and mechanisms of microbiome-induced plant pathogen or herbivore resistance are not broadly understood because these services are primarily studied in model plants and important agricultural species such as Arabidopsis thaliana (Badri et al., 2013), Gossypium (Karban et al., 1987), and Allium cepa (Muvea et al., 2014) and less often in trees such as Populus (Busby et al., 2013). Interestingly, research also points toward the influence of herbivorous arthropods (mites) on the leaf endosphere microbiome and in particular the fungal pathogens Melampsora (Busby et al., 2016) and Drepanopeziza (Busby et al., 2019).

Expanding our understanding of the reciprocal influences of insect and arthropod herbivores and plant host microbiomes could be particularly useful in instances where plants that are especially important to ecosystem health are under threat. For example, eastern hemlock (Tsuga canadensis) is a foundational species in eastern North American forests (Ellison et al., 2005), yet comparatively little is known about the hemlock microbiome and its interactions with pathogens or herbivores, such as the hemlock woolly adelgid (HWA, Adelges tsugae), which is currently devastating native hemlock populations (Eschtruth et al., 2013). The HWA often feeds on young hemlock branches where needles intersect the branch (McClure, 1987), however, the HWA can also be found on the hemlock trunk with unknown consequences for the tree (Oten, 2011; Leppanen et al., 2019). Feeding at the needle base prevents nutrients required for growth from reaching the needles, causing them to discolor and desiccate (Young et al., 1995; McClure et al., 1996; Havill et al., 2016b). The HWA does not appear to harm hemlock species within its native range of Asia and northwestern North America (Oten et al., 2014). However, in the mid 20th century, the HWA arrived in the eastern United States (USA) with the introduction of ornamental hemlocks, and it has since spread from northern Georgia to southern Nova Scotia (Kantola et al., 2019). Hemlock loss can have ecosystem-level effects owing to their foundational role in some eastern mixed hardwood forests. For instance, they provide habitat for many animals (Yamasaki et al., 2000), moderate diel fluctuations in temperature and moisture that improves stream habitats for many invertebrates (Snyder et al., 2002), and slow biogeochemical cycling, preventing stream eutrophication (Jenkins et al., 1999).

The use of chemical control to manage HWA is effective (Silcox, 2002) but not sustainable, and biological control has not yet proven successful in lowering hemlock mortality (Havill et al., 2016b). Resistance or tolerance to HWA in hemlocks from within the native range of the HWA and in some apparently resistant stands in its invasive range is also studied to inform HWA control (Oten et al., 2014; Leppanen et al., 2019; Kinahan et al., 2020). Natural enemies are hypothesized to be at least partially responsible for controlling HWA populations in its native range (Cheah et al., 2004). However, hemlock species from the native range of the HWA (i.e., Asia) introduced to North America (e.g., T. chinensis, T. dumosa) support similarly low HWA populations in eastern North America where these same predators are absent (Bentz et al., 2002, 2007; Leppanen et al., 2019), suggesting a bottom-up resistance to HWA in some hemlock species. This apparent resistance may be conferred through differences in twig tissue chemistry (McKenzie et al., 2014) or cuticle thickness (Oten et al., 2012). Another possibility is that resistance to insect herbivory may, in part, originate from the plant microbiome, as has been demonstrated in some plants (e.g., Mejía et al., 2014; Garbelotto et al., 2019; Hubbard et al., 2019).

Initial investigations of the hemlock microbiome show that the branch microbiome varies across hemlock species, differing between HWA-susceptible and HWA-resistant species (Rogers et al., 2018). However, among HWA-susceptible hemlock species, the microbiome did not differ significantly between HWA population levels. Although these observations suggest that HWA infestation is independent of the plant microbiome, this initial work was limited in replication (n = 3) and investigated only the branch microbiome. It is also possible for outcomes of interactions with the microbiome associated with pest populations to appear in tissues away from the feeding site. For example, the soil microbiome can influence plant secondary metabolism impacting resistance to herbivory (Hubbard et al., 2019). Furthermore, even if the microbiome does not influence HWA populations, HWA feeding and subsequent associated damage still may affect the hemlock microbiome, e.g., because HWA infestation causes a plant immune response and the release of methyl salicylate into the vascular tissue (Pezet et al., 2013). In the rhizosphere of Populus trichocarpa, the concentration of salicylic acid correlated with the abundance of many bacterial and fungal phyla (Veach et al., 2019). Hence, a more systemic evaluation of the hemlock microbiome associated with HWA infestations and resistance is needed to reveal potentially important interactions.

To determine associations between the hemlock microbiome and HWA, we investigated the microbiome of two hemlock species, T. canadensis and T. sieboldii, with different HWA population levels across three plant tissue endospheres (e.g., needle, branch, and root) and their rhizosphere soils. Collectively, we use the term “plant-associated habitats” to describe the plant tissue endospheres and rhizosphere. Tsuga canadensis is native to eastern North America, and T. sieboldii is native to southern Japan (within the native range of HWA, Havill et al., 2016a) but has been introduced throughout the eastern USA (Farjon, 2010). We hypothesize that microbial α-diversity and community composition will differ among plant-associated habitats and between host species as has been shown previously in hemlock (Rogers et al., 2018) and other tree species such as Populus (Cregger et al., 2018), Ginkgo (Leff et al., 2015) and Broussonetia (Chen et al., 2020). However, we also hypothesize that microbial α-diversity and community composition will also differ across HWA population levels. We hypothesize that these differences in microbial community composition will correlate with changes in plant and soil chemistry associated with different host species and HWA population levels. We are also interested in the differences in specific microbial taxa with different hemlock host species and HWA population levels, specifically potential fungal pathogens and mycorrhizal fungi, which are well-characterized in the literature (Nguyen et al., 2016) and are known to affect plant host survivability (Smith and Read, 2008; Dean et al., 2012). We thus hypothesize that the relative abundance of potential fungal pathogens will increase with high HWA population levels due to compromised host defenses (Pezet et al., 2013). Additionally, the relative abundance of mycorrhizal fungi will decrease and the composition of mycorrhizal fungi will differ with high HWA population levels owing to altered resource allocation belowground with HWA infection (Gehring and Whitham, 1994a, b). Our overall goal is to describe the hemlock microbiome across plant tissues and host species and to identify microbial taxa associated with different HWA population levels that might subsequently be considered in HWA control.

Materials and Methods

Site Description and Sample Collection

Hemlock samples were collected on June 23, 2018, at the North Carolina State University Mountain Horticultural Crops Research Station, Mills River, NC, United States (35.420468 N, −82.556092 E, altitude: 643 m). Here, a variety of hemlock species were planted in a mixed-use forested landscape in 2008. Soils are characterized as Hayesville series (clayey, kaolinitic, mesic Typic Kanhapludults). For 2018 and 2019, mean annual precipitation was 176 cm, and temperature between October–March was 7.6°C and April–September was 19.6°C.

We collected samples from 40 hemlock trees across the two hemlock species (T. canadensis and T. sieboldii) and two HWA population levels (high and low) in full factorial design (10 replicate trees per species-HWA population level combination). Trees were characterized as having low HWA population levels when fewer than 10 HWA ovisacs were detected during 10 min censuses of the entire tree, and trees were characterized as having high HWA population levels when ovisacs were detected on >16 of 20 surveyed branches, five in each cardinal direction from approximately 0.1, 0.5, 1.0, 1.5, and 2.0 m above the ground. From each tree, we collected three 10 cm terminal branches from northeast, south–southeast, and west–northwest facing foliage at 1.5 m in height. These samples were composited and frozen on dry ice (−80°C) until sample pre-processing in the laboratory. We also collected fine roots (<2 mm diameter) and the attached soil, which we operationally defined as the rhizosphere, at each tree. Root and rhizosphere collections occurred in the upper 10 cm of soil and within 1 m of the base of the tree. All roots were traced back to the base of the tree. These samples were also frozen on dry ice until sample pre-processing.

Sample Pre-processing and DNA Extraction

Prior to DNA extraction, needles, branches, and fine roots were washed and surface-sterilized as described by Cregger et al. (2018). To increase DNA yield prior to extraction, 50 mg of tissue per sample were cut into ∼5 mm pieces, flash-frozen in liquid nitrogen, and homogenized by bead-beating with a sterile 6 mm steel bead for two 1-minute intervals. Samples went through an additional flash freeze between intervals to prevent thawing. The DNA extractions were performed using the Qiagen PowerPlant Pro DNA Kit (Qiagen, Venlo, Netherlands) following the standard protocol with a slight procedural modification to ensure high-quality, high-concentration DNA yields. This modification consisted of homogenizing in a Precellys 24 (OMNI International, Kennesaw, Georgia, United States) at 3200 g for 3 min at 30 s intervals of pulse and rest. Rhizosphere soil was collected as the pre-sterilized rinsate of the fine roots. Rinsates were centrifuged at 10,000 rcf, and we removed the supernatant. We then used the Qiagen PowerSoil DNA Kit (Qiagen, Venlo, Netherlands) to extract these rinsates, following the standard protocol with the same modification to the procedure as seen above. Extractions were quantified on a NanoDrop 1000 spectrophotometer (NanoDrop Products, Wilmington, DE, United States). We used a Zymo DNA Clean and Concentrator-5 kit (Zymo Research Corporation, Irvine, CA, United States) to purify and concentrate needle, branch, and fine root endosphere extractions prior to polymerase chain reaction (PCR) amplification.

PCR Amplification, Sequencing, and Bioinformatics

Archaeal/bacterial libraries were prepped for 16S rRNA gene sequencing by means of a two-step polymerase chain reaction (PCR) approach with a mixture of custom 515F and 806R primers (Cregger et al., 2018; Rogers et al., 2018), and for fungi using the ITS2 gene region with a custom mixture of primers (Cregger et al., 2018; Rogers et al., 2018; Supplementary Table S1). An adapter sequence was added to each forward and reverse primer to make them compatible with Nextera XT indexes (Illumina). The initial polymerase chain reaction (PCR) consisted of 2 × KAPA HiFi HotStart ReadyMix Taq (Roche, Indianapolis, Indiana, United States), 10 μmol/L total for each forward primer combination, and 10 μmol/L total for each reverse primer combination, with approximately 25 ng DNA. The 16S rRNA and ITS2 PCRs were performed separately. Both reactions consisted of 3 min at 95°C, followed by 25 cycles of 95°C for 30 s, 55°C for 30 s, and 72°C for 30 s, with a final extension at 72°C for 5 min. Successful PCR amplification was confirmed by running 4 μL of PCR product on a 2% agarose gel. The PCR product was then purified by use of AMPure XP beads (Agencourt, Beverly, MA, United States). Nextera XT indexes were then added to the PCR products by use of a second, reduced cycle PCR, such that each sample had a unique combination of forward and reverse indexes. This reduced reaction was the same as the previous reaction but with only eight cycles. The products were purified again using AMPure XP beads. Samples were quantified on a NanoDrop spectrophotometer (Fisher Scientific) and pooled into an archaeal/bacterial pool and a fungal pool to approximately equal concentrations within each pool. Final product sizes and concentrations were confirmed on an Agilent Bioanalyzer (Santa Clara, CA, United States) using the standard sensitivity kit. Both bacterial and fungal libraries were diluted to 4 μmol/L, independently combined with 5% of a 4 μmol/L PhiX adapter-ligated library control, and run paired-end on a v2, 500 cycle flow cell of an Illumina MiSeq sequencer.

Demultiplexed sequences were imported into the QIIME2 environment (Bolyen et al., 2019), and the median Phred quality scores of joined sequences were visualized. Both 16S and ITS2 datasets were denoised and classified into sequence variants (SVs) with the DADA2 algorithm in QIIME2 with reads truncated to 200 bases with the first 25 bases trimmed for 16S and reads truncated to 230 bases with the first 13 bases trimmed for ITS (Callahan et al., 2016). We then assigned representative sequences a taxonomic classification using Naïve Bayes classifier through the sklearn python package for 16S rRNA sequences with the SILVA database (Release 132; Quast et al., 2013), and we assigned taxonomic classifications to ITS rRNA representative sequences using BLAST and the UNITE reference database (version 8.0, Abarenkov et al., 2010). We removed contaminants (unassigned reads, mitochondria, chloroplasts for 16S; Protista, Chromista, Animalia, and Plantae reads for ITS2).

Soil and Plant Chemical Analysis

To determine correlations between microbial community composition and plant and soil chemistry, branch, root, and rhizosphere samples were sent to the University of Georgia Extension Soil, Plant, and Water Laboratory for chemical analyses. Branch and root tissues as well as rhizosphere soils were ground and analyzed for total carbon (C) and nitrogen (N) concentrations by direct combustion using the Elementar vario MAX CNS Element Analyzer (Elementar, Langenselbold, Germany). Additionally, rhizosphere soils were analyzed for pH and lime buffer capacity (LBC). Briefly, pH was measured in a well-mixed 1:1w:v soil:CaCl2 slurry (0.01 M) using a Fisherbrand accuTupH Rugged Double Junction pH Combination Electrode (Waltham, MA, United States). For LBC, pH was measured before, and 30 min after, a 2.7 ml addition of 0.023 M Ca(OH)2 to a 20 g soil and 20 ml 0.01 M CaCl2 slurry using the same pH electrode as above following Kissel et al. (2012).

Statistical Analysis

All statistical analyses were conducted in R (R Development Core Team, 2008) with the phyloseq (McMurdie and Holmes, 2013), DESeq2 (Love et al., 2014), hillR (Li, 2018), nlme (Pinheiro et al., 2017), and vegan (Oksanen et al., 2019) packages.

Differences in α-diversity were compared by means of Hill numbers (Jost, 2006) of the point estimate of samples rarified to 1,000 reads (highest number of reads present across all samples) at orders of q = 0 and q = 1 (full rarefaction curves are presented in Supplementary Figure S1). The parameter q determines the relative weighting of rare species. At q = 0, all species are weighted equally (richness); at q = 1, species are weighted proportionally to their relative abundance (analogous to Shannon’s index). Differences in means of Hill numbers among plant-associated habitats, host species, and HWA population levels were assessed by nested ANOVA with tree identity as a random effect. Because we were primarily interested in a HWA population-level effect, where we found significant interactions between HWA population level and host species and/or plant-associated habitat, we performed individual ANOVAs and corrected the p-values using the Benjamini and Hochberg false discovery rate adjustment (Benjamini and Hochberg, 1995). The resulting ANOVAs did not include a random effect if plant-associated habitat was not a dependent variable because the resulting models would have had one sample per tree. Where independent variables were significant, we assessed multiple comparisons by Tukey’s test of honest significant differences. We used Q-Q plots and scale-location plots to inspect normality and homoscedasticity, respectively.

Differences in the community composition of the archaeal/bacterial and fungal microbiomes among plant-associated habitats, between hemlock host species, and across HWA population levels were assessed by nested distance-based redundancy analysis (dbRDA) constraining permutations within individual tree. We used the varpart() function in vegan (Oksanen et al., 2019) to determine the variance explained by each factor in our dbRDAs. For the dbRDAs, we used quantitative Jaccard (Ružička) distances applied to proportionally normalized data. Similar to our approach for α-diversity, where we found significant interactions, we performed individual dbRDAs for each host species or plant-associated habitat and corrected the p-values as described above. When performing separate dbRDAs for each plant-associated habitat, we similarly did not constrain permutations by tree because we had only one sample per tree.

We also assessed differences in community composition associated with HWA population level by identifying differentially abundant SVs across HWA population levels in each plant-associated habitat across host species. To do this, we first normalized the SV table through variance stabilization, then estimated the fold change of differentially abundant microbial SVs between low and high HWA population levels using Wald tests and shrinkage estimation for dispersions (Love et al., 2014) and similarly adjusted p-values of differentially abundant SVs with the Benjamini and Hochberg false discovery rate adjustment (Benjamini and Hochberg, 1995).

To assess the variation of the microbial community explained solely by plant and soil chemistry, we conducted dbRDAs for each plant-associated habitat individually (and therefore did not have to constrain by tree) using only the plant and soil variables as independent variables in the models. We conducted separate dbRDAs for each plant-associated habitat because we used only proximal plant chemical data (e.g., we did not attempt to correlate root C:N with the branch microbiome because we had chemical data from the tree branch). We determined the variation explained by each predictor with variance partitioning using the varpart() function in vegan (Oksanen et al., 2019).

Differences in fungal potential pathogen and mycorrhizal relative abundance (assessed by FUNGuild, Nguyen et al., 2016) among plant-associated habitats, host species, and HWA population levels were assessed by nested ANOVA and Tukey’s test of honest significant differences similarly to our approach for α-diversity. The resulting models were similarly inspected for normality and homoscedasticity. To satisfy these assumptions, the dependent variable in each model was log-transformed. We also assessed differences in the ectomycorrhizal community composition of the root and rhizosphere separately between host species and across HWA population levels using dbRDA, similarly using quantitative Jaccard (Ružička) distances applied to proportionally normalized data.

Results

Sequencing Results

After quality and taxonomic filtering (i.e., removal of plant and plasmid DNA), we sequenced 6.30 × 106 16S reads across 142 samples [18 samples were removed due to low read depths (<1,000)], with a minimum read depth of 1,448 and a maximum of 244,262. For ITS, we sequenced 5.92 × 106 reads across 160 samples with a minimum read depth of 1,586 and a maximum of 122,712.

Alpha Diversity

Archaeal/Bacterial Community

Archaeal/bacterial α-diversity (at q = 0 and 1) differed across plant-associated habitats (q = 0: F3,88 = 1047.358, p < 0.001; q = 1: F3,88 = 744.668, p < 0.001; Figure 1) with the microbiome of the rhizosphere being more diverse than the plant tissue microbiomes (all comparisons: p < 0.001). Also, at q = 0 (i.e., richness), the branch microbiome was less rich than the needle and root microbiomes (both comparisons: p < 0.01). We also detected greater archaeal/bacterial α-diversity in T. canadensis compared to T. sieboldii, but only at q = 0 (q = 0: F1,37 = 8.539, p = 0.006; q = 1: F1,37 = 1.482, p = 0.231). We failed to detect an effect of HWA population on archaeal/bacterial α-diversity (q = 0: F1,37 = 0.964, p = 0.333; q = 1: F1,37 = 0.759, p = 0.389).

FIGURE 1

Boxplots representing α-diversity based on Hill numbers (Jost, 2006) of archaea/bacteria and fungi across plant-associated habitats, hemlock woolly adelgid (HWA) population levels, and host species at q = 0 (richness) (A) and q = 1 (analogous to Shannon diversity) (B). Note different axis scales.

Fungal Community

Plant-associated habitat and host species interacted in their effect on fungal α-diversity (q = 0: F3,107 = 6.031, p < 0.001; q = 1: F3,107 = 4.023, p = 0.009; Figure 1). Therefore, we analyzed the differences in fungal α-diversity for each plant-associated habitat individually. At q = 0, host species and HWA population level interacted in their effect on needle and root fungal richness (needle: F1,36 = 4.422, p = 0.043; root: F1,36 = 5.457, p = 0.025) such that diversity in T. canadensis exceeded that in T. sieboldii only at low HWA population levels (needle-low: p < 0.001, needle-high: p = 0.386, root-low: p < 0.001, root-high: p = 0.486). At q = 1, T. canadensis had greater needle and root fungal α-diversity compared to T. sieboldii regardless of HWA population level (needle: F1,36 = 7.505, p = 0.010; root: F1,36 = 3.906, p = 0.056). We detected no effect of host species or HWA population level on branch or rhizosphere fungal α-diversity at q = 0 and q = 1 (all: p > 0.1).

Microbial Community Composition

Plant-associated habitat explained 21.3% of the variation in archaea/bacteria community composition (p < 0.001, adj-R2 = 0.203). Because there was a three-way interaction among plant-associated habitat, host species, and HWA population level (p = 0.068), we analyzed each plant-associated habitat separately.

Hemlock woolly adelgid population level explained 2 and 1% of the variation in the needle and branch archaea/bacteria microbiomes, respectively, across both hemlock species (needle: p = 0.002, adj-R2 = 0.019; branch: p = 0.089, adj-R2 = 0.008; Figure 2). Additionally, there was a host species∗HWA population-level interaction in the root archaea/bacteria microbiome (p = 0.011) such that there was a greater HWA population-level effect in T. sieboldii (p = 0.036) than in T. canadensis (p = 0.052). At the phylum level, needles on trees with high HWA populations had greater abundance of Actinobacteria and lower abundance of Proteobacteria compared to needles on trees with low HWA populations, and the branch microbiome of trees with high HWA populations had greater abundance of Bacteroidetes and lower abundance of Actinobacteria compared to the branch microbiomes of trees with low HWA populations (Supplementary Figure S2). At the order level, high HWA population levels corresponded with high levels of Cytophagales in the needle microbiome and high levels of Betaproteobacteriales and Sphingomonadales in the branch microbiome (Supplementary Figure S3). For the needle, branch, and root archaea/bacteria microbiomes, there was a relatively stronger main effect of host species (needle: p < 0.001, adj-R2 = 0.048; branch: p < 0.001, adj-R2 = 0.074; root: p < 0.001, adj-R2 = 0.032), however, differences did not clearly emerge at the phylum level. Instead, these differences emerged at the family and genus level. For instance, we found that T. sieboldii had greater relative abundance of Beijerinckiaceae but a lower relative abundance of the genus Candidatus Uzinura (order: Flavobacteriales) in the branch microbiome compared to T. canadensis (Supplementary Figures S4, S5). We detected no effect of either HWA population level (p = 0.894) or host species (p = 0.182) on the rhizosphere archaea/bacteria microbiome composition. The effect of HWA population level on the microbial community also emerged at the sequence variant (SV) level. Hemlock woolly adelgid population level was associated with four differentially abundant archaeal and 1,057 differentially abundant bacterial SVs across the four plant-associated habitats (some SVs are shared among plant-associated habitats; needle: 104 SVs, 6.4% of SV richness; branch 8 SVs, 1.2% of SV richness; root: 173 SVs, 3.8% of SV richness; rhizosphere: 791 SVs, 2.7% of SV richness; Supplementary Figure S6 and Supplementary Table S2).

FIGURE 2

Distance-based redundancy analysis (dbRDA) ordinations of archaea/bacteria community composition across plant-associated habitats, hemlock woolly adelgid (HWA) population levels, and host species. Note different axis scales.

Plant-associated habitat explained 13% of the variation in fungal community composition (p < 0.001). Because of a three-way interaction among plant-associated habitat, host species, and HWA population level on the fungal microbiome (p < 0.001), we also analyzed each plant-associated habitat separately.

Hemlock woolly adelgid population level explained 1% of the variation in the needle fungal community composition in both host species (p = 0.008, adj-R2 = 0.013, Figure 3). We detected no association between HWA population level and the composition of branch (p = 0.472), root (p = 0.174), and rhizosphere (p = 0.924) fungal communities (Figure 3). Except for the rhizosphere, composition of all other plant microbiomes was influenced by host species (needle: p < 0.001, R2 = 0.082; branch: p < 0.001, R2 = 0.109; root: p < 0.001, adj-R2 = 0.009; rhizosphere: p = 0.204). These differences in host species emerged at the class level with greater relative abundance of Teliomycetes (particularly order Helotiales) in the needles and branches and of Dothideomycetes (particularly order Pleosporales) in the roots of T. canadensis compared to those of T. sieboldii (Supplementary Figures S7, S8). Differences at the family and genus level were more nuanced because taxonomic classification at these levels is for the most part incomplete (Supplementary Figures S9, S10). Hemlock woolly adelgid population level was also associated with 1,481 differentially abundant fungal SVs across the four plant-associated habitats (some SVs are shared among plant-associated habitats; needle: 583 SVs, 19.1% of SV richness; branch: 168 SVs, 17.3% of SV richness; root: 298 SVs, 27.5% of SV richness; rhizosphere: 665 SVs, 19.5% of SV richness; Supplementary Figure S11 and Supplementary Table S3).

FIGURE 3

Distance-based redundancy analysis (dbRDA) ordinations of fungal community composition across plant-associated habitats, hemlock woolly adelgid (HWA) population levels, and host species. Note different axis scales.

Correlation of Microbiomes With Soil and Habitat Characteristics

Total soil C, total soil N, branch C:N, root C:N, pH, and LBC were correlated with the composition of the needle archaeal/bacterial community (p = 0.006, variance explained: 17.3%), the root archaeal/bacterial community (p = 0.003, variance explained: 17.8%), and the needle fungal community (p = 0.046, variance explained: 15.7%; Figure 4). We detected no correlations between the plant and soil chemistry data and the microbial community for all other microbial community × plant-associated habitat combinations (p > 0.1). For the needle archaeal/bacterial community, lime buffer capacity (LBC) explained 1% of the variation in community composition (p = 0.014, adj-R2 = 0.011), and branch C:N explained 2% of the variation in community composition (p = 0.078, adj-R2 = 0.016). Root C:N and soil pH explained 2% and 3%, respectively, of the root archaeal/bacterial microbiome composition (root C:N: p = 0.012, adj-R2 = 0.025, pH: p = 0.026, adj-R2 = 0.020). Branch C:N and LBC explained 2 and 1%, respectively, of the needle fungal microbiome composition (Branch C:N: p = 0.072, adj-R2 = 0.015, LBC: p = 0.073, adj-R2 = 0.006).

FIGURE 4

Distance-based redundancy analysis (dbRDA) ordinations of the needle archaeal/bacterial (16S) community composition, root archaeal/bacterial community composition, and needle fungal (ITS) community composition across hemlock woolly adelgid (HWA) population levels and host species. Microbial community compositions were ordinated along the variables soil total carbon (C), total nitrogen (N), pH (1:2 CaCl2), and lime buffering capacity LBC as well as the C:N of the plant-associated habitat (rhizosphere, root, or branch [same for needle]). Other environmental variable dbRDAs were not significant (p < 0.05). Environmental variables are represented by arrows, and bolded labels represent significant (p < 0.05) variables in the dbRDAs. Note different axis scales.

Fungal Potential Pathogens and Mycorrhizal Fungi

One hundred fungal SVs classified as potential pathogens across the four plant-associated habitats were associated with high HWA populations, and about half of these were found in aboveground plant tissues (Supplementary Table S3). Analyzing the relative abundance of fungal potential pathogens in our samples, we detected a three-way interaction among plant-associated habitat, host species, and HWA population level (F3,107 = 2.801, p = 0.044, Figure 5). Therefore, we analyzed each plant-associated habitat separately. When analyzed separately for each plant-associated habitat, the relative abundance of fungal potential pathogens was comparable overall among HWA population levels and host species, except in specific instances. For example, there was an almost 10-fold greater relative abundance of potential pathogens in the roots of T. sieboldii with a high HWA population level compared to T. sieboldii with a low HWA population level (p = 0.029) (the relative abundance of potential pathogens in T. canadensis roots did not vary among HWA population levels [p = 0.988]). Between host species, the root fungal microbiome had a 5-fold greater relative abundance of potential pathogens in T. canadensis than in T. sieboldii only with low HWA population levels (p = 0.049).

FIGURE 5

Relative abundance of fungal potential pathogen reads as a proportion of all fungal reads across plant-associated habitats, hemlock woolly adelgid (HWA) population levels, and host species.

The relative abundance of ectomycorrhizal (EM) fungi in the root endosphere in T. canadensis significantly exceeded that in T. sieboldii (p = 0.025, Figure 6). However, there was no host species effect on the relative abundance of EM fungi in the rhizosphere (p = 0.693). Similarly, HWA population level was not associated with EM fungal relative abundance in the roots (p = 0.297) or rhizosphere (p = 0.930). The EM community composition did not vary between HWA population levels and host species in both roots (HWA: p = 0.736; host species: p = 0.281) and rhizosphere samples (HWA: p = 0.778; host species: p = 0.107; Supplementary Figure S12).

FIGURE 6

Relative abundance of arbuscular mycorrhizal (AM) and ectomycorrhizal (EM) fungal reads as a proportion of all fungal reads in the roots and rhizosphere across hemlock woolly adelgid (HWA) population levels and host species. Note different axis scales.

Discussion

Consistent with our hypothesis, HWA population level was associated with many specific microbial taxa in the microbiomes of T. canadensis and T. sieboldii across multiple plant tissues and the rhizosphere at the SV level. Such findings, however, are inconsistent with previous research that found no association of HWA population levels with the branch microbiome (Rogers et al., 2018). By increasing the sample size compared to that of Rogers et al. (2018) (10 vs. 3), increasing the scope of the sampling to include other plant-associated habitats, and by investigating HWA-hemlock microbiome associations at multiple scales (e.g., SV level, community level), we were able to detect a significant relationship between HWA population level and the hemlock microbiome.

At the community-level, we detected a significant HWA population level association with the hemlock needle microbiome for both archaea/bacteria and fungi. It is not surprising that the needle microbiome had the strongest association with HWA population level because HWA infestation can affect nutrient delivery to the needles (Havill et al., 2016b). However, we found little effect of altered nutrient status on the needle microbiome in our environmental variable dbRDA, likely because we did not measure needle C or N and used branch C and N instead as a proxy. Future work should measure the nutrient content of the needles, including micronutrients, which may affect microbial community composition as well (Kembel et al., 2014), to test the hypothesis that HWA-induced changes in nutrient content affect the needle microbiome.

Infestation of HWA may also affect plant performance by increasing plant susceptibility to pathogens, either by compromising the plant defense system (e.g., Pezet et al., 2013) or by increasing labile substrate in the affected plant tissues (e.g., Tooker and De Moraes, 2009). Specifically, we found about an 8-fold enrichment of two Gibberella spp. SVs (which could not be classified to species resolution) in the needle microbiome of trees with high HWA populations (Supplementary Figure S11 and Supplementary Table S3). Gibberella species are globally widespread plant pathogens associated with many plant hosts, and they have multiple modes of pathogenesis (Desjardins, 2003). In an agricultural study, Gibberella ear rot severity in corn (Zea mays) was linked with the western bean cutworm (Striacosta albicosta) infestation (Parker et al., 2017), highlighting the interaction between plant pests and fungal pathogens. These Gibberella spp. SVs and the other 98 SVs classified as potential pathogens that were associated with high HWA population levels should be prioritized for future study of the interaction between HWA infestation and the hemlock microbiome.

We present preliminary evidence that hemlocks with high HWA population levels are selecting for microorganisms that may improve plant defense. For instance, we detected an eight-fold (on average) enrichment of two Mycobacterium and two Pseudomonas SVs (Supplementary Figure S6 and Supplementary Table S2) in the root endosphere; in some cases, these are known to produce salicylic acid (Ratledge and Winder, 1962; Visca et al., 1993; Lemanceau et al., 2017). Salicylic acid is an important plant defense compound (Pieterse et al., 2012), and by selecting for microorganisms with the capability to produce salicylic acid, plants may be better equipped to defend against pathogenesis (Lebeis et al., 2015). However, such evidence is highly speculative, and further metabolomic and transcriptomic work is necessary to determine if these taxa increase salicylic acid production in HWA-infested plants.

Lack of a mycorrhizal response to HWA population level is surprising in light of the fact that HWA not only reduces photosynthetic capacity (Nelson et al., 2014) and presumably C allocation belowground but also increases the nutrient supply in litter through increased throughfall (Stadler et al., 2006), both of which can decrease mycorrhizal colonization (Gehring and Whitham, 1994a, b). Resistance to HWA could be supported through mycorrhizal networks where mycorrhizae colonize multiple trees (Simard et al., 2012), altered growth strategies (e.g., mycorrhizal to saprotrophic, Johnson et al., 1997), or a delayed signal from the plant. Also, it was interesting that arbuscular mycorrhizal (AM) fungi were in such high abundance in hemlock roots, especially those of T. canadensis with low HWA population levels. Many of these AM fungi could not be identified beyond the family level (Glomeraceae). Hemlock species (family: Pinaceae) are not normally associated with AM fungi (Smith and Read, 2008), but in a greenhouse bioassay experiment, 25% of T. heterophylla seedlings were colonized by AM fungi (Cázares and Smith, 1995). Such findings counter the traditional paradigm that members of Pinaceae associate exclusively with EM fungi and promote the idea of mycorrhizal co-occurrence in Pinaceae (Wagg et al., 2008). Because the relative dominance of mycorrhizal types can potentially affect ecosystem-level processes (Phillips et al., 2013), the impact of HWA infestation on co-occurrence of AM and EM fungi in hemlock warrants detailed research.

Consistent with earlier work (Rogers et al., 2018), we found a greater percent of the variation in the microbiome composition explained by host species than by HWA population level. The effect of host species on the microbiome composition was strongest in the needles, branches, and roots, where the plant has a relatively stronger control over the microbiome environment (Kembel et al., 2014). Indeed, root C:N, which was, on average, about 12% lower in T. canadensis compared to T. sieboldii, was a significant determinant of microbiome composition across microbial domains. However, the differences in the microbiome composition among host species generally did not affect the relationship between the microbiome composition and HWA population level (i.e., no host species × HWA population level interaction). Therefore, we conclude that the microbiome compositions of these two HWA-susceptible species correlate with HWA population level in much the same way.

As with other studies (e.g., Beckers et al., 2017; Rossmann et al., 2017; Cregger et al., 2018), we found large differences in the composition of microbial communities among the different plant-associated habitats. Aboveground plant tissues were dominated by Alphaproteobacteria and Ascomycota, specifically two fungal classes: Dothideomycetes and Eurotiomycetes. Roots were dominated by Actinobacteria, and rhizospheres were enriched in Acidobacteria. Belowground habitats also had a greater proportion of Basidiomycota reads, specifically Agaricomycetes. These broad taxonomic patterns among different plant-associated habitats resemble those found in other temperate tree species such as Magnolia kwangtungensis (Qian et al., 2019), Populus trichocarpa and P. deltoides × P. trichocarpa hybrids (Cregger et al., 2018), and Picea abies (Kovalchuk et al., 2018; Ren et al., 2019; Terhonen et al., 2019), suggesting that, at higher taxonomic levels, microbiomes are fairly consistent among tree species. However, as our study and others show, at more specific taxonomic levels, microbiomes diverge among closely related host species (Cregger et al., 2018; Rogers et al., 2018) and even among different genotypes of the same host species (Bálint et al., 2013; Veach et al., 2019).

An important consideration of this work is that these results were obtained during a single sampling date. Indeed, microbiomes change seasonally and interannually (Redford and Fierer, 2009; Shade et al., 2013), and these temporal dynamics of the microbiome may increase or decrease our ability to distinguish ecological phenomena (Grady et al., 2019; Dove et al., 2020). Nevertheless, our results suggest modest to strong variations in the microbiome among HWA population levels, host species, and plant-associated habitats. Future work should determine the temporal robustness of these trends.

Investigating interactions among pests, microbial communities, and plant genetics contributes to a holistic understanding of the plant system that can be leveraged to promote plant health. Using 16S rRNA and ITS gene amplicon sequencing, we found a relatively modest relationship between HWA population level and the hemlock microbiome composition in two species. Nevertheless, even modest dissimilarities in the overall microbiome can result in functional consequences when specific driving taxa are differentially abundant (Agler et al., 2016). Future work should specifically investigate interactions between HWA infestation and the differentially abundant taxa highlighted in this study, especially those classified as potential plant pathogens.

Author’s Note

This manuscript has been authored by UT-Battelle, LLC under Contract No. DE-AC05-00OR22725 with the United States Department of Energy.

Data Availability Statement

The datasets presented in this study can be found in the Sequence Read Archive (https://www.ncbi.nlm.nih.gov/bioproject/622552) under project ID 622552. The SV, taxonomy, and sample data tables used in this analysis can be found in Supplementary Tables S4–S9, and the R code for all statistics and figures can be found at https://github.com/nicholascdove/Hemlock_microbiome.

Author Contributions

MC, JF, CL, TJR, and DS designed the study. MC, AL, CL, TGR, TJR, and DS collected the sample. AL and TGR established and maintained the study site. VB, CL, and TJR processed the sample and analyzed the laboratory. ND led the data analysis. ND and MC wrote the manuscript with critical input from all coauthors. All authors contributed to the article and approved the submitted version.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.