Pathogen and Endophyte Assemblages Co-vary With Beech Bark Disease Progression, Tree Decline, and Regional Climate

Morrison, Eric W.; Kasson, Matt T.; Heath, Jeremy J.; Garnas, Jeff R.

doi:10.3389/ffgc.2021.673099

ORIGINAL RESEARCH article

Front. For. Glob. Change, 24 May 2021

Sec. Pests, Pathogens and Invasions

Volume 4 - 2021 | https://doi.org/10.3389/ffgc.2021.673099

This article is part of the Research Topic Forest Pathology in Changing Climate View all 10 articles

Pathogen and Endophyte Assemblages Co-vary With Beech Bark Disease Progression, Tree Decline, and Regional Climate

$\r\nEric W. Morrison*$ Eric W. Morrison^1*

Matt T. Kasson²

Jeremy J. Heath¹

Jeff R. Garnas¹

¹Department of Natural Resources and the Environment, University of New Hampshire, Durham, NH, United States
²Division of Plant and Soil Sciences, West Virginia University, Morgantown, WV, United States

Plant–pathogen interactions are often considered in a pairwise manner with minimal consideration of the impacts of the broader endophytic community on disease progression and/or outcomes for disease agents and hosts. Community interactions may be especially relevant in the context of disease complexes (i.e., interacting or functionally redundant causal agents) and decline diseases (where saprobes and weak pathogens synergize the effects of primary infections and hasten host mortality). Here we describe the bark endophyte communities associated with a widespread decline disease of American beech, beech bark disease (BBD), caused by an invasive scale insect (Cryptococcus fagisuga) and two fungal pathogens, Neonectria faginata and N. ditissima. We show that the two primary fungal disease agents co-occur more broadly than previously understood (35.5% of infected trees), including within the same 1-cm diameter phloem samples. The two species appear to have contrasting associations with climate and stages of tree decline, wherein N. faginata was associated with warmer and N. ditissima with cooler temperatures. Neonectria ditissima showed a positive association with tree crown dieback – no such association was observed for N. faginata. Further, we identify fungal endophytes that may modulate disease progression as entomopathogens, mycoparasites, saprotrophs, and/or additional pathogens, including Clonostachys rosea and Fusarium babinda. These fungi may alter the trajectory of disease via feedbacks with the primary disease agents or by altering symptom expression or rates of tree decline across the range of BBD.

Introduction

Plant–microbe or plant–insect interactions are often considered in a pairwise manner with minimal consideration of the impacts of the broader community on the nature and outcomes of herbivory or pathogen attack. In some systems, tri- or even multipartite interactions have been shown to be important, particularly where fungus-insect or fungus-insect-mite symbioses are involved (Wingfield et al., 2010, 2016). Rarely, however, are broader communities considered, despite the fact that pathogens and insects attacking forest trees are embedded in variable and/or, in the case of non-native species, novel communities. In the most aggressive and well-studied examples (e.g., Chestnut blight, Dutch elm disease, or the more recent Emerald ash borer) disease agent aggressiveness and ensuing host mortality may be sufficiently rapid such that co-occurrence with other organisms is of minimal relevance to system dynamics. Though, even in these disease systems, colonization with certain endophytic microbes can influence host susceptibility (Feau and Hamelin, 2017) or pathogen aggressiveness (Kolp et al., 2020) in subtle or complex ways. Often, however, interactions with host trees are embedded in a diverse community that varies both spatially and temporally, and disease “complexes” (diseases with multiple causal agents that act in concert to produce symptoms and mediate host decline), are increasingly recognized as important (Desprez-Loustau et al., 2016). Endophytes may influence disease outcome through a variety of mechanisms including competition with pathogens for resources (Oliva et al., 2020) and production of toxins and antifungal compounds (Rodriguez et al., 2009). The role of the community in determining host fate is extremely difficult to ascertain and could proceed via multiple, potentially interacting mechanisms (Table 1).

TABLE 1

Table 1. Summary of ecological guilds in addition to primary pests and pathogens that are potentially involved in tree decline associated with multi-species disease complexes.

Beech bark disease (BBD) in North America is a canker disease of American beech (Fagus grandifolia) caused by an invasive scale insect, Cryptococcus fagisuga Lind., and one of two presumptively native fungal pathogens, Neonectria faginata and N. ditissima. BBD is referred to as a disease complex because it involves both insects and fungi with at least some degree of ecological redundancy with respect to disease etiology and symptom development. The invasive felted beech scale (C. fagisuga) is recognized as the primary initiating agent of BBD symptoms, but only inasmuch as it facilitates the initial infection of beech trees by the fungal BBD pathogens Neonectria faginata and N. ditissima (Ehrlich, 1934) particularly in stands without previous disease exposure (Cale et al., 2012). Both N. faginata and N. ditissima cause similar cankering infections and bole defects resulting in host callus tissue formation (Cotter and Blanchard, 1981). In addition to the primary disease agents, at least two fungal mycoparasites are known from the BBD system – Clonostachys rosea (the anamorph and preferred name of Bionectria ochroleuca; Houston et al., 1987; Schroers et al., 1999; Rossman et al., 2013; Stauder et al., 2020b) and Nematogonum ferrugineum (Houston, 1983). Other fungi are regularly isolated from infected trees (i.e., Fusarium babinda; Stauder et al., 2020b) though their roles are less clear. Further, saprotrophic fungi [e.g., Inonotus glomeratus, Phellinus igniarius (Cale et al., 2015a)] are likely involved in late stages of disease wherein host stem tissue is weakened to the point of mechanical failure (“beech snap,” Houston, 1994), though the presence or importance of wood-rot fungi on disease progression remains unclear.

In addition to questions of the role of the broader community in disease dynamics, the relative frequency and ecological importance of the two primary (Neonectria) pathogens has long been the subject of study and debate. While apparently ecologically similar within the BBD system, the life histories of these fungi differ in important ways. For example, Neonectria ditissima is a generalist that infects many diverse tree hosts including species of birch, maple, walnut, mountain ash, and holly, among others (Castlebury et al., 2006; Stauder et al., 2020b). The species is also an important pathogen of apple (Gómez-Cortecero et al., 2016). The diversity and abundance of alternative hosts could plausibly influence ecological and evolutionary dynamics (Houston, 1994; Kasson and Livingston, 2009), though this question has not been adequately studied to date. In contrast, N. faginata has never been observed outside of the BBD complex in North America (Castlebury et al., 2006) despite recent surveys focused on uncovering possible cryptic native reservoirs for this pathogen (Stauder et al., 2020b).

