Inter-annual Persistence of Canopy Fungi Driven by Abundance Despite High Spatial Turnover

While it is now well established that fungal community composition varies spatially at a variety of scales, temporal turnover of fungi is less well understood. Here we studied inter-annual community compositional changes of fungi in a rainforest tree canopy environment. We tracked fungal community shifts over 3 years in three substrate types (live bryophytes, dead bryophytes, and host tree bark) and compared these changes to amounts of community turnover seen at small spatial scales in the same system. The effect of substrate type on fungal community composition was stronger than that of sampling year, which was very small but significant. Although levels of temporal turnover varied among substrates, with greater turnover in live bryophytes than other substrates, the amount of turnover from year to year was comparable to what is seen at spatial distances between 5 and 9 cm for the same substrate. Stability of communities was largely driven by a few fungi with high relative abundances. A majority of fungal occurrences were at low relative abundances (≤ 0.1%). These fungi tended to be short lived and persisted to following years ≤ 50% of the time, depending on substrate. Their presence and persistence are likely impacted by stochastic processes like dispersal limitation and disturbance. Most samples contained only one or a few fungi at high relative abundance (≥ 10%) that persisted half or more of the time. These more abundant and persistent fungi are expected to have sustained functional interactions within the canopy ecosystem.


Introduction
Molecular surveys of fungi in the environment have revealed diverse communities with high beta diversity and heterogeneity at a range of spatial scales, from the local to continental [1][2][3][4]. While spatial turnover tends to be high, it is largely unknown how stable these communities are over time and how temporal turnover compares to changes across space. Studies of fungal communities typically only sample at a single time point and do not capture temporal variation or consider how it may impact community dynamics or its ecological functions. It is also unclear if a sample taken at one time point can be used to predict the composition of fungal communities into the future.
Fungal genets, in the form of continuous vegetative mycelia, have the potential to live for centuries [5], which could promote long-term community stability. Turnover of fungal structures like rhizomorphs and hyphae, however, can occur over the span of days to months [6][7][8][9]. Because of this, individual fungi might become extirpated from local areas over short time scales. Mycelial death coupled with 1 3 continuous arrival of new spores, which themselves vary with time [10], or hyphae from the surroundings create the possibility for dynamic fungal communities that change rapidly with little consistency over time. The level of stability in fungal communities has implications for their ecological functions. High community turnover could result in continuously changing functions that cannot be predicted into the future. While functional redundancy among fungi has been demonstrated for some traits [11], different taxa also appear to often have unique functions. These include enzymatic activity to mobilize nutrients from complex organic material [12] or host-specific interactions of pathogens [13] or mycorrhizal fungi [14][15][16][17]. Alternatively, stable communities would allow ecological functions to persist over time, creating localized "hot-spots" for particular functions.
Longevity and temporal persistence vary among fungal taxa [18]. Recent research has found different life history strategies among fungi within the same guild, such as in wood decay fungi [19,20] and mycorrhizal fungi [21][22][23], in which fungi differentially invest in biomass, stress-tolerance, or reproduction. These strategies and traits impact how fungal taxa differently respond to and tolerate environmental conditions like disturbance [22] and inter-specific competition [24]. These biological differences can allow some fungal taxa to persist long term, while others appear and disappear rapidly. Most studies of fungal temporal stability and turnover have focused on seasonal trends in soils. Fungal community composition [25][26][27], biomass [28], and enzymatic activity [28] can all vary intra-annually among seasons. While seasonal comparisons have generally represented temperate or boreal regions with pronounced seasonality, differences in fungal biomass and community composition can also be seen in tropical areas between wetter and drier parts of the year [29,30]. When these studies are conducted over multiple years, they typically find that intra-annual seasonal changes in community composition exceed changes seen among years in the same season within a site [26,31].
Another way to understand temporal turnover is to compare it to turnover in space. While this is seldom done, it allows insights into how spatial and temporal dynamics separately and together influence fungal communities. When assessing seasonal changes in fungal communities at the continental scale, Averill et al. [31] found that the difference between seasons at one site could be equivalent to the differences between sites thousands of kilometers apart, suggesting strong seasonal effects on fungi and high community turnover between seasons. At a smaller spatial scale, Izzo et al. [32] found inter-annual temporal changes in temperate forest ectomycorrhizal fungi colonizing roots were small compared to differences in the community seen across space at the same site. Plots were more similar to each other between years than to nearby plots within the same year.
Dominant fungi tended to be more constant among years, suggesting that temporal stability could be driven by only a subset of the fungal community, while other fungi appear and disappear erratically.
In the present study, we investigated inter-annual turnover of hyperdiverse fungal communities in an epiphytic tree canopy environment. In previous work in this system, we found high spatial turnover at very small, sub-meter distances and among three closely intertwined substrates [33]. We resampled at the same locations for two additional years, with 3 years of collections in total, to quantify temporal changes in fungal communities in three different substrates. We statistically examined the contributions of sampling year, substrate type, and spatial location to explaining community composition. We also compared the observed inter-annual temporal turnover to intra-annual spatial turnover to determine how much spatial distance is needed to reach equivalent levels of community change as is seen with time. Additionally, we examined how abundances of individual fungi are related to their persistence in subsequent years.

