Contrasting microbial assembly patterns in the woody endosphere of hybrid and non-hybrid Populus trees

Sampling design

Twig samples were collected from natural tree stands located across two study sites in southern Quebec: Mont-Saint-Hilaire and Oka National Park (collection approved by the Société des établissements de plein air du Québec; authorisation number: PNO-2020-008), located ∼70 km from each other. At each site, twigs were collected from two host taxa: Populus deltoides and the naturally occurring hybrid, P. × jackii (P. balsamifera × P. deltoides). Five trees of each host were sampled between July 20th and July 28th, 2020, at each of the two sites, resulting in a total of 10 trees per host. Five shaded branches measuring ∼24 inches were excised from the phyllosphere of each tree at a height of ∼10–20 ft from the forest floor. Cuttings were placed in plastic bags and transported to the laboratory on ice, where they were stored at 4 °C until processing.

Sample processing

We processed all stored samples within 48 h of collection. Processing involved several steps: Asymptomatic twigs were first cut from tree branches using sterilised loppers and then rinsed with sterile double-distilled water and several drops of Tween 20 before surface sterilization. The twigs were then submerged and stirred vigorously in 70% ethanol for one minute, followed by 4% sodium hypochlorite for 10 min. We rinsed the samples between wash steps by submerging them in sterilised double-distilled water, with an additional five rinses at the end of sterilisation to remove any residual chemicals. To assess the effectiveness of surface sterilisation, twigs were pressed into culture plates containing potato-dextrose agar and incubated at room temperature for 48 h. The sterilised twig samples were then stored in sterile falcon tubes at −80 °C until DNA extraction. We repeated these steps for each of the 20 twig samples.

DNA extraction, library preparation, & sequencing

To homogenise each of the samples, a sterile scalpel was used to cut approximately 12 randomly sampled sterilised twigs into small pieces between 1.0 to 2.0 mm in size. The cut pieces were mixed, and 150 mg of tissue was collected in tubes for bead beating. Samples were bead-beaten for 90 s prior to DNA extraction. We extracted genomic DNA from samples using a DNeasy Plant Mini Kit (Qiagen) following the manufactures instructions.

To reduce contamination from host DNA, we chose to amplify the V5–V6 region of the 16S rDNA gene using the chloroplast discriminating primers, 799F/1115R. For fungal DNA amplification, we targeted the ITS1 rDNA gene using the universal fungal primers, ITS1F/ITS2. A two-step PCR approach was used to amplify both gene regions. PCR amplification of the V5–V6 region was performed in a 25 µL mixture containing three µL of DNA extract, 0.2 µM of both forward and reverse primer, 0.5 U of Phusion Hot Start II High-Fidelity DNA Polymerase (ThermoFisher), 1x of Phusion HF Buffer, 0.2 µM of dNTPs, and 3% DMSO. We used the same mixture for ITS1 amplification with the exception that four µL of DNA extract was included in the mixture.

PCR cycles for amplification of the V5-V6 region were performed for 30 s at 98 °C for initial denaturation, followed by 15 s at 98 °C, 36 cycles of 30s at 64 °C, 30 s at 72 °C, and 10 min at 72 °C for final extension. The same intermediate steps were used for ITS1 amplification, except 37 cycles were performed at 55 °C. A 2% agarose gel was used to visualise PCR products, which were normalised using a Just-a-Plate Normalisation Kit, following the manufacturer’s protocols (CharmBiotech). Samples were pooled into V5-V6 and ITS1, and V5-V6 samples were purified with AMPure XP beads at a ratio of 0.75:1 (bead-to-DNA ratio; Beckman Coulter). ITS1 pools were purified three times with a 0.8:1 bead to DNA ratio to eliminate primer dimers. A PhiX control library (Illumina) was spiked into the amplicon pool, and Illumina MiSeq sequencing was performed using a MiSeq reagent kit V3 (paired-end 300 bp; Illumina).

Bioinformatics & data processing