These fungi exhibit spatiotemporal trends with respect to the timing of site-level infestation with the felted scale. The current range of BBD is defined by the range of the felted beech scale, which, unlike many forest pests, has spread slowly (∼13 km per year; Morin et al., 2007) from the site of initial introduction in Halifax, Nova Scotia in 1890 (Hewitt, 1914). This progression, together with a handful of long-distance dispersal events (i.e., to North Carolina, West Virginia, and Michigan) has resulted in a gradient of duration of infestation ranging from very recent (∼9 years in Wisconsin) to more than eight decades (86 years in Maine) (Houston, 1994; Cale et al., 2017). Surveys of Neonectria species distribution have generally found N. ditissima to be more prevalent in the killing front of the disease (i.e., 10–20 years post scale insect arrival; Cale et al., 2017 and references therein). Neonectria faginata appears to dominate aftermath forests to the degree that researchers have suggested near replacement of N. ditissima with N. faginata as early as 7 years after pathogen attack becomes apparent (Houston, 1994; Cale et al., 2017). However, N. ditissima can maintain a presence in stands dominated by N. faginata, including within the same tree (Kasson and Livingston, 2009). Persistence of N. ditissima in areas that are previously BBD-affected may be attributed to reinfection from reservoirs of this fungus in non-beech hardwood tree hosts (Kasson and Livingston, 2009), and/or the development of secondary killing fronts when climatic conditions allow beech scale to colonize areas where it was previously inhibited (e.g., release from killing winter temperatures by warm periods; Kasson and Livingston, 2012). These species – while morphologically indistinguishable in the field – can be separated using culture morphology and spore size measurements, though the latter process is tedious and is dependent on the presence of sexually produced ascospores. Further, spore size comparisons can only detect co-infection if either both species are simultaneously producing perithecia (sexual spore structures) or multiple isolates are collected, cultured, and induced to mate and produce sexual structures (Cotter and Blanchard, 1981; Stauder et al., 2020a). Partly because of these challenges, it is yet unclear whether trends in species dominance with infection duration are consistent, whether N. ditissima plays an important and/or predictable role in aftermath forests, and how the prevalence and distribution of each species reflects climate, disease stage, and other tree host and environmental conditions.

The objectives of this study are twofold. First, we examine patterns of occurrence (and co-occurrence) of N. faginata and N. ditissima across the current range of BBD and use joint species distribution modeling to evaluate hypothesized biotic and abiotic drivers of the prevalence and relative dominance of these species. Second, we characterize the bark mycobiome of American beech using Illumina-based metabarcoding on bark samples collected from across the range of BBD to ask how disease-associated communities vary geographically and to assess the fidelity and potential role of key species within the BBD system, whether as direct or indirect drivers or as indicators of disease state. Based on previously observed patterns with respect to disease dynamics across the range of BBD, we evaluated the relative contribution of hypothesized drivers of pathogen and associated community distribution. These included the number of years since regional infestation with scale insect (hereafter referred to as duration of BBD infection), disease severity (i.e., tree condition as well as beech scale and Neonectria perithecia density), and climate.

Materials and Methods

Site Description and Sample Collection

Bark disks including phloem tissue were collected from American beech (Fagus grandifolia) from 10 sites across the range of BBD. Sites ranged from northern Maine, western North Carolina, and eastern Wisconsin, representing latitudinal and longitudinal transects across the current range of BBD with a range of infection duration (Table 2). Sampling was performed from December 2017 through January 2019. At each site American beech trees were surveyed for levels of Neonectria perithecium density, beech scale density, crown dieback, and amount of cankering. Scale insect and Neonectria perithecium density were scored on a 0–5 ordinal scale (Houston et al., 2005; Garnas et al., 2011b) and tree condition on a 0–4 scale. Distinct cankering types were pooled and measured as roughly corresponding to 20% bins by bole coverage. Trees were sampled along 100 × 5 m transects in a random direction from starting point up to a maximum of 50 trees or 400 m. All trees were measured along the transects to facilitate estimation of site-level tree size distribution and mean disease severity. Where possible, stratified random sampling was performed for bark plug collection so as to obtain an unbiased sample across a range of tree conditions. Stratified sampling levels were tree size (three levels 1st, 2nd/3rd, and 4th quartile of diameter at breast height [DBH]), four levels of Neonectria perithecium density (0, 1, 2–3, 4–5), and two levels of beech scale density (0–1 vs. 2–5), yielding 24 possible stratification levels. It was not possible to collect all combinations at all sites, but most sites had representative trees in most categories. In particular, in one site (Wisconsin) there were no visible Neonectria perithecia and we instead stratified within DBH and the available levels of beech scale density (2–3, 4–5) with four replicates per stratum (n = 24). Bark plugs were collected using a flame-sterilized 1-cm diameter hollow leather punch and stored on ice in sterile 24-well plates. The punch was flame-sterilized by first wiping it free of surface debris using a clean Kim wipe and then applying a butane torch flame for 20 s. Multiple plugs were taken from a random subset of trees with number of plugs ranging from 1–6. Where perithecia were present we targeted plugs to include perithecia that were accessible from the ground (0–2 m). Samples were stored on ice and then frozen within 48 h of sampling and stored at –20°C until processed for DNA extraction.

TABLE 2

Table 2. Site locations, sampling effort, and duration of BBD infection in terms of years since beech scale first observed.

We used the PRISM dataset (4 km²-daily resolution, PRISM Climate Group, 2020) to calculate climate variables for each site. We first determined the start and end dates of the growing season in each year based on empirical values describing heat accumulation for American beech leaf out and leaf drop (Richardson et al., 2006), wherein bud break for a given site and year was estimated as the date when a site had accumulated 100 cumulative GDD₄ (base 4°C) from January 1. Leaf drop was defined by 500 cumulative chilling degree days (below 20°C) from August 15. These calculations resulted in leaf out estimates ranging from March 15 to May 12 and leaf fall estimates from October 20 to November 8 along the natural climate gradient among sites. We considered non-growing season climate because fungi are likely to grow in periods where minimum temperatures are non-limiting, and growth during periods of tree host dormancy may be important for fungal establishment, growth, and/or aggressiveness. We then summed GDD₄, daily precipitation, and freeze-thaw frequency (the number of days that temperatures crossed 0°C) for both the growing season and non-growing season.

DNA Extraction, PCR, and Sequencing

Phloem plugs were prepared for DNA extraction by removing surface debris, Neonectria perithecia, and the periderm layer using a sterile scalpel that was cleaned with 70% ethanol and flamed-sterilized between samples. Plugs were then washed under a steady stream of 1 ml sterile 1x PBS pH 7.2 buffer. Phloem plugs were freeze-dried for 24 h then crushed and homogenized, and a subsample of phloem (mean 56 mg ± 18 mg s.d.) was subjected to bead bashing. DNA was extracted from ground phloem samples using a QIAGEN DNeasy Plant Mini Kit following the factory protocol. Extracted DNA was then purified using a Zymo OneStep PCR Inhibitor Removal Kit following the factory protocol and diluted 1:10 in PCR grade water. Positive and negative controls for Neonectria detection were created by collecting phloem from an asymptomatic beech tree in Durham, NH. The phloem was prepared following the protocol above, but with the addition of three cycles of autoclaving at 121°C for 30 min prior to the freeze-drying step. Negative controls were processed through DNA extraction as described above, whereas for positive controls 1–3 mg of mycelium from axenic cultures of N. ditissima and N. faginata was added to phloem subsamples prior to the bead-bashing step. Two replicates of each positive and three replicates of negative controls were included as separate samples in downstream PCR and sequencing.