Sample Collection
We performed this study in Parque Nacional Tapantí in Cartago Province, Costa Rica. The sampling site was located along the east bank of the Rio Orosi (9.742°N, 83.784°W, 1300 m elevation) and consisted of five branches on three Saurauia montana (Seem., Actinidaceae) trees. In 2015, we sampled 135 locations along these branches with a 9-mm diameter borer in a study of fungal community spatial turnover and substrate specificity. After collection, each core was dissected into up to three substrate types: live bryophytes, dead bryophytes, and host tree bark. We rinsed substrate samples in sterile water and removed any vascular plant materials. Sampling locations were situated between 1 cm and over 8 m apart from each other to capture community turnover at a range of small spatial scales. The sampling design is more thoroughly described in [33]. We returned in 2016 and 2017 and resampled immediately adjacent (~ 1 cm away) to 59 of these locations. Sampling in all years took place between July 8 and 16 during the rainy season. This period is characterized by a mean daily relative humidity of > 90% and mean temperatures between 16 and 19.5 C (Fig. S1). All resampled locations were spaced at least 15 cm apart to minimize autocorrelation among points due to spatial proximity. Some sampling locations were lost or destroyed due to flooding and branch breakage and were only sampled for 2 out of the 3 years. In 2015, samples were preserved in RNALater (Ambion, ThermoFisher) and later lyophilized upon return to the US, while in 2016 and 2017 they were frozen after dissection and lyophilized within a week of collection.

Bioinformatics
We merged paired end reads with USEARCH v9.2.64 [35], excluding merged sequences less than 150 bp and sequences that did not merge. Primer sequences were removed with cutadapt [36]. We then quality filtered and dereplicated sequences with USEARCH using program defaults. We clustered OTUs and removed chimeric OTUs with the clus-ter_otus command. Only OTUs with at least two sequences were retained. We constructed a table of OTUs clustered at 97% with the usearch_global command. We made taxonomic assignments using the RDP Classifier [37] and the UNITE database (v7.2 QIIME release 2017-10-10 [38]), keeping only assignments at 80% or above. All OTUs not assigned to the kingdom Fungi were removed from the dataset as likely non-fungal taxa, and samples with fewer than 1000 remaining reads were also excluded. Only 27.7% of remaining OTUs could be identified to the genus level.

Effects of Year, Space, and Substrate on Fungal Communities
To test the effects of sampling year, spatial location, and substrate type on fungal community composition, we used redundancy analysis (RDA). These analyses were done with both fungal presence-absence and transformed relative abundance data. We included spatial relationships among samples by creating spatial vectors using principal coordinates of neighbor matrices (PCNM) as implemented in the vegan package v2.5-4 [40]. In the matrix of distances, different substrates collected at the same sampling location at the same time were treated as being 0.1 cm apart, and samples collected at the same location but in different years were treated as being 1 cm apart. The spatial vectors were forward-selected with the Blanchet et al. [41] stopping criterion. We partitioned the variance in fungal community composition among year, substrate type, and spatial vectors, and tested significance using permutation tests with 999 permutations.
We visualized the effects of sampling year and substrate type on fungal community using two-dimensional nonmetric multidimensional scaling (NMDS) with relative abundance data and the Bray-Curtis dissimilarity metric. We also performed NMDS for each of the three substrate types separately to examine inter-annual variation within a given substrate. NMDS was done in the phyloseq package v1.24.2 [42] with 100 iterations. To determine if OTUs that were common (found in many samples) in 1 year tended to be common in other years, we tested Spearman correlations between percent of samples each OTU was found in between pairs of years.