We utilised the Quantitative Insights Into Microbial Ecology pipeline (QIIME2 v.2020.11.1) for all data processing steps (Bolyen et al., 2019). Cutadapt was first used to trim primer regions from demultiplexed reads prior to denoising steps (Martin, 2011). Fungal and bacterial-derived datasets were denoised independently through the DADA2 algorithm to generate amplicon sequence variants (ASVs), remove chimeric and low-quality sequences, and merge paired-end sequences (Callahan et al., 2016). We used the default settings of DADA2 with the following exceptions: (1) low-quality tails were first removed from sequences (Q score <30) by trimming the ITS1 forward reads at 276 bp and the reverse reads at 239 bp. The 16S rRNA reads were trimmed at 281 and 217 bp, respectively; (2) We used the pseudo-pooling parameter to increase the detection sensitivity of low-frequency ASVs.

Taxonomic assignment was then performed using the naive Bayes classifier available in the sklearn Python package using the SILVA database for 16S rRNA (Quast et al., 2012) and the UNITE database for ITS1 sequences (Abarenkov et al., 2010). Classified sequences were then filtered to remove contamination from non-target DNA, including chloroplast and plant DNA. Unassigned ASVs were also removed for downstream analysis. We detected three reads belonging to the genus Modestobacter and five reads belonging to an unidentified fungal ASV within bacterial and fungal negative control samples, respectively. The low number of reads detected indicates that our samples were free from contamination. We, therefore, removed the negative controls from subsequent analyses. Finally, we removed potentially spurious low-abundance ASVs that contained fewer than 10 sequences. Although there is no consensus on the most appropriate filtering threshold for phyllosphere amplicon sequencing studies, we chose a conservative filtering approach, as our analysis relied heavily on diversity comparisons, which may be sensitive to low-abundance taxa (Nikodemova et al., 2023). For example, Nikodemova et al. (2023) found that most sporadically detected operational taxonomic units generated through 16S rRNA sequencing occur at less than 10 copies, and removing such microbiota may improve the reliability of our microbial diversity estimates.

Phylogenetic inference

Phylogenetic trees were constructed for bacterial and fungal communities using the FastTree 2 algorithm in QIIME2 with a midpoint root (Price, Dehal & Arkin, 2010); Figs. S1 & S2). However, the high variability of ITS sequences can lead to poor sequence alignments among disparate fungal lineages, resulting in unresolved phylogenies that may bias metrics of community phylogenetic structure (Molina-Venegas & Roquet, 2014). We, therefore, ran the fungal community phylogenetic analysis using two trees: the original tree generated with FastTree 2 and a second tree generated using a ‘ghost tree-like’ method (Fouquier et al., 2016), whereby extension trees inferred from ITS sequences were grafted onto a resolved backbone phylogeny (Fig. S3).

For the ‘ghost tree-like’ method, we used Shen et al.’s (2020) phylogeny of Ascomycota as a backbone tree, accessed through TreeHouse (Steenwyk & Rokas, 2019). To generate extension trees, we first grouped ASVs by Class and generated an ITS sequence alignment for each group using the MUSCLE algorithm (Edgar, 2004). We then generated an extension tree for each alignment, using the IQTREE algorithm (Nguyen et al., 2014) in QIIME2, and grafted the resulting trees onto the backbone tree at their respective crown node (i.e., their Class node). An ITS extension tree was also generated for Basidiomycete sequences and grafted to the root of the backbone tree.

We focused our analysis on the fungal community results obtained with the FastTree 2 method, as we were only able to run the ghost tree analysis with ASVs identified to the class level, which excluded ∼36% of sequences; however, we report the results obtained with the ghost tree (Figs. S4–S7, Tables S1 & S2) and where they differ from those of the FastTree 2 method.

Community analysis