The ITS2 region was amplified in duplicate PCR reactions with the primers 5.8S-Fun and ITS4-Fun using Phusion High Fidelity polymerase, a 58°C annealing temperature, and 30 PCR cycles (Taylor et al., 2016). Primers included the Illumina TruSeq adapters and sample identification tags were added to amplicons via a second round PCR at the University of New Hampshire Hubbard Center for Genome Studies. A dual-unique indexing strategy was used such that each sample had a matching pair of indices on forward and reverse reads in order to reduce potential for sample misassignment due to index switching or index bleed. We included PCR negatives (one per PCR plate) and DNA extraction negatives (one per kit) in the sequencing run.

Bioinformatic Analyses

Amplicon sequence variant (ASV) calling was performed using the DADA2 (Callahan et al., 2016) protocol v1.8 for ITS sequences in R 3.6.2 (DADA2 version 1.14.0), with an additional step to extract the ITS2 region from 5.8S-Fun–ITS4-Fun amplicon sequences using itsxpress (Rivers et al., 2018). The core DADA2 algorithm was run with option ‘pool = TRUE’ to allow greater sensitivity for ASVs that were rare in a single sample but more abundant across the entire dataset, and putative chimeras were removed using DADA2::removeBimeraDenovo. Sequences from trees with multiple plugs were pooled at the tree level before ASV calling. Taxonomy was assigned to ASV representative sequences by comparison to the UNITE dynamically clustered database (release date 04/02/2020; Nilsson et al., 2019) using the DADA2::assignTaxonomy algorithm.

We used the LULU post-processing algorithm to group putatively erroneous ASVs with their parent ASVs based on sequence similarity and co-occurrence patterns (Frøslev et al., 2017), which has been shown to improve reconstruction of biological taxa for fungi relative to ASV denoising alone (Pauvert et al., 2019). We note that this is not a clustering algorithm but rather uses minimum sequence similarity as the first of three criteria for culling of likely error variants. After comparing a range of minimum sequence similarity thresholds (84, 90, 93, 95%) we selected a 93% threshold, which maximized the identification and removal of likely erroneous ASVs while minimizing taxonomic reassignment.

One of the goals of this study was to determine the distribution of N. faginata and N. ditissima across the range of BBD. In order to increase sensitivity for the detection, after calling ASVs we subsequently mapped the original quality filtered reads to representative sequences for ASVs identified as N. faginata and N. ditissima using the ‘vsearch usearch_global‘ algorithm (Rognes et al., 2016). This approach increases detection sensitivity by including reads that may have otherwise been discarded or misassigned during ASV calling steps (Edgar, 2013; Pauvert et al., 2019). Samples were scored as containing N. faginata or N. ditissima if the species were discovered in a sample using either the ASV calling or mapping-based approaches. Singleton ASVs were removed from the dataset and samples with less than 1000 sequences after singleton filtering were also excluded.

Statistical Analyses

We first performed pairwise correlations of site characteristics data to examine relationships between disease severity and climate variables using the R package Hmisc (Harrell et al., 2019). These were either means of disease severity variables collected from random transects, or 10-year mean climate variables for each site extracted from the PRISM database. We first tested for normality using a Shapiro–Wilk test (Shapiro and Wilk, 1965). We report Spearman rank correlation (ρ) where one or both variables were non-normal (including all ordinal variables) – for other variables we report Pearson correlation (r).

To test for deviations from random patterns of co-occurrence of the two Neonectria species we used a probabilistic model (Veech, 2013; Griffith et al., 2016). In brief, all possible permutations of species occurrence were determined based on the number of samples and observed species frequencies. A probability distribution was then calculated wherein the probability of a given co-occurrence frequency was equal to the number of permutations with that co-occurrence frequency as a proportion of the total possible permutations. Significant deviations from random were then assessed by comparing the observed co-occurrence to the probability distribution with P equal to the sum of probabilities for co-occurrence in less than (P_lt) or greater than (P_gt) the observed number of samples. We excluded samples from the Wisconsin site, which at 9 years post-beech scale colonization had a considerable number of uninfected trees that would have skewed the analysis toward detecting aggregation.

To examine the effects of environmental covariates on N. faginata and N. ditissima occurrence, we applied spatially explicit joint species distribution modeling implemented in the R HMSC package (Ovaskainen et al., 2016, 2017; Tikhonov et al., 2017). We used the N. faginata and N. ditissima presence-absence matrix as the dependent variables and employed a probit link function. Pairwise site geographic distance was included as a random effect to control for spatial correlations between predictors and Neonectria occurrence. Tree was included as nested random effect (within site); natural log of sequence count was included as a fixed effect to account for sampling depth effects on species detection probability (Ovaskainen et al., 2017). Non-growing season climate variables were stronger predictors of Neonectria occurrence in exploratory analyses, and we therefore only considered non-growing season variables in subsequent models. The primary model included climate variables and disease severity variables in order to determine the best predictors of N. faginata and N. ditissima occurrence. Independent variables were standardized to their respective means and standard deviations to obtain comparable slope estimates. We report Tjur’s R², a coefficient of determination for logistic regression (Tjur, 2009; Ovaskainen et al., 2017) along with slope coefficients. To minimize issues with multicollinearity, precipitation was excluded from the model due to strong correlations with beech scale density (ρ = –0.88, P < 0.001) and DBH (ρ = –0.73, P = 0.02, Supplementary Table 1).

We used PERMANOVA (vegan::adonis function, Oksanen et al., 2019) to determine whether Neonectria species occurrences were associated with differences in composition of the remaining community. We used the presence-absence of each ASV with greater than 10% frequency as predictors of community composition, by first randomly subsampling the dataset to 1000 sequences per sample, iteratively removing each predictor ASV from the dependent variable matrix, and then performing PERMANOVA on transformed sequence counts (log₁₀ + 1) with the removed ASV as a categorical predictor. We tested for a relationship between ASV frequency and its strength as a predictor of community composition by regressing PERMANOVA R² against ASV sample incidence.