Alpha Diversity
We analyzed the Shannon index and log-transformed expected OTU richness using a linear mixed-effects model in the lme4 package v1.19-1 [43]. We calculated the expected OTU richness of each sample at 1000 reads of sequencing depth with the rarefy function in vegan [40]. For the Shannon index, we rarefied each sample down to 1000 reads, calculated the index, repeated this process 1000 times, and took the average Shannon diversity for each sample. The model included substrate type, sampling year, and their interaction as fixed effects. Sampling location and substrate sample nested within sampling location (ex. bark substrate at sampling point 1) were included as random effects. We tested the significance of each fixed effect with Wald chi-square tests and did pairwise comparisons between years and substrates using the emmeans package v1.3.3 [44].

Temporal Turnover
We assessed overall inter-annual change in fungal community within each sampling location by calculating the Bray-Curtis and Jaccard dissimilarities between pairs of years for each substrate sample collected at a given location. Jaccard dissimilarity was partitioned into nestedness and turnover components according to Baselga [45] in the betapart package v1.5.1 [46]. We used a linear mixedeffects model as described above for alpha diversity to test the effects of year, substrate, and their interaction on interannual community dissimilarity. For each substrate type and pair of years, we calculated the average Bray-Curtis and Jaccard dissimilarities among pairs of samples from the same sampling location collected in different years. We compared these dissimilarities to those seen in intraannual spatial turnover and found the spatial distance at which there were comparable dissimilarities. We used 2015 as our baseline for spatial structure because many more samples were collected over a wider range of distances in 2015 than in subsequent years. Spatial distances and dissimilarities were derived from the linear models fitted to the first 10 cm of distance decay curves described in [33].

OTU Reoccurrence by Relative Abundance
To assess the effect of OTU relative abundance in 1 year on likelihood of that OTU being found again in a subsequent year in the same location, we divided relative abundance values into five levels: > 0 to 0.01%, > 0.01 to 0.1%, > 0.1 to 1%, > 1 to 10%, and > 10%. For each relative abundance category, we calculated the percentage of OTUs that were present 1 and 2 years later at the same sampling location and substrate. This process was repeated for each substrate across all three pairs of years, 2015 to 2016, 2016 to 2017, and 2015 to 2017. To account for differences in sequencing depth among samples, we repeated this analysis with data rarefied to 1000 reads per sample. Rarefactions were performed 1000 times, and the resulting percentages from each rarefaction were averaged. The lowest relative abundance class, > 0 to 0.01%, was excluded from the analysis with rarefied data.

Effect of Sampling Year
In the RDA and variance partitioning analyses, substrate type, sampling year, and PCNM spatial vectors all had significant effects (p < 0.001) on community composition for both presence-absence and relative abundance data. All predictors combined explained 8.2% and 10% of the variation in community composition for presence-absence and relative abundance, respectively. In both cases, year alone explained the least variation, at less than 1% (Fig. 1, S2).
NMDS ordination (Fig. 2) showed that the same patterns in fungal community composition initially seen in 2015 persisted over two subsequent years. When all years and substrates are plotted together, samples cluster by substrate type rather than by sampling year. When the data are separated by substrate and each substrate is analyzed separately, samples do not cluster by year (Fig. 3), suggesting there is no strong effect of sampling year on community composition. There were also significant

