High-Throughput DNA Metabarcoding of Stem Sections from Trees with Cavities Describes Fungal Communities Associated with Variable Wood Decay, Position on Stem and Tree Species

: The presence of unculturable fungi, variability in timing and frequency of fungal fruiting, hyper-rich fungal communities, and genetic and environmental variability explains the difﬁculty in adopting ideal sampling schemes and fungal identiﬁcation approaches in studies of fungal communities in wood at variable stages of decay. Here, we use intensive within-tree sampling of ﬁve standing trees with cavities paired with high-throughput DNA metabarcoding, to study fungal communities in decayed and healthy wood of trees from two Populus species in British Columbia, Canada. The ampliﬁcation of over 3000 fungal DNA sequence variants shows the presence of a hyper-rich wood fungal community that not only varied depending on PCR primers, tree species, tree stem portion and wood decay stage. but also included a large number of taxa unassignable to any known sub-kingdom taxonomic order based on published DNA sequences. Our data show that the use of two different primer sets greatly increases the power of the metabarcoding analysis. By testing three alternative models of fungal community composition, we identify the model that best explains fungal community by considering the position on the stem and distance from the cavity. We suggest this model may be used to design optimal sampling schemes to describe fungal communities in trees experiencing discrete decay pockets or cavities.


Introduction
Studies of fungal communities associated with wood in trees have been undertaken for decades [1,2], and they continue to be undertaken for purposes often as disparate as understanding the role of fungi in the nutrient cycle and C sequestration [3][4][5], fungal endophytism in trees [6], and investigating the associations between wood decay and arthropods [7] or between cavity-nesting birds and fungi [8,9].
What is known about natural wood communities in forests comes from threedimensional detailed studies of fungi in trees or logs to identify fungal species [10][11][12][13][14], from the identification of fungi from the wood of trees obtained with a coarse, less-detailed sampling [15][16][17], and from studies of logs on the forest floor [18][19][20][21]. Notwithstanding the sampling approach and its intensity, the identification of fungal species in wood has been based on a variety of approaches, including sporocarp surveys [22][23][24][25]; fungal culturing, alone or coupled with DNA sequencing of cultures [26][27][28][29][30]; direct amplification of environmental fungal DNA from wood [6]; or cloning followed by traditional Sanger sequencing [8]. The fungal identification methodologies mentioned above are all valid, but they are limited by that fact that the majority of fungi are known to produce visible fruiting bodies only rarely [31], may produce ephemeral fruiting bodies depending on climatic cues [32], or, even if theoretically culturable, they may not be always easily cultured in can improve the explanatory power of permutational analysis of variance (PERMANOVA) when describing the composition of wood fungal communities. The results of our study can describe in detail fungal communities in wood and assist in the design of optimal sampling and analytical strategies in the study of wood decay from standing trees.

Tree Selection and Sampling Approach
Five cavity-harboring live trees uninhabited by wildlife and belonging to the genus Populus (three P. tremuloides and two P. balsamifera) in sub-boreal mainland British Columbia were selected for intensive destructive sampling ( Figure 1). We overcame the difficulty of selecting trees with a truly comparable level of decay by selecting trees with cavities that were judged by field biologists as having the specific characteristics required for being "precursors" to dens for fishers (Pekania pennanti), but not yet large enough to be inhabited [49], and without any sign of stable wildlife inhabitation. This approach ensured that our selected trees were rather uniform in size with similarly sized cavities in comparable states of cavity formation.