We next used indicator species analysis to explore whether certain species were associated with Neonectria species occurrence, or ordinal measures of Neonectria perithecia density, beech scale density, crown dieback, or cankering using the indicspecies::multipatt function (De Caceres and Legendre, 2009). We chose these variables because in combination they describe various stages of tree decline and aspects of disease and are amenable to transformation to categorical predictors (and therefore appropriate for ISA analysis). We used the sample-ASV matrix, randomly selecting 1000 sequences per sample for this analysis. We then identified ASVs that were indicators of at least two measures of disease severity. The goal of filtering out ASVs that were indicators of only one disease severity measure was to increase interpretability and reduce misclassification of ecological roles due to spurious correlations. We then performed functional classifications of ASVs using the FungalTraits database (Põlme et al., 2021). We noted primary and secondary lifestyles, and where additional functional potential was noted, such as endophytic capacity, we recorded this as tertiary lifestyle. We also noted “wood saprotroph” as a separate category, but collapsed other saprotrophic categories (i.e., “unspecified saprotroph,” “litter saprotroph,” “soil saprotroph”) into a single “saprotroph” designation.

Data Accessibility

Raw sequence data and associated sample metadata are archived at the National Center for Biotechnology Information Sequence Read Archive under BioProject accession PRJNA701888.

Results

Disease Severity and Site Characteristics

We performed a pairwise cross-correlation analysis using site-level means of disease severity variables and climactic variables (Supplementary Table 1). As expected, there were strong correlations between climate variables and latitude. Specifically, latitude was negatively correlated with heat accumulation (GDD₄) in both the growing (r = –0.64, P = 0.047) and non-growing season (r = –0.95, P < 0.001) and with precipitation in the growing season (Spearman’s ρ = –0.96, P < 0.001). Latitude was also negatively correlated with elevation (r = –0.73, P = 0.016) with more southerly sites tending to be at higher elevation. Longitude was not significantly correlated with the climate variables tested but was strongly positively correlated with duration of BBD infection (r = 0.97, P < 0.001) and cankering (ρ = 0.83, P = 0.003) and negatively with elevation (r = –0.69, P = 0.03). Various climate parameters correlated with disease severity. Non-growing season precipitation was negatively correlated with both scale insect density (ρ = –0.88, P = 0.001) and DBH (ρ = –0.73, P = 0.02), while growing season precipitation was negatively correlated with scale insect density (ρ = –0.81, P = 0.005). We also found correlations among some of the disease severity indicators we measured. Duration of infection correlated positively with cankering (ρ = 0.81, P = 0.004), while scale insect density and DBH were also positively correlated (ρ = 0.65, P = 0.043). No other correlations with disease severity indicators were found. We note, however, that the youngest site (Wisconsin, infested in 2010) had high mean wax density but no visible Neonectria perithecia at the time of sampling. We performed pairwise partial-correlation analysis using non-growing season precipitation, DBH, and beech scale while controlling for the effect of the third of these variables in each pairwise combination. Partial correlation analysis showed that only non-growing season precipitation and beech scale retained a significant correlation after controlling for the effect of DBH (ρ_{precip,scale.DBH} = –0.78, P = 0.013), whereas the other relationships were no longer significant (ρ_{scale,DBH.precip} = 0.12, P = 0.97; ρ_{precip,DBH.scale} = –0.45, P = 0.22).

Sequencing Results

After sequence processing and ASV denoising the number of sequences per sample (tree) ranged from <100 to 331,624 (median 21,315) across 117 samples. The LULU post-processing algorithm resulted in the identification of 149 ASVs (13% of 1132 original ASVs) that were better interpreted as variants of existing “parent” ASV’s (either due to minor sequence variation or sequencing error; Supplementary Table 2). Over 90% of the ASVs culled in this way shared full taxonomic identity with their respective parent ASVs. For example, the dominant N. faginata and N. ditissima ASVs each subsumed 8 and 3 “children” variants respectively with no taxonomic reassignment. Fourteen ASVs were taxonomically reassigned, but in all cases this involved reassignment to a higher taxonomic rank (i.e., species designations were removed but the genus [13 ASVs] or family [1 ASV] was retained). After LULU post-processing and removing samples with fewer than 1000 sequences 796 ASVs remained across 102 samples retained of the 117 original samples, with the number of ASVs per sample ranging from 12 to 172 (median 60). Rarefying (to 1000 sequences per sample) led to a further drop in ASV retention; ASVs richness ranged from 7 to 67 per sample (median 25). The mean ASV richness per site ranged from 13.2 ± 4.4 s.d. to 46.8 ± 16.4 s.d. and total site richness ranged from 66 to 412 ASVs. Further exploration of the relationships between disease severity, climate, and broad trends in community composition and diversity are possible using this dataset and are currently underway. We focus here on the primary fungal disease agents and description of species that may play a role in disease progression (Table 1).

We did not observe sequences in our DNA extraction negative or PCR negative controls (Supplementary Figure 1). Only one of the three negative controls composed of autoclaved phloem produced sequences – the three total sequences from that sample included two sequences from a Neocucurbitaria ASV and one sequence from a Corynespora ASV. Sequence counts from the nine samples of PBS buffer that had been used to wash plugs were likewise extremely low (1–8 sequences). Three of the PBS buffer wash samples contained N. faginata whereas N. ditissima was not detected in any of the samples. We included two positive controls for each of N. ditissima and N. faginata and observed the expected species and no other taxa in each of those samples, though one N. faginata positive control produced only a single sequence.

Neonectria Species Distribution and Its Drivers

Neonectria faginata was present in all 10 sites and N. ditissima was found in all but the southernmost site in North Carolina (Figure 1). Importantly, the two species often co-occurred both in sites and within the same tree, including within the same 1-cm phloem disc. Overall, N. faginata was present in 135 of 170 phloem plugs (79.4%) and 71 of 102 trees (69.6%), while N. ditissima was present in 46 of 170 plugs (27.1%) and 32 of 102 trees (31.4%). The two species co-occurred in 38 of 170 plugs (22.4%) and 27 of 102 trees (26.5%) including 35.5% of the 76 trees infected with at least one species. Both species were absent in 27 of 170 plugs (15.9%) and 26 of 102 trees (25.5%). At least one of the two species was detected in all 109 plugs where perithecia were present on the plug periderm surface prior to processing for DNA extraction (64.1% of plugs). Neonectria faginata was detected in 107 (98.2%) of the plugs with fruiting structures present while N. ditissima was detected in 31 plugs (29.4%). The species co-occurred in 29 of 109 plugs (26.6%). Sixty plugs had no perithecia present and one plug of the 170 total had degraded periderm such that it was not possible to record perithecia presence-absence. Of the 60 plugs with no perithecia at least one species was detected in 33 plugs (55%), wherein N. faginata was detected in 27 (45%) and N. ditissima in 15 (25%). The two species co-occurred in nine plugs without perithecia (15%), and were both absent in 27 plugs (45%).

FIGURE 1