Rarefaction curves were generated for bacterial and fungal samples independently (Fig. S8). To standardise sampling effort, we rarefied bacterial samples to 1,628 reads and fungal samples to 6,342 reads. These values represented the minimum resampling depths that would allow us to maintain all samples for subsequent analysis. To ensure adequate sequencing depth was achieved, we used Good’s coverage index, which indicated that all samples had greater than 99% coverage (Fig. S9). Community analysis was performed with the rarefied dataset using the phyloseq (McMurdie & Holmes, 2013), picante (Kembel et al., 2010), and vegan (Oksanen et al., 2020) packages in R v.4.3.2 (R Core Team, 2023). We performed all community analyses independently on bacterial and fungal community data.

We first explored whether the taxonomic or phylogenetic diversity within samples (alpha diversity) differed between host taxa. Taxonomic measures included ASV richness, Simpson’s dominance, Chao1, and Pilous evenness. We utilised picante (Kembel et al., 2010) to assess the phylogenetic dispersion of microbial communities (i.e., if assemblages consist of taxa that are closer or more distant relatives than expected) by calculating the standardised effect size of mean pairwise distance (ses.MPD) and mean nearest taxon distance (ses.MNTD). Both indices represent complementary measures of dispersion, with ses.MPD placing greater emphasis on deeper branches connecting taxa in the phylogeny (basal patterns) and ses.MNTD shallower branches (terminal patterns; Mazel et al., 2016). We utilised the abundance-weighted versions of ses.MPD (ses.MPD_ab) and ses.MNTD (ses.MNTD_ab), which emphasise microbial taxa that are common within samples and, therefore, account for bias that may result from rare microbial ASVs. Both dispersion measures are independent of the species richness of samples and are expressed as z-scores relative to a null distribution of random community assembly. We computed ses.MPD_ab and ses.MNTD_ab using the null model argument taxa.labels and 1,000 randomisations. To assess potential differences in alpha diversity between hosts, an analysis of variance (ANOVA) was performed.

We then calculated the between-sample (beta diversity) equivalent of ses.MPD_ab and ses.MNTD_ab—that is, ses.βMPD_ab and ses.βMNTD_ab—following methods adapted from Dini-Andreote et al. (2015). ses.βMPD_ab and ses.βMNTD_ab measure the phylogenetic distance separating ASVs in two assemblages (Kembel et al., 2010) and are expressed as z-scores relative to a null distribution of phylogenetic distances generated from random communities. The null distribution thus reflects the expected βMPD_ab or βMNTD_ab between assemblages when stochastic ecological processes dominate community assembly (Dini-Andreote et al., 2015). To calculate ses.βMPD_ab and ses.βMNTD_ab, we first Hellinger transformed the community matrix (Oksanen et al., 2020) and calculated the observed βMPD_ab and βMNTD_ab between samples using the comdistnt and comdist functions, respectively. We then constructed null distributions using 1,000 randomly assembled communities generated by randomising the community matrix and ASV labels on the phylogeny using the randomiseMatrix function (with the null.model argument set to richness) and the tipShuffle function, respectively. Each random matrix was Hellinger transformed before calculating their respective βMPD_ab or βMNTD_ab pairwise distances. ses.βMPD_ab and ses.βMNTD_ab were then calculated as follows: ${\displaystyle ses.\beta ={\beta }_{obs}-{\beta }_{null}/sd\left({\beta }_{null}\right)}$ where ses.β is ses.βMPD_ab or ses.βMNTD_ab, β_obs is the observed pairwise distance between samples (βMPD_ab or βMNTD_ab), and β_null and sd(β _null) are the mean and the standard deviation of the corresponding null distribution of each pairwise comparison, respectively. When ses.βMPD_ab or ses.βMNTD_ab diverge from the mean of the null distribution by more than two standard deviations (z-score <−2 or >+2), they significantly differ from βMPD_ab or βMNTD_ab expected under stochastic community assembly, reflecting the role of deterministic ecological processes. Conversely, ses.βMPD_ab or ses.βMNTD_ab z-scores within two standard deviations of the mean are consistent with stochastic processes dominating community assembly (Stegen et al., 2013; Dini-Andreote et al., 2015). Positive z-scores indicate that assemblages are more phylogenetically dissimilar than expected, reflecting variable selection (z-score >+2), while negative z-scores indicate that assemblages are more similar than expected, reflecting homogenising selection (z-score <−2; see Dini-Andreote et al., 2015).