Alpha Diversity
We observed 11,410 OTUs across all substrates and sampling years. For both OTU richness and the Shannon index, sampling year (χ 2 = 108.1 and 68.5 respectively, p < 0.001 for both), substrate type (χ 2 = 52.4 and 38.8, p < 0.001), and their interaction (χ 2 = 35.0 and 13.9, p < 0.001 and p = 0.0076) all had significant effects on diversity. OTU richness and Shannon diversity increased from 2015 to 2016 (p < 0.001) but remained steady from 2016 to 2017 (Fig. S4). Live bryophytes had greater OTU richness than other substrates in 2015 (p < 0.02) and greater richness than bark in 2016 and 2017. Communities in live bryophytes also had greater Shannon diversity than other substrates in 2016 and 2017 (p < 0.05), but all substrates were same in 2015. OTU richness in outer host tree bark changed the least between years. Overall and per sample OTU richness were high and most OTUs were rare: 58.8% of OTUs were observed three or fewer times over the course of the study.
Changes in fungal community composition over the course of 1 year, as quantified by the Bray-Curtis index, were comparable to variation seen over 5.7 to 8.3 cm within a single year (Table S1, Fig. S5), distances at which distance decay is apparent and fungal communities are more similar to each other than expected by chance. While substrates differed in their average temporal turnover, when this inter-annual change is compared to variation with space, the spatial distance needed to reach the same level of community turnover seen over 1 year is roughly the same across substrates. Fungal communities in live bryophyte samples had the highest average temporal turnover, and they also had the highest levels of intra-annual spatial turnover and dissimilarity among samples (Fig. S5, Table S2). Temporal turnover using presence-absence and the Jaccard index was more variable when compared to space. In one case, with live bryophyte fungal community change from 2016 to 2017, inter-annual dissimilarity was lower than the y-intercept value for the spatial distance decay curve; that is, there was less dissimilarity among years than would be expected at a distance of 0 cm. Spatial distances comparable to the temporal dissimilarities varied from − 2 to 79 cm (Table S1).

Relative Abundance and Persistence
OTUs at a higher relative abundance persisted to the next year more often than those at lower relative abundances in both raw data and data rarefied to 1000 reads per sample (Fig. 6, S6 Table S3). OTUs at most relative abundance levels were most likely to persist in the dead bryophyte substrate and least likely to persist in live bryophytes. OTU persistence rates were higher for 2015-2016 than in other pairs of years. For most substrate/year combinations, an OTU needed to make up at least 1% of the fungal reads in a sample to be present again the next year at least half of the time. When present, OTUs typically had very low relative abundances and seldom reached this 1% threshold (Fig.   Fig. 4 Boxplot of paired Bray-Curtis dissimilarities by substrate and year pair. Each data point is the dissimilarity between samples at the same sampling location in different years. Live bryophyte samples were the most dissimilar over time Fig. 5 Boxplot of paired partitioned Jaccard dissimilarities by substrate and year pair. Each data point is the dissimilarity between samples at the same sampling location in different years. Nestedness was low, while turnover and total Jaccard dissimilarity were high S7a). Individual OTUs made up 1% or more of a sample only 5.5% of the time and reached 10% of a sample 0.69% of the time. The most abundant OTU in a given sample, however, usually made up a large portion of said sample, with a mean relative abundance of 26% (Fig. S7b). These most abundant OTUs reached or exceeded 10% of reads 95.6% of the time.