Figure 1. Neonectria species occurrence across 10 sites. Proportions in pie charts indicate the number of trees in a site where a species was detected using metabarcoding. The number of trees per site is indicated by pie chart diameter as indicated by the inset scale legend.

Neonectria ditissima was only found in isolation at the tree level (i.e., without N. faginata) at the three northernmost sites we sampled (MI, WI, and northern ME). In the remaining seven of the 10 sites, N. ditissima was only detected in trees where N. faginata was also present. In terms of species prevalence, N. faginata occurred in a greater number of trees than N. ditissima in all but the three northernmost sites (90% mean occurrence versus 25% mean occurrence for N. faginata and N. ditissima, respectively, in the seven more southerly sites). The two species each occurred in isolation in two trees in Wisconsin and co-occurred in one tree, and in Michigan each species occurred in isolation in one tree and co-occurred in six trees. Neonectria ditissima was more prevalent in our northern Maine site (70% versus 60% of trees for N. ditissima and N. faginata, respectively, including 50% of trees where the species co-occurred). In sites within the killing front and aftermath zone (i.e., excluding the Wisconsin site within “advancing front” of the disease), N. faginata was detected in 84% of all trees across our sites, whereas N. ditissima was detected in 36% of trees. Given these observed occurrence frequencies, co-occurrence did not differ from a random distribution [32% observed vs. 30% expected; P_gt = 0.24 sensu Veech (2013)].

We next examined correlations between Neonectria species incidence across sites (i.e., presence-absence at the tree level) and indices of disease severity and climate using spatially explicit joint species distribution modeling (HMSC, Ovaskainen et al., 2017). We first tested for effects of growing season versus non-growing season climate parameters and found that climate during the non-growing season was overall a stronger predictor of patterns of Neonectria occurrence (Supplementary Figure 2). Specifically, heat accumulation (GDD₄) during the non-growing season was significantly associated with incidence of both species, while growing season GDD₄ was generally a poor predictor of incidence. Non-growing season freeze-thaw was significantly associated with N. faginata (posterior support P = 0.99) but not N. ditissima incidence (P = 0.86). In the growing season neither of these variables was correlated with Neonectria species occurrence.

After controlling for sampling effort and spatial structure, our models explained 54% and 37% of variance in N. faginata and N. ditissima incidence, respectively, based on Tjur R² (Tjur, 2009). The full model indicated that N. faginata had a significant, positive relationship with non-growing season heat accumulation (R² = 0.19), duration of infection (R² = 0.14) and DBH (R² = 0.04), and a negative association with beech scale density (R² = 0.06) (Figure 2). Heat accumulation was the strongest predictor of N. faginata incidence (R² = 0.19). Neonectria ditissima incidence was positively associated with DBH (R² = 0.07), freeze-thaw cycle frequency (R² = 0.07), and crown dieback (R² = 0.04), but negatively associated with non-growing season heat accumulation (R² = 0.07) and beech scale density (R² = 0.05). Results were unchanged when non-growing season precipitation was included in the model, except there was no significant correlation between N. faginata and beech scale or infection duration (Supplementary Figure 3). Neither species showed a significant correlation with non-growing season precipitation. The HMSC approach also allows examination of residual correlation between dependent variables (i.e., N. faginata and N. ditissima occurrence) after accounting for the effect of independent predictors. We found no residual correlation between the two species after accounting for the effects of disease severity and climate, and also found no correlation between the species using a reduced model that only controlled for sampling effort (i.e., sequence count) and spatial structure. This result suggests that distribution of the two species is not strongly structured by inter-species interactions (e.g., competition or facilitation).

FIGURE 2

Figure 2. Effects of disease severity and climate on Neonectria species occurrence within trees. Points with solid outlines represent significant effects at posterior probability P > 0.95 (analogous to a p-value of P < 0.05). Box fill indicates the strength and direction of the relationship with blue indicating a negative slope, and red indicating a positive slope. Point size indicates magnitude of R² with overall fits for N. faginata and N. ditissima of Tjur R² = 0.54 and 0.37, respectively.

In light of the significant positive association of N. ditissima with crown dieback indicated in HMSC modeling, we visually explored patterns of species incidence across levels of crown dieback and cankering (Figure 3). Neonectria ditissima occurrence increased across crown dieback classes (crown dieback level 0 = 21%, 1 = 24%, 2 = 38%, 3 = 58%), whereas N. faginata occurrence was relatively stable or modestly declined (crown dieback level 0 = 79%, 1 = 71%, 2 = 62%, 3 = 67%). The increase in N. ditissima with increasing crown dieback resulted in an increase in co-infection by the two species at the tree level. Neonectria ditissima occurred in isolation in 2–8% of trees in each dieback class, whereas trees co-infected with both Neonectria species increased across dieback classes (crown dieback level 0 = 16%, 1 = 21%, 2 = 31%, 3 = 50%). When discrete cankers were absent both species were at relatively low occurrence (38 and 26% for N. faginata and N. ditissima, respectively). However, N. faginata more than doubled in trees with cankers compared to no cankers (85% versus 38%, respectively; χ² = 21.6, P < 0.001), and N. ditissima doubled in the highest cankering level compared to levels 0–2 (50% versus 25%, respectively; χ² = 5.1, P = 0.02) resulting in elevated co-infection at the highest cankering category (46% of trees compared to 17–23% in remaining categories).

FIGURE 3

Figure 3. Occurrence of N. faginata and N. ditissima across levels of crown dieback (A) and cankering (B). The proportion of trees in which each species was detected is indicated for N. faginata (dark red) and N. ditissima (light red). The proportion of co-occurrence is indicated by overlap of the bars and the remaining white space indicates the proportion of trees where neither species was detected. The number of trees per grouping (n) is indicated above respective bars.

In this work we were also interested in whether N. faginata or N. ditissima was structuring or responding to different fungal endophyte communities. Of the 62 most commonly detected ASVs (minimum 10% incidence across all samples) N. faginata was the fifth best predictor of mycobiome composition (PERMANOVA R² = 0.076, P = 0.001) while N. ditissima was the 42nd best predictor (R² = 0.023, P = 0.007; Supplementary Figure 4A). We tested for a relationship between PERMANOVA R² and sample incidence to examine the possibility that species presence-absence effects on community composition were primarily driven by the frequency of ASV occurrence. There was a weak but significant relationship between incidence and PERMANOVA R² (R² = 0.165, P = 0.001; Supplementary Figure 4B) with the top five most predictive ASVs occurring in between 13 and 61% of samples.

Ecological Roles of Fungi in the BBD System