To assess differences in community composition between hosts, we utilised the abundance-weighted taxonomic measure Bray–Curtis dissimilarity (Bray & Curtis, 1957) and its phylogenetic analog weighted UniFrac (Lozupone & Knight, 2005; Lozupone, Hamady & Knight, 2006). Distance matrices were generated using Hellinger-transformed community data and visualised through principal coordinates analysis (PCoA) to examine differences in the composition of host-associated microbial communities. We performed a permutational analysis of variance (PERMANOVA) on distance matrices to test for differences in community composition between host taxa and sites using the adonis2 function in vegan. We also utilised a permutational analysis of multivariate dispersions (PERMDISP), available through the betadisper function, to test for differences in variance between the microbial communities of host taxa and sites using the same distance matrices (Oksanen et al., 2020).

To test for differences in the abundance of microbial taxa between host taxa (i.e., differential abundance analysis), we used the Analysis of Compositions of Microbiomes with Bias Correction (ANCOM-BC) R package (Lin & Peddada, 2020). Given that ANCOM-BC accounts for differences in sequencing depth, we ran the analysis separately on filtered non-rarefied bacterial and fungal data. We used default arguments with the exception that we set the alpha threshold to 0.001 and struc_zero argument to TRUE. This argument performs a presence-absence test for microbial taxa that are only present within a single host. By default, the ancombc function excludes taxa from the analysis that are detected in fewer than 10% of samples, which for our analysis excluded taxa that were detected in only one sample.

Sequencing & filtering

Illumina sequencing generated ∼1.56 million raw reads (680,153 16S rRNA and 874,920 ITS1 reads). Denoising and quality filtering with DADA2 resulted in the retention of 356,641 and 360,050 reads from 16S rRNA and ITS1 libraries, respectively. 16S rRNA reads contained a high proportion of contamination (50.2% of reads), resulting from the presence of mitochondria and chloroplast DNA, which accounted for 28.5% and 21.7% of reads, respectively. ITS1 reads contained low amounts of Populus DNA contamination (∼0.6% of reads) and lepidopteran DNA belonging to a single ASV (0.02%). Filtering out contaminants, unclassified reads, and low-abundance ASVs (<10 sequences; Nikodemova et al., 2023) resulted in 355,245 fungal and 176,528 bacterial reads, belonging to 227 and 667 ASVs, respectively. The number of sequences per sample ranged from 1,628 to 35,003 for bacterial samples and 6,342 to 39,711 for fungal samples. We rarefied samples to a sequencing depth of 1,628 and 6,342 for bacterial and fungal samples, respectively (Fig. S8). These rarefied data were used for all subsequent steps in our analysis. No differences in the qualitative patterns or significance of our results were detected between resampled runs of the analysis. We present the results of a single run with rarefied data.

Alpha & beta diversity analysis

The bacterial endophytic communities of twigs were, on average, more diverse and exhibited more even species distributions than fungal communities, although no differences in taxonomic diversity were detected between host taxa (Fig. S10). Measures of basal and terminal phylogenetic dispersion indicated that both bacterial and fungal communities trended towards phylogenetic clustering (i.e., negative z-scores), although some individual host plants harboured microbial assemblages with weak to random phylogenetic structure (i.e., z-scores ≈ 0). We detected no significant differences in phylogenetic dispersion between the microbial communities of host taxa, except for the fungal community of P. × jackii, which was less terminally clustered (ses.MNTD_ab) than P. deltoides (Fig. 1).