Discussion
While fungi present at individual sampling locations varied among years, the effect of sampling year on the overall fungal community was modest and smaller than the effects of space or substrate type. Patterns of fungal specificity among substrates observed in a single year persisted over the duration of the study (Fig. 2), and fungi that were widely distributed in 1 year tended to also be common in other years (Fig. S3). Year explained less than 1% of community variation (Fig. 1), so while the effect was statistically significant, its biological significance is dubious. These findings indicate that the fungal community in this epiphytic environment was largely stable at the broad, site-level scale. The same set of taxa are present from year to year at similar levels of abundance and within the same substrates, though their exact locations within a site may shift. This is consistent with findings in soil systems that inter-annual changes tend to be small compared to intra-annual seasonal ones [26,31].
In contrast to community composition, alpha diversity did change among years. There was an increase in per sample OTU richness from 2015 to 2016, but it remained stable from 2016 to 2017 (Fig. S4). The reason for this increase is unknown. It could have been related to an environmental change or a change in our sample preservation or PCR cleanup methods. Longer term monitoring is needed to better understand the dynamics and drivers of fungal biodiversity in the epiphytic environment.
If fungal communities varied randomly from year to year, we would expect temporal dissimilarities to match those seen at larger spatial distances, in which there is no detectable spatial structure and distance decay curves have reached an asymptote (Fig. S5, Table S2, generally > 1 m in this system [33]). Instead, inter-annual dissimilarities matched those seen at much shorter distances, indicating that communities within individual sampling locations were more similar to each other over time than expected due to chance, and that there is continuity of fungi over the course of multiple years. When Jaccard dissimilarities were partitioned into turnover and nestedness components, turnover predominated (Fig. 5). This indicates there was not just a simple gain or loss of OTUs and fungal diversity with time, but that both were happening simultaneously, resulting in turnover of community members.
Inter-annual turnover within sample locations, as measured by Bray-Curtis and Jaccard dissimilarity varied among substrates (Figs. 4,5). When compared to spatial turnover in the same substrate, however, the spatial distances needed to see the same amount of Bray-Curtis dissimilarity were remarkably similar among substrates, ranging from 5.7 to 8.3 cm (Table S1). For example, while live bryophytes had the greatest OTU turnover with time, they also had the highest levels of change across space. Live bryophytes may have the highest inter-annual dissimilarity and OTU turnover because living bryophyte tissues are actively growing, creating new substrate for fungi to colonize, unlike the underlying dead bryophytes or tree bark. The growing bryophytes may also be producing anti-fungal compounds [47,48] which could limit the growth rate or longevity of fungi. Fig. 6 OTU reoccurrence rates at different initial relative abundances by substrate and year pair. OTUs with a higher relative abundance in the first year were more likely to be present in the second year. OTUs in live bryophytes tended to reoccur less often than in other substrates Live bryophytes, when present, were also the top-most substrate layer, leaving them most exposed to incoming fungal spores and disturbance, both of which could increase fungal community turnover. Dead bryophytes and host tree bark are less exposed and, we expect, more physically and chemically stable, allowing for greater fungal community stability.
Individual fungal OTUs were more likely to persist to the next year if they were at high relative abundance in the first year (Fig. 6, S6, Table S3). The least abundant OTUs making up 0.1% or less of a sample seldom persisted to the next year. In contrast, the most abundant OTUs usually did persist. This suggests that fungi able to establish themselves and accumulate biomass can maintain their presence for at least a year. Greater persistence in higher abundance fungi also implies stability in ecological functions. Fungi that are more abundant and have greater biomass presumably have more functional significance than low abundance ones represented by only a few reads and that might only be inactive spores. The larger ecological impacts of these more abundant fungi, possibly including decomposition and nutrient cycling, inhibiting plant growth via disease, or promoting plant growth through mycorrhizal symbiosis, can persist and remain stable over multiple years. This temporal stability of high abundance fungi coupled with the high spatial turnover seen in this system [33] could create small but persistent patches with greater or lesser suitability for germination and growth of epiphytes, for example.
These dominant fungi likely drive the patterns of temporal stability in this system. The spatial distance needed to see the amounts of Jaccard dissimilarity observed among years was typically much higher than for Bray-Curtis. The Jaccard metric only considers presence-absence and is blind to relative abundance. All fungi, regardless of relative abundance, are given the same weight, and the many low abundance, but non-persistent, taxa create high dissimilarity values. With Bray-Curtis, the more abundant taxa, which tend to be more persistent, drive down the dissimilarity values.
The dominance of certain fungi could be due to their early arrival at that location compared to other fungi [49,50]. Early colonizing microbes, however, do not always persist long term, as has been shown in decaying logs [51] and biofilms [52]. The host tree branches used in this study were also several years old, as we avoided sampling branch tips or twigs, making it unlikely that we were observing early stages of succession. Wide variation in relative abundances could also be due to environmental factors, with particular fungi being better suited to different environments, as seen with fungi assorting among different substrate types (Fig. 2). High and low relative abundances may also be a function of the traits and life histories of different fungi, with some being short-lived ruderal species while others are larger with greater longevity. It is unclear how much of the gain and loss of OTUs is due to life history strategies, particularly short-lived r-strategists which might live less than 1 year even in favorable circumstances, versus failure of fungi to establish in the first place, regardless of their life history. Sampling at shorter time intervals would allow for greater insights into the dynamics of these diverse and everchanging communities.
The number of dominant fungi was low, with the overwhelming majority of OTU observations being less than 1% of a sample. Most fungi were also rare, with the majority of OTUs being observed in three or fewer samples over the course of the study. This, combined with high spatial turnover, implies that most fungi in the epiphytic environment are short lived, due to either an unsuitable environment or their own biology, and likely rely on stochastic spore dispersal to spread within and between branches, rather than continuous mycelial growth. Despite this high temporal and spatial turnover, with many fungi appearing and disappearing erratically, there is an underlying stability, at both the site and sample scale, of common and high abundance fungi that are able to persist. These taxa presumably have sustained functional impact on their environment, but further research using techniques like transcriptomics or enzyme assays are needed to better understand temporal stability of fungal activity and function.