We used a combination of Indicator Species Analysis (ISA; De Caceres and Legendre, 2009) and literature-based functional classification (Põlme et al., 2021) to discern degrees of statistical association and potential ecological roles for fungal species that appear as beech bark endophytes across the range of disease and tree decline levels sampled. Overall, 38 ASVs were identified as indicator species of at least two disease categories. All but two of these ASVs were indicators of different levels of crown dieback or cankering (Table 3). Nine of the 38 ASVs were indicators of the absence of crown dieback and either low levels of cankering or the absence of scale insect (see Table 3, “Healthy beech” indicators). Of these nine, only three were taxonomically identified to the genus level, with six other ASVs identified to order, phylum or kingdom. Another five ASVs were indicators of both low levels of cankering and absence of scale insect (Table 3, “Minor cankering, scale absent”), with two being taxonomically identified to the genus level.

TABLE 3

Table 3. Amplicon sequence variant (ASVs) associated with at least two indicators of disease according to ISA [columns “Disease indicators (ISA)”].

In total, 16 ASVs were associated with intermediate to high levels of either cankering or crown dieback. Six ASVs were indicators of intermediate to high levels of cankering, low scale density, and/or presence of Neonectria species (Table 3, “Intermediate BBD severity”). These included N. faginata, as well as ASVs annotated as animal pathogens/entomopathogens, mycoparasites, plant pathogens, and saprotrophs or wood saprotrophs, the latter two functional groupings accounting for five of the six ASVs. Ten ASVs were indicators of intermediate to high levels of crown dieback (Table 3, “High BBD severity”). Four of these 10 were also indicators of high scale density, and another four were indicators of high levels of cankering. Eight of the 10 ASVs associated with high levels of crown dieback were annotated as either saprotrophs or wood saprotrophs, and five were annotated as plant pathogens. All five of the ASVs assigned a plant pathogen function mapped to multiple functional groups, with saprotrophic lifestyle assigned as primary or secondary functions. Three ASVs associated with high levels of crown dieback were also among the top five predictors of community composition (Supplementary Figure 3), and were annotated to entomopathogen (animal pathogen), plant pathogen-saprotroph-mycoparasite, or saprotroph-plant pathogen-endophyte functions, respectively. We note that N. ditissima was also an indicator of high levels of crown dieback (crown dieback 2–3) but was not formally included in this analysis due to lack of association with additional disease indicators.

Eight ASVs were associated with presence or absence of the primary disease agents (Table 3, “Scale and Neonectria associates”), including ASVs associated with high beech scale density and Neonectria absence (two ASVs), high density of both beech scale and Neonectria perithecia (one ASV), or absence of both beech scale and Neonectria (one ASV). Six of the eight ASVs were also associated with absence of cankering. Five of these eight ASVs were annotated to saprotroph or wood saprotroph functions, three as entomopathogens/animal pathogen, and two as endophytes, including three ASVs with multiple functional mappings.

We also specifically explored the distributions of Clonostachys rosea, Nematogonum ferrugineum, and Fusarium babinda given their previously described roles as BBD associates as either mycoparasites of Neonectria (C. rosea and N. ferrugineum; Barnett and Lilly, 1962; Houston, 1983; Stauder et al., 2020b) or potential entomopathogens or secondary beech pathogens (Stauder et al., 2020b). Two ASVs were annotated as C. rosea and together they occurred in four of our 10 sites ranging from 10 to 27% of trees in respective sites (Figure 4A). One of these ASVs was also an indicator of the highest level of Neonectria perithecia. One ASV of putative importance as a mycoparasite of the BBD fungi, N. ferrugineum, only occurred at two sites at low frequency (9 and 20% of trees Figure 4B) and was not a statistical indicator of any disease categories. Another ASV (ASV 19) was identified as Fusarium babinda, which has been previously identified in association with BBD and is a suspected beech scale associate (Stauder et al., 2020b). This ASV was associated with high wax density and high crown dieback (Table 3) and occurred in all 10 of our sites (Figure 4C). Of these three species only C. rosea showed a statistical association with either Neonectria species. Clonostachys rosea was positively associated with N. ditissima (Fisher’s exact test odds ratio = 6.3, P = 0.03) and only occurred in one sample where both species were absent (Supplementary Figure 5).

FIGURE 4

Figure 4. Occurrence of ASVs 16 and 21, Clonostachys rosea (A), ASV 807, Nematogonum ferrugineum (B), and ASV 19, Fusarium babinda (C), across 10 sites. Proportions in pie charts indicate the number of trees in a site where the species was detected using metabarcoding.

Discussion

In the present study we explored the distribution of N. faginata and N. ditissima, the primary pathogens involved in BBD, in relation to disease severity and climate characteristics in 10 sites across the range of BBD. We further explored patterns of association with other fungal taxa in the beech bark endophytic community in relation to tree disease state, and used these relationships to highlight and hypothesize ecological roles of potential relevance to BBD severity and progression. We show that N. faginata and N. ditissima have divergent correlations with climate and may be associated with different stages of tree decline. Further, we show that fungal taxa occurring in association with the primary BBD agents (i.e., the host, F. grandifolia; the scale insect initiating agent, C. fagisuga, and the fungal pathogens, N. faginata and N. ditissima) may contribute to disease outcomes in complex and interacting ways in this system. We suspect that these types of feedbacks are not unique to the BBD system, but rather are likely to generalize to other complex diseases where variation in the occurrence, relative densities, or interactions among ecologically redundant disease agents influence the trajectory of disease (Cobb and Metz, 2017). Outcomes are likely to be especially hard to predict where component species respond uniquely to key biotic and biophysical drivers.

Correlations Among Disease Agents and Climate

We found that precipitation in both growing and non-growing seasons was negatively correlated with scale insect density. Non-growing season precipitation also correlated negatively with DBH. Negative correlations between beech scale and precipitation are consistent with a hypothesized causal relationship whereby precipitation reduces scale insect population density by washing colonies off of tree boles (Houston and Valentine, 1988; Dukes et al., 2009; Garnas et al., 2011b; Kasson and Livingston, 2012). The latter relationship between precipitation and DBH is likely driven by multicollinearity between infection duration (and thus disease severity indicators), geographic distribution, and climate. Sites with older infections tend to be dominated by small diameter trees and also occur in lower precipitation sites in our dataset. That said cankering (which constitutes evidence of past, non-lethal infection) was the only index of disease status that was related statistically (positively in this case) with duration of infection. Our dataset included sites ranging from nine to 68 years from the arrival of beech scale based on county-level data (Cale et al., 2017), but only one site was in the “advancing front” stage of infection as typically defined (i.e., ≤10 years post scale insect arrival; Houston et al., 2005). At least one further site, and possibly up to three sites, occurred in the “killing front” stage of infection (i.e., 5–10 years after advancing front conditions, or 15–20 years after scale insect arrival; Houston et al., 2005), whereas the remaining sites were sampled in an aftermath forest infection stage. Interactions between disease agents are variable in aftermath forests despite density dependent growth within populations of insects or fungi (Garnas et al., 2011b). Our sampling design, which was weighted toward aftermath forest stands, may have obscured the typically observed patterns across disease stages in beech scale and fungal pathogen abundance. However, there were signals of duration of infection-driven patterns in our dataset. For example, our advancing front site had high mean wax density but no visible Neonectria perithecia production, as is typical of advancing front stage forests (Shigo, 1972). Further, scale insect density and DBH were positively correlated, likely reflecting both a positive relationship at the tree scale (Latty et al., 2003; Garnas et al., 2011a) or an effect of our sampling design. Specifically, since we included stands within the advancing front, sites with higher beech scale density, as is common in recently infested stands, also tended to have larger trees prior to experiencing BBD-induced mortality. Thus, inclusion of both aftermath and advancing front sites is likely to over-represent the scale insect-tree size relationship, which is considerably weaker within the aftermath forest alone (Garnas et al., 2013).