Phylogenetic beta diversity patterns revealed that microbial assembly was driven by a combination of ecological processes. Stochastic processes tended to be more prevalent in bacterial assembly; however, deterministic processes were more common within the bacterial community of P. × jackii than P. deltoides when measured through ses.βMNTD_ab (53.3% versus 24.5% of pairwise comparisons, respectively; Fig. 2). Conversely, fungal assembly was characterised mainly by deterministic processes, although stochastic processes were more dominant in the fungal community of P. × jackii than P. deltoides when measured through ses.βMNTD_ab (80.0% versus 20.0% of pairwise comparisons, respectively) and ses.βMPD_ab (44.4% versus 4.4%, respectively; Fig. 2). The fungal community of P. deltoides displayed a strong signature of determinism in the form of homogenising selection (80.0% of ses.βMNTD_ab and 95.6% of ses.βMPD_ab pairwise comparisons; Fig. 2).

Visualisation of community composition through PCoA plots revealed that host-associated communities formed distinct clusters within multivariate ordination space through both weighted UniFrac (Fig. 3) and Bray–Curtis distance measures computed on Hellinger-transformed community data (Fig. S11). Results from the PERMANOVA (using the same distance matrices) indicated that variance in the fungal and bacterial community was explained by host identity, while sampling site and the interactions between sampling site and host identity were not significant predictors (Table 1). The fungal community of P. × jackii had greater variance in taxonomic and phylogenetic community composition (i.e., weighted UniFrac and Bray–Curtis distances) than P. deltoides, as indicated by PERMDISP (Fig. 3 & Table 2). However, we did not detect significant differences in variance through PERMDISP analysis with the ghost tree, suggesting that fungal community variance among P. × jackii is primarily driven by ASVs that are not identified to the class level (Table S2). No significant differences in bacterial community variance were detected between host taxa (Fig. 3 & Table 2).

Bioindicator detection

Fungal endophytic samples were dominated by ASVs belonging to Ascomycota (92.4%), with Basidiomycota present at lower abundances (0.2%). The remaining ASVs were unassigned fungal sequences (7.4%). At the class level, Dothideomycetes were the most abundant ASVs (49.1%), followed by Eurotiomycetes (30.2%), Orbiliomycetes (1.0%), Leotiomycetes (0.7%), Saccharomycetes (0.6%), and Agaricomycetes (0.2%), with ∼36% of sequences not identified to the class level (Fig. 4). Through ANCOM-BC analysis, we detected a greater relative abundance of phylum Basidiomycota in P. × jackii samples (which were absent in P. deltoides), as well as the classes Agaricomycetes, Eurotiomycetes, and Orbiliomycetes, while the class Sordariomycetes was more abundant in P. deltoides samples (Fig. 4 & Table S3). We detected many differentially abundant fungal ASVs between host taxa; the most abundant of which belonged to genus Neoconiothyrium (Dothideomycetes), which was detected in greater abundances in P. deltoides, while ASVs belonging to genus Xenocylindrosporium (Eurotiomycetes), Orbilia (Orbiliomycetes), and Knufia (Eurotiomycetes) were more abundant in P.× jackii. Additional differentially abundant genera and ASVs were detected; however, their relative abundances were low (Table 3, Tables S3 & S4).

The most dominant bacteria of endophytic twig communities belonged to phylum Actinobacteriota (47.6%) and Proteobacteria (44.9%), with smaller abundances of Bacteroidota (5.3%). ASVs belonging to Abditibacteriota, Acidobacteriota, Armatimonadota, Bdellovibrionota, Deinococcota, Firmicutes, Fusobacteriota, Gemmatimonadota, and Myxococcota were also detected at relative abundances of less than 1%. Through ANCOM-BC analysis, we detected a greater relative abundance of phylum Gemmatimonadota in P. deltoides samples, while Acidobacteriota were in greater abundances in P. × jackii samples. The classes Acidimicrobiia, Longimicrobia, and Polyangia were more abundant within P. deltoides, while Acidobacteriae and Armatimonadia were more abundant in P. × jackii (Fig. 4 & Table S5). We detected many additional differentially abundant ASVs and taxa at lower taxonomic levels (Table 4, Tables S5 & S6).