Metabarcoding Analysis
DNA was extracted with a Qiagen QiaAMP Fast DNA Stool Minikit following the manufacturer's protocol and using 0.125 g of wood tissue, pulverized with an MP Biomedicals FastPrep 24 Instrument Homogenizer. We employed a two-step phased amplicon sequencing for our polymerase chain reaction (PCR) amplifications. The first step consisted of the amplification of two subsamples from each extraction using the ITS1-F and ITS-2 primer combination to target the ITS-1 region [50] and the 5.8S-Fun and ITS4-Fun primer combination to target the ITS-2 region [51]. These primers consisted of the published primer sequences plus one, two, three, four, or five random nucleotides to act as spacers; and stubs, which complimented our Illumina adapters. This phased approach has been shown to solve the problem of low diversity in the earliest cycles of Illumina sequencing (the variable length of the primers creates a sequencing frame shift that creates base diversity), improve base read quality, increase raw sequence output by approximately 15% (phasing eliminates the need to add PhiX to the sequencing run), and reduce sequencing errors [52]. The PCR cycling profile used was as follows: a 2 min hot-start at 95 °C followed by 35 cycles of 30 s at 95 °C for denaturing, 30 s at 50 °C (for ITS-1) or 55 °C (for ITS-2) for primer annealing, and 1 min at 72 °C for primer extension, followed by 5 min at 72 °C. The PCR products were verified using gel electrophoresis on 1% agarose gels and visualized using UV light. The second step consisted of attaching Illumina adaptors and barcodes to the PCR products from the first step. Dual-indexed adapters were provided by the Vincent J. Coates Genomics Sequencing Laboratory, California Institute for Quantitative Biosciences (QB3), University of California, Berkeley, CA, from their custom set of 960 compatible dual-index pairs (https://qb3.berkeley.edu/facility/genomics/ last accessed on 13 April 2023). The PCR protocol used included a 2 min hot-start at 95 °C followed by 15 cycles of 30 s at 95 °C for denaturing, 30 s at 58 °C for primer annealing, and 50 s at 72 °C for primer extension, followed by 10 min at 72 °C. AMPure XP beads at a ratio of 1.8 were used for clean-up after each PCR amplification. PCR products were quantified using a NanoDrop 2000 spectrophotometer and normalized through dilution in PCR water. Fragment analysis of normalized pools was completed using an Agilent Bioanalyzer, and, finally, all samples were pooled and were sequenced by QB3 using the Illumina MiSeq platform at 2 × 300 cycles.
Prior to OTU clustering, all sequences were filtered and truncated so that all reads would have a minimum length of 100 base pairs and all base pairs would have a minimum Phred Quality score or Q-score of 30 (probability of error in the call = 0.001). We used the USEARCH Ultra-fast sequence analysis pipeline [53], and denoising/ZOTU clustering (at We sampled each tree in evenly spaced planes perpendicular to the long axis of the tree bole. Four sampling planes were situated at 15 cm intervals above the cavity and four were situated below the cavity, plus an additional sampling plane included the center of the cavity itself, resulting in a total of nine sampling planes per tree. At each sampling plane, a 1-inch-thick tree section, or "cookie", was removed from the bole of each felled tree. Before cutting each stem section, the bark areas subject to cutting were surface-sterilized by spraying a 10% solution of commercial bleach and wiping off any residual bleach with a clean paper towel. Before and in between cuts, the chainsaw was wiped clean of any wood debris using a metal brush and surface-sterilized using a sprayed Lysol application. As soon as they were cut, wood cookies were placed and sealed in a double plastic bag. Bags were kept refrigerated and shipped to U.C. Berkeley with next day service within 48 h from collection. Within 24 h from the samples' arrival at U.C. Berkeley, wood cookies were placed in a sterile laminar flow hood and at least four samples of drill shavings, representative of the maximum decay level present in the cookie, were collected from each sampling plane and pooled to create one composite sample per cookie. A total of 168 samples were collected and lyophilized and kept at −20 • C until further processing. All tree sections were photographed, and photos were then used to assign decay stage to individual samples. In the event of no differences between the fungal communities of P. tremuloides and P. balsamifera trees, we wanted to be able to discern whether lack of differentiation may have been a result of the limitation of our methods. In order to do this, two trees from a location clearly distinct from the Populus study site ( Figure 1) and belonging to the phylogenetically distant tree species Cupressus nootkatensis (yellow cedar) were included as outgroups. These trees had large basal cavities and were non-destructively sampled by collecting and pooling three drill shavings of the cavities at breast height per tree.

Metabarcoding Analysis
DNA was extracted with a Qiagen QiaAMP Fast DNA Stool Minikit following the manufacturer's protocol and using 0.125 g of wood tissue, pulverized with an MP Biomedicals FastPrep 24 Instrument Homogenizer. We employed a two-step phased amplicon sequencing for our polymerase chain reaction (PCR) amplifications. The first step consisted of the amplification of two subsamples from each extraction using the ITS1-F and ITS-2 primer combination to target the ITS-1 region [50] and the 5.8S-Fun and ITS4-Fun primer combination to target the ITS-2 region [51]. These primers consisted of the published primer sequences plus one, two, three, four, or five random nucleotides to act as spacers; and stubs, which complimented our Illumina adapters. This phased approach has been shown to solve the problem of low diversity in the earliest cycles of Illumina sequencing (the variable length of the primers creates a sequencing frame shift that creates base diversity), improve base read quality, increase raw sequence output by approximately 15% (phasing eliminates the need to add PhiX to the sequencing run), and reduce sequencing errors [52]. The PCR cycling profile used was as follows: a 2 min hot-start at 95 • C followed by 35 cycles of 30 s at 95 • C for denaturing, 30 s at 50 • C (for ITS-1) or 55 • C (for ITS-2) for primer annealing, and 1 min at 72 • C for primer extension, followed by 5 min at 72 • C. The PCR products were verified using gel electrophoresis on 1% agarose gels and visualized using UV light. The second step consisted of attaching Illumina adaptors and barcodes to the PCR products from the first step. Dual-indexed adapters were provided by the Vincent J. Coates Genomics Sequencing Laboratory, California Institute for Quantitative Biosciences (QB3), University of California, Berkeley, CA, from their custom set of 960 compatible dual-index pairs (https://qb3.berkeley.edu/facility/genomics/ last accessed on 13 April 2023). The PCR protocol used included a 2 min hot-start at 95 • C followed by 15 cycles of 30 s at 95 • C for denaturing, 30 s at 58 • C for primer annealing, and 50 s at 72 • C for primer extension, followed by 10 min at 72 • C. AMPure XP beads at a ratio of 1.8 were used for clean-up after each PCR amplification. PCR products were quantified using a Nan-oDrop 2000 spectrophotometer and normalized through dilution in PCR water. Fragment analysis of normalized pools was completed using an Agilent Bioanalyzer, and, finally, all samples were pooled and were sequenced by QB3 using the Illumina MiSeq platform at 2 × 300 cycles.
Prior to OTU clustering, all sequences were filtered and truncated so that all reads would have a minimum length of 100 base pairs and all base pairs would have a minimum Phred Quality score or Q-score of 30 (probability of error in the call = 0.001). We used the USEARCH Ultra-fast sequence analysis pipeline [53], and denoising/ZOTU clustering (at 100% read similarity with a minimum of eight reads) was completed using the UNOISE algorithm with singleton reads of at least 97% similarity then added to ZOTUs [54,55]. SINTAX searches of Warcup Fungal ITS Training Set Vol. 2, a database of the full ITS region [56], were used to assign OTUs to taxa. A maximum of 5% error rate was applied to taxonomic assignments at each taxonomic level.
Statistical analyses were carried out using the R Statistical Computing package Bio-diversityR [57]. Pairwise analysis of variance was used to explore differences in fungal community based on tree genera, tree species, tree ID, and choice of primer set used to identify community members. Samples from different stem sections within the same tree were analyzed as nested within the tree they came from. An incremental R 2 approach was combined with permutational analysis of variance (PERMANOVA) to further explore independent and interactive effects of each variable, providing a base model, and to explore the effectiveness of our de novo models of community succession.
As we did not know the time of stem injury, and each tree cavity was slightly different in size (and, in all probability, start date), we used the inner upper and lower limits of the cavity as our lower and upper boundaries to take the first section above and below the cavity, thus homogenizing decay level in our samples. In other words, our first sections above and below the cavity may have been at different distances from the center of the cavity, but they all were taken at the upper and lower limits of the hollow inside the tree.
We then incorporated three competing models about fungal community structure in the tree stem before decay to determine which one of the three resulted in the greatest increase in explanatory power, i.e., produced the largest R 2 increase. The three models are here described: (a) tree stem fungal community is uniform throughout the tree prior to cavity formation; (b) tree stem fungal community in the lower section of the stem differs from the one in the upper section, but it is homogeneous within each section; (c) tree stem fungal community is highly heterogeneous, with different communities, each inhabiting a different section of the tree. These three starting points or scenarios will be further referred to as models A, B and C, respectively. In any of these scenarios (Figure 2), each community from different tree samples along the stem was assigned to a different categorical variable, depending on the model. Wood from the cavity center was always assigned to its own category in all three models, but assignment of communities from stem sections above and below the cavity varied depending on the model. For instance, in model A, communities from multiple stem sections of undecayed wood in the lower and upper part of the stem would be placed in the same category, while in model B, communities from above and below the cavity would be placed in two different categories depending on whether they were above or below the cavity. In model C, the community from each stem section was placed in its own different category, irrespective of whether it was above or below the cavity. Once that was performed, we asked PERMANOVA to determine how much of the variance was explained by each categorical variable from the different models and used an incremental R 2 approach to select the best model, i.e., the model with the best explanation of the variance.  Putative de novo models were created to consider three starting points of wood-inhabiting fungal communities within tree stems prior and during cavity formation: (A) that tree stem fungal community is uniform throughout the tree prior to cavity formation; (B) that prior to cavity formation, tree stem fungal community in the lower section of the stem differs from the one in the upper section, but is homogeneous within each section; (C) that tree stem fungal community is highly heterogeneous, with different communities, each inhabiting a different section of the tree. Different colors indicate theoretical different fungal communities and labels of individual samples are not supposed to be readable.
Additionally, we examined the photographs of the tree cookies that provided our wood samples independently from our de novo succession models. Each cookie was assigned a decay class based on Maser et al.'s five-class system of decay in fallen logs [47]. Due to our low sample size, we modified this five-class system into a three-class system to ensure adequate representation in each class ( Figure 3). Fungal communities for P. tremuloides and P. balsamifera, individually, as well as combined, were then compared using PERMANOVA analyses of models that included only the decay classes and of models that included decay class, tree species, and individual tree identification. All community differences were further explored using pairwise comparisons. Indicator Species/OTUs for all factors that had significant effects on fungal community were examined using the Figure 2. Putative de novo models were created to consider three starting points of wood-inhabiting fungal communities within tree stems prior and during cavity formation: (A) that tree stem fungal community is uniform throughout the tree prior to cavity formation; (B) that prior to cavity formation, tree stem fungal community in the lower section of the stem differs from the one in the upper section, but is homogeneous within each section; (C) that tree stem fungal community is highly heterogeneous, with different communities, each inhabiting a different section of the tree. Different colors indicate theoretical different fungal communities and labels of individual samples are not supposed to be readable.
Additionally, we examined the photographs of the tree cookies that provided our wood samples independently from our de novo succession models. Each cookie was assigned a decay class based on Maser et al.'s five-class system of decay in fallen logs [47]. Due to our low sample size, we modified this five-class system into a three-class system to ensure adequate representation in each class ( Figure 3). Fungal communities for P. tremuloides and P. balsamifera, individually, as well as combined, were then compared using PERMANOVA analyses of models that included only the decay classes and of models that included decay class, tree species, and individual tree identification. All community differences were further explored using pairwise comparisons. Indicator Species/OTUs for all factors that had significant effects on fungal community were examined using the Dufrêne and Legendre Indicator Species tests [58], and differences were visualized using non-metric multidimensional scaling (NMDS) with a Bray-Curtis dissimilarity metric. Dufrêne and Legendre Indicator Species tests [58], and differences were visualized using non-metric multidimensional scaling (NMDS) with a Bray-Curtis dissimilarity metric.

Results
Our high-throughput fungal ITS amplicon sequencing project yielded over 25.87 million reads, from which 3467 fungal OTUs were identified, combining the results of both primer sets and assigning each OTU to the lowest taxonomic rank that could be achieved while maintaining an error rate of <5% (Tables 1 and S1). Not surprisingly, a large number of OTUs (Table 2) could not be confidently matched to sequences of known species (2660 had no match; 410 had a low-confidence match), genera (1595 had no match; 867 had a low-confidence match), or family (1213 had no match; 930 had a low-confidence match). However, the lack of matches was more surprising at the higher taxonomic levels. In fact, 751 OTUs could not be assigned to any known order (1160 had low-confidence assignments), 525 OTUs could not be assigned to any known class (780 had low-confidence

Results
Our high-throughput fungal ITS amplicon sequencing project yielded over 25.87 million reads, from which 3467 fungal OTUs were identified, combining the results of both primer sets and assigning each OTU to the lowest taxonomic rank that could be achieved while maintaining an error rate of <5% (Tables 1 and S1). Not surprisingly, a large number of OTUs (Table 2) could not be confidently matched to sequences of known species (2660 had no match; 410 had a low-confidence match), genera (1595 had no match; 867 had a low-confidence match), or family (1213 had no match; 930 had a low-confidence match).
However, the lack of matches was more surprising at the higher taxonomic levels. In fact, 751 OTUs could not be assigned to any known order (1160 had low-confidence assignments), 525 OTUs could not be assigned to any known class (780 had low-confidence assignments), and 295 OTUs could not even be assigned to a known phylum (338 had low-confidence assignments), even though they were confidently assigned to the Domain Kingdom Fungi. Table 1. Distribution of distinct OTUs by phylum (unknown taxa in parenthesis) of caulosphere fungi found in the seven sampled trees from the genera Cupressus and Populus. The upper portion of the table was generated using all OTUs produced by the dataset, while the lower portion was generated using only OTUs confidently (p < 0.05) identified to the species level.  The top 20 families of Ascomycota and Basidiomycota, based on the number of OTUs that could be confidently assigned to them, is shown in Table 3, while Table 4 lists the top 20 OTUs in these phyla by read abundance. It should be noted that our pipeline designated OTUs as Sequence Variants (SVs) based on a precise 100% similarity of sequence (followed by clustering of singleton reads within a 3% "cloud" of similarity), hence there were OTUs that, being highly similar but not identical, were treated separately even though they blasted to the same species (i.e., five or more distinct OTUs were confidently assigned to each of six known species within the database; four distinct OTUs were confidently assigned to each of eighteen known species; three distinct OTUs were assigned to each of seventeen species; and two distinct OTUs were assigned to each of sixty-nine known species within the database). Fungal taxa identified as indicator species (see below) of primer pair, tree genus, tree species, and decay class are provided in Table S2.  Accumulation curves were created based on these four factors of interest: primer pair, tree genus, tree species, and individual tree ( Figure 4). These curves display the total number of species observed during data collection as additional samples were added to the pool of all previously analyzed samples. These curves allow for a comparison of the alpha-diversity of wood-inhabiting fungal communities associated with different groupings. When slopes of these curves flatten, as if they had reached an asymptote, this result suggests the scale of sampling was adequate to properly describe the fungal communities in the sample. The curves in Figure 4 indicated that the 5.8S-Fun&ITS4-Fun PCR primer combination consistently amplified 30%-50% more fungal OTUs than the ITS1F&2 combination did, likely due to its affinity for the more species-rich phylum, the Ascomycota. Of the 397 OTUs confidently assigned to species, 99 were only amplified by one primer set or the other. Pairwise comparisons of fungal community amplified by each primer set show significant differences (p-value = 0.027), and the NMDS analysis clearly indicates OTUs detected by each primer pair are largely not overlapping ( Figure 5A).  Populus had three times as many fungi per sample than Cupressus; however, this was an expected result given the difference in sample size, and OTU accumulation curves (Figure 4) reveal that this is likely due to the small number of Cupressus samples. Notwithstanding the difference in the total number of OTUs and the analytical limitations of such difference, the NMDS analysis revealed that the fungal community of Cupressus nootkatensis was clearly distinct from that of Populus as a whole ( Figure 5B). Additionally, even if specific findings regarding the community composition of wood fungi from C. nootkatensis are not discussed in detail here, because of the small sample size, an initial PERMANOVA , not only do we show accumulation curves for the data subset that was identified at the species level, but we can also draw confidence intervals, given that the two curves are clearly distinct. For all other curves in panels (B,C), curves are overlapping and so are their confidence intervals.
Populus had three times as many fungi per sample than Cupressus; however, this was an expected result given the difference in sample size, and OTU accumulation curves (Figure 4) reveal that this is likely due to the small number of Cupressus samples. Notwithstanding the difference in the total number of OTUs and the analytical limitations of such difference, the NMDS analysis revealed that the fungal community of Cupressus nootkatensis was clearly distinct from that of Populus as a whole ( Figure 5B). Additionally, even if specific findings regarding the community composition of wood fungi from C. nootkatensis are not discussed in detail here, because of the small sample size, an initial PERMANOVA analysis identified genus as a highly significant variable (p = 0.0009). Populus tremuloides had approximately 20% higher species richness than P. balsamifera, and the NMDS analysis partitioned wood fungal communities from individual trees in two distinct and non-overlapping clusters based on trees species ( Figure 5C). Based on the NMDS analysis ( Figure 5D), individual trees had wood fungal communities that were distinct, even if partially overlapping. Finally, it is interesting to note that even in individual trees with as many as nine pooled samples per tree, rarefaction curves were close to, but had not yet reached, an asymptote ( Figure 5D).
Pairwise analysis of variance of the fungal communities using all distinct OTUs within the genus Populus found significant effects on observed fungal communities at the levels of tree species and individual tree. Our PERMANOVA model (Table 5) found significant independent effects of tree species, tree, and primer pair (which also had interactive effects with each of the other two variables, resulting in primer selection alone accounting for ~21% of the total variation in fungal communities of our three tree species), giving our base model a total R 2 of 0.358. Populus tremuloides had approximately 20% higher species richness than P. balsamifera, and the NMDS analysis partitioned wood fungal communities from individual trees in two distinct and non-overlapping clusters based on trees species ( Figure 5C). Based on the NMDS analysis ( Figure 5D), individual trees had wood fungal communities that were distinct, even if partially overlapping. Finally, it is interesting to note that even in individual trees with as many as nine pooled samples per tree, rarefaction curves were close to, but had not yet reached, an asymptote ( Figure 5D).
Pairwise analysis of variance of the fungal communities using all distinct OTUs within the genus Populus found significant effects on observed fungal communities at the levels of tree species and individual tree. Our PERMANOVA model (Table 5) found significant independent effects of tree species, tree, and primer pair (which also had interactive effects with each of the other two variables, resulting in primer selection alone accounting for 21% of the total variation in fungal communities of our three tree species), giving our base model a total R 2 of 0.358. A more focused exploration on the genus Populus that included the location of the sample on each tree with respect to the cavity, using PERMANOVA and an incremental R 2 approach, indicated that, although pairwise comparisons did not show significant differences between any of the sampling points along the tree trunk, the absolute distance of the sample from the center of the cavity had independent effects as well as significant interactions with tree species and tree, explaining an additional 15% of fungal community variation (Table 6). Table 6. Comparison of models of community response using additive and interaction effects of (a) tree species (Species); tree ID (Tree); distance of sample from cavity (Distance); and location of sample in relation to cavity (above vs. below or Ab vs. Be), or (b) tree species (Species), tree ID (Tree), best de novo model of decay spread (Model B), and location of sample in relation to cavity (above vs. below or Ab vs. Be).
(a) Y = Species * Tree * Distance * Ab vs. Be When our de novo models of succession (i.e., models A, B, and C) were added to the above base model, the primary models A and B had significant independent additive effects and increased R 2 by~0.025 and~0.045, respectively; additionally, both also had interactive effects that increased R 2 by~0.072 and~0.197, respectively. The primary model C had no additive affects, and interactive effects could not be determined as the similarity between the community structure factor and the distance from the cavity resulted in 0 as the sum of squares and an infinitely negative mean squares value. The diverse, abundant community structure presumed in model C provided no improvement when compared to the basic additive model (tree species, tree ID, and distance from cavity). In the end, the B iteration of our models of succession held the most explanatory power. Table 6 shows that in comparison with a model that considered only the distance of the sample from the cavity and the sample's placement above or below the cavity, our model explained 5.8% more variability in fungal communities, with a total R 2 of 0.815.

Sum of Sqrs
The Indicator Species tests [58] were also conducted to explore trends in fungi that were amplified by ITS1F&2 and 5.8S-Fun&ITS4-Fun. Of the 3467 OTUs, 739 were found to be Indicator Species of a primer set, 358 of ITS1F&2, and 381 of 5.8S-Fun&ITS4-Fun (Table 7). When tree species were compared (excluding the outgroup C. nootkatenis) combining both primer sets, there were 44 P. tremuloides indicators, while P. balsamifera had 210 (Table S2). Pairwise comparisons of fungal community based on sample decay class found significant differences between all three decay classes in our modified decay classification system for P. tremuloides, P. balsamifera and for both species combined ( Table 8). The PERMANOVA analysis of a model including tree species, individual tree, and decay class confirmed the effect of decay class on the fungal community (Table 9), and this was visualized using NMDS analysis ( Figure 6).   (panels A,B) and then with data from both tree species combined (panel C).

A Hyper-Rich Fungal Community Is Found in the Wood of Two Populus Species
The number of reads in this experiment was high, and so was the number of sequence  panels (A,B)) and then with data from both tree species combined (panel (C)). Table 9. Models of community response using additive and interaction effects of tree ID (Tree), and wood decay class (Decay Class). Models are presented for each Populus species independently (top two panels) and for the two Populus species combined (bottom panel). Overall, and in ranking order in terms of number of reads, the top five Ascomycota families identified in wood by this study were the Xylariaceae, Saccharomycetaceae, Herpotrichiellaceae, Helotiaceae, and Diatrypaceae. The combined number of reads of these five families (~2.6 million) is comparable to the number of reads from the top five ranking Basdiomycota families (~2.1 million), indicating both groups of fungi are well represented in Populus wood. However, the overall number of reads of ascomycetes (~5.8 million) is basically double the number of reads of Basidiomycota (~2.9 million) obtained from the same substrate. Overall, and in ranking order in terms of reads, the top five families of Basidiomycota were the Hymenochaetaceae, the Pleurotaceae, the Hydnaceae, the Sporidiobolaceae, and the Auriscalpiaceae. Although the Tremellaceae ranked first in number of OTUs, they were not included in the top five families in terms of number of reads.

Sum of Squares
When comparing fungal communities from Populus samples in the three distinct classes of decay based on the photographs (1 or healthy, 3 or intermediate, and 5 or very advanced; (Figure 3), communities were significantly different among the three classes. This was true when comparing samples within species (P. tremuloides or P. balsamifera) or when comparing the two species together (Tables 8 and 9). P. tremuloides had the largest sample size and thus fungal communities associated with different stages of decay in this species were analyzed more in-depth by looking at the top-ranking 10 OTUs in terms of the number of reads in each wood decay stage. Table 10 lists the top ten fungal OTUs (by number of reads) for each of the decay classes analyzed above for P. tremuloides. Although the same OTUs may be present in samples with different decay levels, the three different wood decay stages showed almost no overlap in terms of the 10 OTUs with the largest number of reads, and minimal overlap in terms of fungal families represented. In decay class 1, i.e., healthy or minimally stained wood, the largest number of reads belonged to OTUs that, although fungi, could not even be assigned to any known fungal class, order, or family. A comparable high number of OTUs in decay stage 1 belonged to the Malasezziaceae (Malasezzia spp.) and Saccharomycetaceae (unknown spp.), two families of ascomycetes both known to include a vast number of yeast species found in the most diverse substrates in the world, including wood [59,60]. Finally, the Leotiomycetes and the Agaricomycetes were also very well represented in decay class 1 samples. Although the lifestyles of species in these two fungal classes can be extremely varied, they both include many species that can be endophytic in at least one stage of their life.
In decay stage 3 samples, i.e., samples with intermediate decay, the dominant fungal community structure changes dramatically, and the pathogenic genus Phellinus (Hymemochetaceae, Agaricomycetes) was clearly dominant. Once again, yeasts were important and ranked second in terms of reads. However, in decay stage 3, dominant yeast OTUs were different from those detected from in the stage 1 decay samples. Here, Candida and Pichia spp. (Saccharomycetacae) dominated. Finally, a Neobulgaria sp. (Helotiaceae) was also abundant in this decay stage.
In decay stage 5, i.e., samples with the most advanced decay, once again, we observed a shift in the composition of the dominant fungal community, with different dominant OTUs ranking as the top ten most abundant species. Here, OTUs belonging to the Xylariaceae (Xylaria spp., Hypoxylon macrocarpum, unknown species) and the Hydnaceae (e.g., Sistotrema raduloides) were dominant. Among the most dominant OTUs, we found the yeast Capronia pilosella (Herpotrichiellaceae).
The NMDS ordination performed to visualize differences among the three decay stages communities in P. tremuloides ( Figure 6) shows a clear separation of samples in decay stage 1 and 5, while samples in decay stage 3 overlapped with the other two stages. Samples from decay stage 1 were included within decay stage 3 samples, and two decay stage 5 samples appeared to be disjunct from all other decay stage 5 samples.

A Hyper-Rich Fungal Community Is Found in the Wood of Two Populus Species
The number of reads in this experiment was high, and so was the number of sequence based fungal OTUs, well exceeding 3000. This high number, generated by combining OTUs obtained with both primer pairs, is obviously an inflated estimate of the total number of species for two reasons. First, OTUs with highly similar but not identical DNA sequences were regarded as different Sequence Variants (SVs) [61]. The presence of highly similar yet different SVs may indicate that classic morphology-based fungal taxonomy still has to catch up with the phylogenetic-based species concept [62], or it may indicate the presence of significant intraspecific genetic variation [63,64]. Second, multiple OTUs generated by different primer pairs and assigned to the same taxonomic level, higher than species (e.g., both sequences were matched to the same family), could, in fact, originate from the same individual, and thus should be counted as the same OTU. The determination of how many of these 3467 OTUs may actually be duplicates is virtually impossible, given the analytical tools employed in this study, with the exception of OTUs that could be assigned to species, in which case the ITS1 and the ITS2 would have been matched to the same complete ITS sequence. Although the number of OTUs may be halved to 1733 if we had complete duplication of results between the two primer sets, complete OTU duplication is highly unlikely given that 28.5% of 347 OTUs that were confidently identified to the species level were unambiguously amplified only by one primer set or by the other. Accordingly, by applying these percentages to the entire OTU dataset, we estimate that the maximum percentage of full duplication may only occur for approximately 72% of the total number of OTUs, resulting in a maximum 36% drop (76/2 = 36) in the total number of OTUs from 3347 to 2142.
Given that the purpose of this study was not to describe in-depth the number of fungal species within a single substrate, but rather to compare fungal communities from different hosts, stem locations, and decay stages, the inclusion of fungal intraspecific genetic variability actually increases the power of the analysis. In fact, any intraspecific variability associated with host, stem location, or decay stage would be lost if similar but non-identical sequences were clumped together in individual OTUs. Of course, using 100% sequence similarity has its own drawbacks, and it implies limited sequencing error. While random errors are not likely to affect the results, systematic reaction-or substrate-based errors could bias the results. In order to increase the confidence of the dataset, we used replicated biological filters and sequences with high Phred Q-scores to ensure that neither chance nor sequencing ambiguities may be responsible for the analytical outcomes. It is also encouraging that a recent study has shown high correspondence between haplotypes generated by HTS and those obtained by Sanger sequencing [64]. Finally, while using 97% similarity among sequences to generate a single OTU would have allowed for a better estimate of the number of species [64,65], that approach has its own limitations. First, the 97% threshold for species identification is not uniform across the fungi [63]; second, pooling multiple sequences in a single OTU inevitably decreases the estimation of the abundance of the species, genus, family, or order that the OTU has been assigned to. If abundance and substrate preference is of interest, then using a 100% similarity approach may be a better alternative, especially when many sequences in the database cannot be identified with confidence at the species level [56]. This approach is echoed by many ecologists and epidemiologists, and a call has been made to use sequence variants or SVs rather than a priori OTUs [66]. Our approach effectively is based on an analysis of SVs.

Wood Fungal Communities Vary Depending on PCR Primers Used in High-Throughput DNA Metabarcoding
Based on the PERMANOVA, NMDS, and Indicator Species analyses ( Figure 5, Table 5), our data show that primer set selection can have strong effects on the description of fungal communities. This has been shown when primers amplify different loci [17], but also when they amplify the same locus [67]. This result indicates that using multiple primer pairs may provide a more complete and less biased description of fungal communities in wood than using a single primer pair [68]. The debate has been intense about the suitability of the ITS region as the genetic barcode for fungi; however, the ability of ITS to discriminate a large number of fungal species [69], including thousands of wood-inhabiting agarics [65], and the large number of fungal taxa for which the ITS has already been sequenced [63], all but guarantee the likelihood that it will remain a cornerstone of fungal ecology. Tedersoo et al. [68] echo this sentiment; however, until high-throughput amplicon sequencing technology improves to the point of capturing the entire ITS region-researchers are faced with having to choose ITS1 vs. ITS2, each of them with its own pros and cons [51,65,70]. We suggest that one way to minimize primer bias and maximize captured diversity, while working around high-throughput sequencing limitations on read length, is to employ multiple primer sets which target ITS1 and ITS2 independently, as we have done by using both ITS1F&2 and 5.8S-Fun&ITS4-Fun.

Wood Fungal Communities Vary Depending on Tree Species
Tree genus and species were found to have a significant effect on fungal community ( Figure 5; Tables 5 and 6) and are likely to be drivers of fungal biodiversity on the forest floor. Accordingly, they should be taken into consideration in any forest management plan that aims to conserve forest biodiversity. Each tree species had a sizeable number of Indicator Species [58], showing that some wood-inhabiting fungi are host-specific, even when comparing closely related congeneric Populus host species. While differences in woodinhabiting fungal communities have been reported for distantly related tree species [5], and host specificity of many wood-inhabiting fungi is well known [71], this is one of the first reports on a clear distinction in the entire fungal community composition of wood-inhabiting fungi when comparing closely related congeneric tree species.

The Composition of Wood Fungal Communities Varies Depending on Wood Decay Stage
The PERMANOVA, NMDS, and Indicator Species analyses all showed that fungal communities associated with different wood decay stages were different ( Figure 6; Tables 8, 9 and S2). Likewise, dominant fungi differed when comparing various wood decay stages. In samples at decay class 1 (the class described as having little to no visible decay), large read numbers were found for unknown taxa in the Agaricomycetes and Leotiomycetes. The lack of precise species-level identification for many sequences in this case indicates there are taxa within known fungal groups that either have still to be identified or that have yet to be sequenced. Both fungal classes above, however, are known to include wood endophytes, which is consistent with our detection of these fungi in healthy wood. Malasezziaceae were the dominant yeasts in healthy wood. Their dominant presence in wood is interesting, given that Malasezzia spp. are increasingly being found in a variety of terrestrial, marine, animal, and plant substrates and have been identified as having an unexpectedly rich array of catabolic and metabolic genes [60], a result that could explain their ability to live in healthy wood.
Phelllinus spp. and Phellinus tremulae, in particular, were dominant in the intermediate decay class 3; these are fungi that can persist as endophytes before shifting to a pathogenic lifestyle. It makes sense that decay would be initiated by this genus and by P. tremulae, in particular, given that this is a well-known pathogen of Populus tremuloides, and its fruitbodies were observed in the study sites, although not on the trees we sampled. Phellinus spp. were dominant in decay class 3 wood, but this does not exclude that they may have been present in samples displaying lesser decay. The dominance of a Neobulgaria sp. in decay stage 3 is also interesting, because species within the genus are known to cause wood discoloration in healthy trees and to be frequently found in decayed wood [72]. Hence, based on the literature, the presence of species belonging to this genus might be expected in early to intermediate stages of decay, although their dominance is probably a novel and relatively unexplored result.
It is interesting to note that yeasts in the Saccharomycetaceae were dominant both in decay class 1 and 3. Although the OTUs belonging to this family were different in different decay classes: this finding suggests that these yeasts are capable of an endophytic lifestyle in healthy wood and may later contribute to decay wood. Many Saccharomycetaceae have also been associated with the formation of soft rot, a secondary decay process that can persist even when conditions become unfavorable to other wood decay fungi [73,74]. Given the dominance of these yeasts in wood and given the shift in yeast OTUs that were dominant in different wood decay classes, we believe the role of yeasts as endophytes in healthy trees, as well as their role in the wood decay process, needs to be studied in-depth.
The dominance of Hypoxylon macrocarpum and Sistotrema raduloides in decay class 5 wood samples from P. tremuloides is also of great interest. Both species are known to be weak pathogens and excellent saprobes, hence, their presence in the later stages of decay is to be expected, although we believe this to be one of the first studies reporting their dominance in late-stage decay of Populus spp. Both species have the ability to initiate wood decay in live trees when tree health is compromised by biotic (pathogens, insects) and/or abiotic (drought, flooding, cold spells, mechanical damage, soil compaction, etc.) factors. Hypoxylon macrocarpum is a cosmopolitan sapwood rotter, commonly found in plants belonging to the Salicaceae family, as well as in other plant families usually adapted to riparian or mesic ecosystems. It is able to degrade both lignin and cellulose, further assisting in the breakdown of wood elements previously initiated by pathogens and wood primary decay agents. Sistotrema raduloides, instead, is possibly regarded as a hypersaprobe, given that it is unclear whether it possesses the enzymatic ability to degrade wood [75]. The black Herpotrichiellaceae are the dominant yeasts in decay class 5 samples. We note that Capronia pillosula, the dominant yeast species in class 5, is unable to degrade plantderived polymers and chitin [76]; hence, similar to S. raduloides above, it can be defined as a hypersaprobe inhabiting wood previously degraded by primary saprobic organisms.

Fungal Communities Change Progressively as Wood Transitions from Healthy to Decayed, with Some Fungi Being Present in Different Wood Decay Stages
Although we have clearly shown here that (a) entire fungal communities are significantly different when comparing different decay stages and (b) dominant taxa are different in wood characterized by different decay stages, we should also note that decay represents a continuum, with progressive transitions from one stage to the next. Consequently, a further expectation was that some overlap in fungal species may be present when comparing P. tremuloides wood in the three distinct decay stages. The NMDS ordination performed on fungal community data from P. tremuloides samples assigned to three distinct decay stages ( Figure 6) shows such overlap by placing wood samples at decay stage 3 as overlapping both samples belonging to decay class 1 and decay class 5. Additionally, it appears that the community in decay class 1 (undecayed wood) may be mostly included in the cluster of samples with intermediate decay class 3, while the fungal community identified in decay class 3 comprises a significant number of taxa that are well outside the cluster of samples from decay class 1. This result suggests that a large portion of fungi present in wood with incipient to intermediate decay were already present in healthy wood, further reinforcing the concept of a continuum in which endophytic fungi experience a series of lifestyle shifts (e.g., sometimes innocuous, sometimes beneficial, sometimes pathogenic, and sometimes saprobic) throughout their relationship with their hosts in response to changes in their environment.

Decay and Location on Stem Interact to Drive Composition of Fungal Communities in the Wood of Tree Stems
Even if this study is preliminary, we have shown here that investigating the composition of fungal communities in trees displaying cavities allows to differentiate between fungal communities present in decayed and healthy portions of the trunk. Moreover, the high explanatory power of our de novo PERMANOVA model (including tree species, tree ID, a group factor from our theoretical model of succession B, and placement of sample above or below the cavity) outperformed a model that used only tree species, tree ID, distance from cavity, and placement of sample above or below the cavity. However, when the three-class decay assessment was added to our complete de novo model, it did not increase the overall explanatory power of the model, indicating that our de novo model of decay community succession was able to either fully explain all the variation caused by differences in decay class or to explain the variation in the healthy upper and lower parts of the stem that had not yet experienced decay. In other words, the fungal community in the stem area around the cavity is better described when considering distance from the cavity, independent of whether the stem section analyzed is above or below the cavity itself. This suggests decay has a strong effect on the composition of fungal communities, an effect that is dominant on the effect of stem location (upper vs. lower stem) on fungal composition. However, as distance increases from the decay-generated cavity, the effect of decay on fungal community composition wanes and, instead, location on the stem, i.e., upper vs. lower stem, becomes the best predictor of fungal community composition. In summary, the overall tree stem fungal community was best described by looking at the upper stem in an area unaffected by decay, at the lower stem in an area unaffected by decay, and at various distances in the proximity of the cavity generated by decay.

Conclusions
Here, we have proposed that studying wood at various stages of decay in standing live trees with cavities using high throughput metabarcoding may represent a valuable addition to the study of decay in standing trees using culturing or traditional PCR and of decay in dead and downed logs or snags using SST approaches. We acknowledge that the various approaches are complementary and not fully interchangeable, given that the biological and physical processes causing wood decay are different in standing and downed trees. Nonetheless, given that we were able to differentiate fungal communities from: (a) wood of closely related tree species, (b) wood at different decay stages, and (c) wood in different portions of the stem, we believe the approach here described, including the use of a priori models of within-tree community composition, has the potential to facilitate our general understanding of the dynamics of fungal communities inhabiting healthy and decaying wood. Our approach is capable of describing hyper-rich communities, including thousands of taxa, and provides the significant benefits of a homogeneous genetic and spatial background while reducing the potential for an abundance of airborne and soilborne contaminants. The analyses here presented suggest that the use of more than one PCR primer set and a sampling scheme that includes both samples in the higher and lower parts of the healthy portion of a tree stem and samples at variable distances from decay pockets or cavities in the stem has the potential to best describe fungal communities in the caulosphere of trees.