Neonectria Species Distribution and Its Drivers

Both Neonectria species were widely distributed geographically and in terms of disease stage (i.e., infection duration). Neonectria faginata was detected in all 10 sites, and N. ditissima in nine of 10 sites, with both occurring in advancing front, killing front, and aftermath forests. Importantly, the two species regularly co-occurred, not only within the same tree (26% co-occurrence), but within the same 1-cm phloem disk (22% co-occurrence). The prevailing understanding of the dynamics of BBD is that infected stands undergo a progression from initial infection by N. ditissima to near-total replacement by N. faginata in later disease stages in most cases (Houston, 1994; Cale et al., 2017). This idea persists despite evidence of N. ditissima persistence and co-occurrence of the two pathogens (Kasson and Livingston, 2009). Here we show that while N. faginata is the dominant species throughout killing front and aftermath forests (84% of all trees), N. ditissima maintained a substantial foothold throughout these stands (36% of all trees). Such relatively high incidence indicates that spillover of N. ditissima to non-beech trees could be an important mechanism by which this disease impacts forest structure, function, and diversity. Further, given random co-occurrence patterns between N. faginata and N. ditissima, we found no evidence of strong facilitation such that co-infection frequency would be elevated nor obvious competitive exclusion within sites or trees.

Despite approximately random co-occurrence between the two species we did observe patterns with regard to species distributions, stand characteristics, and climate. For example, in accordance with previous research (Houston, 1994) we found evidence that N. faginata becomes more prevalent as disease progresses over decades. Generally, the two species appear to have divergent climate associations, where N. faginata was associated with warmer climates, and N. ditissima with colder climates. It is possible that these patterns arise from infection duration dynamics that are obscured by our use of coarse-grained county level data. For example, northern Maine has experienced secondary killing fronts after release of beech scale populations from suppression by winter killing temperatures as temperature warms (Kasson and Livingston, 2012). However, the two warmest sites in our dataset were also of intermediate infection duration (19 and 20 years) and so infection duration is unlikely to fully explain the climate-driven distribution patterns we observed. Indeed, non-growing season heat accumulation was the strongest predictor of N. faginata occurrence while infection duration was the second strongest predictor, suggesting climate as an important influence on prevalence of these species within the BBD system.

Both Neonectria species exhibited significant associations with various disease severity metrics. For example, both species were negatively correlated with C. fagisuga density. Despite the fact that scale insects appear to be obligate initiating agents of BBD in the advancing front, correlations between insect and fungal components of BBD are weak or negative in aftermath forests (Cale et al., 2012; Garnas et al., 2013). Our data support this pattern. Alternative hypotheses have been offered for the priming of trees for fungal infection, including a role for the native birch margodid scale insect, Xylococculus betulae, as a predisposing agent specifically for N. ditissima (Cale et al., 2015b). While this insect can be locally abundant in certain sites and years, we did not encounter it in sufficient densities to warrant inclusion in our models. Nor does it seem likely that X. betulae is an important driver of fungal dynamics across the range of BBD.

Occurrence of Neonectria was most strongly associated with tree condition and with aspects of climate. Neonectria ditissima but not N. faginata was positively correlated with crown dieback class, a result supported by both HMSC modeling and indicator species analysis (ISA). This pattern resulted in increased co-infection by the two species in higher crown-dieback classes. Higher rates of co-infection associated with increasing crown dieback may indicate that a synergistic attack by the two species contributes to tree decline. For example, N. ditissima showed higher prevalence in sites with lower temperatures during the non-growing season. Maintaining growth during the non-growing season, however slow, could result in higher apparent aggressiveness if attacking dormant hosts limits host defensive response or wound compartmentalization (Manion, 2003; Dukes et al., 2009; Copini et al., 2014). Alternatively, N. ditissima may be more prevalent in later stages of tree decline as a secondary pathogen that is favored by weakened host tissue (Houston, 1981; Manion, 1981).

There is some evidence that trees infected with heart-rot decay fungi have higher density of Neonectria-induced lesions, potentially indicating that modification of tree tissues by either Neonectria pathogens or decay fungi facilitates host colonization (Cale et al., 2015a). However, the directionality or strength of this relationship is unclear. Indeed, N. faginata was among the top five predictors of community composition of the bark endophyte community overall, whereas N. ditissima was a relatively poor predictor, suggesting that colonization by N. faginata may be a predisposing factor for colonization of bark by a suite of other disease-associated fungi. However, the two species produced similar sized cankers on American beech in inoculation trials (Stauder et al., 2020b) suggesting that N. ditissima is not restricted to tissues already weakened by primary N. faginata infection.

Since our results suggest a much more central role for N. ditissima in the aftermath forest than has been previous indicated, it is important to attempt to reconcile differences between our study and previous work on BBD-associated Neonectria species. For example, in a recent study, N. ditissima represented just 4.2% of perithecial isolates from American beech across the central Appalachian Mountains and was recovered from only two of 13 sampled locations (Stauder et al., 2020b). This is considerably lower than the 27% of bark samples we estimate here. Most studies to date rely on culturing from ascospores or spore measurement as direct isolation from bark tissue is challenging. Higher rates of perithecium production, seasonal or within-tree differences in sporocarp production and/or spore viability could favor the collection or successful rearing of N. faginata over N. ditissima, particular in cases of co-infection. While all sampling methodologies likely suffer from some form of bias, our metabarcoding approach using DNA extracted directly from bark tissue likely represents an improvement on historical detection of BBD fungal agents and associates.

Ecological Roles of Fungi in the BBD System

We used a combination of ISA and literature searches to describe the ecological roles of fungi occurring in the context of the BBD complex. ISA delineated clear groups of fungi associated with different stages of disease or disease agents, including healthy beech associates, fungi associated with intermediate or high BBD severities and associated tree decline, and fungi associated with presence and/or absence of the primary disease agents (i.e., Neonectria and beech scale).

The majority of healthy beech associates were not taxonomically identified past the order level (67%). In addition, only 40% of ASVs associated with low levels of cankering and absence of beech scale were identified past the order level. Of the 24 remaining indicator ASVs 96% (23 ASVs) were taxonomically identified to at least the genus level, suggesting that bark endophytic communities of healthy American beech are a relatively unexplored reservoir of fungal diversity given comparatively low taxonomic identification. Bark endophytes have historically received less attention than foliar endophytes in temperate forests (Unterseher, 2011), and further exploration of bark endophyte diversity and functioning is warranted.

We observed a shift in both taxonomic composition and functional potential in the indicators of intermediate to high BBD severity compared to healthy beech. Overall, 13 of the 16 ASVs (81%) associated with intermediate to high BBD severity were annotated to saprotrophic functional groups, and five of the putatively saprotrophic ASVs associated with high BBD severity were further identified as facultative plant pathogens. Four ASVs in intermediate- and high BBD severity categories, including N. faginata, were also among the top five predictors of community composition. Many endophytes and plant pathogens function as facultative saprotrophs (Frankland, 1998; Stone et al., 2004); a shift in community function toward saprotrophy may be an important indicator of later stages of tree decline. Decay fungi may have negative and weak positive associations with beech scale and Neonectria, respectively (Cale et al., 2015a), suggesting complex interactions between disease organisms and co-occurring fungi potentially mediated by host tissue modification. The fungal communities associated with late stages of tree decline in particular, as indicated by high levels of crown dieback, may contribute to tree death by weakening tissues to the point of mechanical failure (i.e., “beech snap”; Houston, 1981; Manion, 1981). Together, enrichment of saprotrophs and plant pathogens along with an apparent consistent shift in community composition indicate that the endophytic fungal community may play an important role in disease progression beyond the direct action of the primary disease agents.

We observed eight ASVs that were associated with presence or absence of the primary disease agents (Neonectria and beech scale). It is difficult to assign ecological roles or the nature of interactions based on pairwise species associations. For example, a positive correlation between species may indicate facilitation, overlapping habitat or microclimatic preferences, or may indicate deadlocked competition or mycoparasitic relationships (Maynard et al., 2018). Indeed, some of the ASVs in this group were associated with high beech scale density and absence of Neonectria, or vice versa. Two of these (ASVs 27 and 234) were restricted geographically – occurring primarily in our Wisconsin site, which, at 9 years infection duration, was unique in our dataset with high average beech scale density and low Neonectria incidence. As such, these two ASVs may be indicators of early stage BBD infection or geographic location rather than of interactions with beech scale or Neonectria, per se. However, two other ASVs (ASVs 69 and 217) that were geographically widespread, occurring in 4-5 sites, were both indicators of beech scale absence, with one also being an indicator of high Neonectria perithecium density and the other an indicator of Neonectria absence. It is possible that these taxa function as facultative entomopathogens and/or mycoparasites depending on BBD disease stage and available hosts.

We also examined distribution of two species previously reported as mycoparasites (N. ferrugineum and C. rosea) and one reported as a potential additional plant pathogen and/or an entomopathogen (F. babinda) in the BBD system. One of these, N. ferrugineum (ASV 807) occurred at only two sites and was not a statistical indicator of any disease categories. Despite its prevalence in visual surveys (Houston, 1983) our data suggest that given its relative rarity, this fungus may not play a primary role in limiting growth of Neonectria pathogens involved in BBD. In contrast, two ASVs were identified as C. rosea, which together were detected at four of our 10 sites including recently infested sites in Wisconsin and Michigan, suggesting this fungus is present in BBD-affected forest stands even at the earliest stages of Neonectria establishment. In addition, one of the C. rosea ASVs was an indicator of high Neonectria perithecia density, consistent with its hypothesized mycoparasitic status in this system (Stauder et al., 2020b). This species has long been used as a biocontrol agent against plant pathogenic fungi (Schroers et al., 1999). Genomic mechanisms of mycoparasitism have also been described in this species (Karlsson et al., 2015) making this an intriguing candidate for further study in terms of interactions with the primary Neonectria disease agents. Fusarium babinda (ASV 19) was an indicator of high BBD severity and was found at all 10 of our sites, suggesting this fungus is a geographically widespread and consistent member of late-stage decline communities associated with BBD. In particular, F. babinda was associated with high wax density supporting its hypothesized association with beech scale (Stauder et al., 2020b). Whether this fungus is parasitizing the scale insects remains unclear. Previous literature supports a possible entomopathogenic lifestyle, having been previously recovered from other non-native forest insect pests (Lymantria dispar and Adelges tsugae) in the eastern U.S. (Jacobs-Venter et al., 2018).

Conclusion

This study explores new aspects of the fungal community ecology of the BBD system. We combined structured sampling at a landscape scale and across a known gradient of duration of infestation with beech scale with robust characterization of the beech bark mycobiome using high-throughput, metabarcoding sequencing. We found that one of the primary BBD agents, N. ditissima, occurs far more widely than previously known, co-occurring with N. faginata in nearly all of the sites and regularly within the same tree and even the same 1-cm phloem disk. The two species have apparently contrasting climate associations with N. faginata being associated with warmer temperatures and N. ditissima with cooler temperatures. Further, the two species are non-randomly associated with different stages of disease. While N. faginata is indeed the dominant pathogen, N. ditissima becomes most prevalent both very early during stand BBD-level infection and also at later stages of tree decline within the aftermath forest. The role of this pathogen, both in conjunction N. faginata and as trees progress through stages of the decline cycle, should be further explored. We also identified fungi that regularly occur in association with BBD with the potential to influence disease trajectory by functioning as entomopathogens, mycoparasites, saprotrophs and/or alternate or additional pathogens, and resulting in downstream shifts in community composition in the fungal communities of beech bark.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm.nih.gov/, PRJNA701888.

Author Contributions

JG and EM designed the study. JH and EM collected and processed the samples. EM and JG designed the data analysis approach. EM performed the data analysis and wrote the first draft of the manuscript. EM, JG, and MK wrote sections of the manuscript. All the authors contributed to manuscript revision, read, and approved the submitted version.

Funding

This work was funded in part by the New Hampshire Agricultural Experiment Station including support for open-access publishing fees. This is Scientific Contribution Number 2895. This work was supported by the USDA National Institute of Food and Agriculture McIntire-Stennis Project 1023443.

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.

Acknowledgments

We thank Paul Super and Bambi Teague at the National Parks Service in North Carolina as well as Keith Kanoti at the University of Maine for their help in obtaining collection permits. We appreciate the cooperation of the New York Department of Environmental Conservation, and the Michigan, Pennsylvania, and Wisconsin Departments of Natural Resources for granting collection permits.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/ffgc.2021.673099/full#supplementary-material

References

Barnett, H. L., and Lilly, V. G. (1962). A Destructive Mycoparasite, Gliocladium roseum. Mycologia 54, 72–77. doi: 10.1080/00275514.1962.12024980