Soil Bacterial and Fungal Communities Show Distinct Recovery Patterns during Forest Ecosystem Restoration

The processing of raw sequencing reads resulted in 16.8 million high-quality V4 sequences, averaging about 280,000 sequences per sample, and 5.1 million high-quality internal transcribed spacer 1 (ITS1) sequences, averaging about 86,000 sequences per sample. The rarefaction curves are shown in Fig. S2 in the supplemental material. Of the initial 60 samples for ITS1 sequencing, 59 samples were used for downstream analyses. Sequences were then clustered into 17,017 bacterial and 5,979 fungal operational taxonomic units (OTUs), with a 97% identity threshold (35). For all downstream analyses, 97,000 reads were randomly selected for bacteria, and 28,000 reads were randomly selected for fungi per sample to correct for differences in sequencing depth.

The number of bacterial 16S rRNA copies decreased significantly from May and July to September and October, while the number of fungal ITS copies increased significantly over the same period (P < 0.05, Tukey's honestly significant difference [HSD] test) (Fig. 1e and g). However, along chronosequence ages, there were no obvious patterns in the 16S rRNA or ITS copy numbers, except that bacterial abundance was significantly lower at 5 years (P < 0.05, Tukey's HSD test) (Fig. 1f and h). Overall, the monthly shifts in bacterial and fungal abundance were equal to or greater than the shifts across chronosequence ages. Bacterial richness, estimated by the Chao1 index, was significantly higher in May and July than in September (but not October), while fungal community richness was significantly higher in May than in September only (P < 0.05, Tukey's HSD test) (Fig. 1a and c). Both bacterial and fungal richness increased generally with increasing chronosequence age up to 30 years, but they were lower in the unmined sites (Fig. 1b and d). There were no significant interactions between chronosequence age and sampling month for either richness or abundance of either bacteria or fungi.

The patterns of bacterial and fungal β-diversity were visualized with nonmetric multidimensional scaling (NMDS) plots and tested for significant differences using analysis of similarity (ANOSIM). Bacterial β-diversity varied significantly within the study year, with May and July harboring significantly different communities from those of September and October, regardless of chronosequence age (P = 0.001, ANOSIM with 999 permutations) (Fig. 2a and Table S1a). The differences between May and July or between September and October were not significant (P = 0.437 to 0.846, ANOSIM with 999 permutations) (Table S1a). All pairwise comparisons of bacterial community structure among chronosequence ages were significantly different (P = 0.001 to 0.037, ANOSIM with 999 permutations) (Table S1b). Interestingly, the overall pattern of β-diversity across the chronosequence was maintained in each of the two clusters (Fig. 2a), with bacterial community structure becoming more similar to those of reference sites with increasing chronosequence age. The magnitude of the shift in bacterial β-diversity within a year was comparable to that observed across decades of succession, as evidenced by separation on the NMDS plot (Fig. 2a) and average Bray-Curtis similarity values among samples (data not shown). In contrast, fungal communities from different chronosequence ages clustered more tightly, without any overlap across ages (Fig. 2b), and communities in each age were significantly different from all other ages (P = 0.001, ANOSIM with 999 permutations). In contrast, intra-annual variability across months in each chronosequence age was not significant (P = 0.146 to 0.996, ANOSIM with 999 permutations) (Table S2). In addition, multivariate dispersion analysis showed that the dispersion of bacterial communities by month and chronosequence age was not significantly different (Welch's two-sample t test, t = 1.629, P = 0.147), while for fungi, the dispersion of fungal communities in each age class was significantly lower (i.e., tighter clustering) than when tested by month (Welch's two-sample t test, t = 3.2691, P = 0.023).

Patterns of β-diversity were also related to soil and environmental data collected by Avera et al. (36) at each site. Changes in temperature, moisture, total microbial biomass (TMB), substrate-induced respiration (SIR), and NH₄⁺ levels were significantly correlated with bacterial community shifts (P = 0.001 to 0.041) (Table S3). Moisture and NH₄⁺ levels were mainly correlated with the intra-annual bacterial community shifts, with higher moisture in May and July and higher NH₄⁺ in September and October (Fig. 2a). In contrast, SIR was mainly correlated with bacterial community shifts that occurred with chronosequence age, with an increase along the chronosequence. For fungi, changes in temperature, moisture, TMB, SIR, and NO₃⁻ levels were significantly correlated with community shifts (P = 0.001 to 0.004) (Table S3). SIR, TMB, and temperature were most strongly correlated with fungal shifts in β-diversity (Fig. 2b), with higher SIR and TMB and lower soil temperature associated with fungal communities from older chronosequence ages than from newer ones. It is important to remember, however, that associated environmental factors could also be autocorrelated, and additional experimentation would be required to determine which might be directly driving or resulting from changes in soil microbiota.

The dominant bacterial phyla and proteobacterial classes across all samples were Acidobacteria and Alphaproteobacteria, respectively, representing 21.5% and 18.0% of all sequences assigned to the domain Bacteria (Fig. 3a). Sampling month had significant influences on some taxa, with the relative abundances of Bacteroidetes, Deltaproteobacteria, Gemmatimonadetes, Betaproteobacteria, and Gammaproteobacteria significantly higher in May and July and the relative abundances of Chloroflexi and Planctomycetes significantly higher in September and October (P < 0.05, Tukey's HSD test) (Fig. 3b and Table S4). Chronosequence age also had significant influence on some taxa (P < 0.05, analysis of variance [ANOVA]) (Fig. 3c and Table S4). The relative abundances of Actinobacteria, Chloroflexi, and Gemmatimonadetes were significantly higher in the 5-year sites than other sites, while the relative abundances of Nitrospirae, Alphaproteobacteria, and Verrucomicrobia were significantly higher in later chronosequence ages. Nitrospirae were higher in the 30-year and unmined (UM) sites than in the 5-year and 11-year sites, Alphaproteobacteria were higher in the UM sites than in other sites, and Verrucomicrobia were higher in UM sites than in the 5-, 11-, and 21-year sites (P < 0.05, Tukey's HSD test). A summary of the most abundant OTUs that changed significantly across the chronosequence ages is presented in Fig. S6A.

Ascomycota and Basidiomycota were the dominant fungal phyla at all sites, representing 46% and 30% of classified OTUs, respectively, but their relative abundances were not clearly related to chronosequence age or sampling month. The dominant classes across all samples were Agaricomycetes (phylum Basidiomycota; 28.7%) and Sordariomycetes (phylum Ascomycota; 20.0%) (Fig. 3d). Sampling month had significant effects on some fungal classes (P < 0.05, ANOVA) (Fig. 3e and Table S4). Dothideomycetes was significantly different between July and October, and incertae sedis within Zygomycota was significantly different between May and September. Chronosequence ages had significant effects on more taxa, including Dothideomycetes, Eurotiomycetes, Leotiomycetes, Sordariomycetes, Rozellomycota, and incertae sedis within Zygomycota (P < 0.05, ANOVA) (Fig. 3f and Table S4). The relative abundances of Dothideomycetes and Sordariomycetes were significantly higher in the 5-year sites than in other sites, while the relative abundance of Agaricomycetes was significantly higher in the UM sites than in the 5-year sites. The relative abundance of incertae sedis in Zygomycota was significantly higher in the 30-year and UM sites than in the 5- and 11-year sites. All related P values are shown in Table S4, and the most abundant fungal OTUs that changed significantly across the chronosequence ages are summarized in Fig. S6B.

Cooccurrence networks were generated separately for bacteria and fungi at each chronosequence age, and topological properties were calculated to characterize changes over time. The numbers of nodes and edges in the bacterial networks typically increased with increasing chronosequence age, with the unmined reference network having the highest number of nodes and edges (Table 1). Moreover, the average degree also increased from early chronosequence ages (5 and 11 years) to later chronosequence ages (21 and 30 years) (Table 1 and Fig. 4). The average degree of the unmined reference bacterial network (at 11.66) is about 1.5 times those at 21 and 30 years (average degrees, 7.53 and 7.16, respectively). The networks in fungal communities, however, were markedly different from those in bacterial communities (Table 1 and Fig. 5). The numbers of fungal nodes and edges both peak at 21 years (Table 1), with the fungal networks at 5 years having the smallest number of nodes and edges and the lowest average degree. The average degree of the unmined reference network is larger than that at 5 years but smaller than those at 11 years and 21 years. Overall, the nodes, edges, and average degree of fungal networks typically peaked in the middle chronosequence ages.

Networks in consecutive chronosequence ages were compared to identify overlapping nodes and edges. On average, 76% of bacterial nodes and 8% of bacterial edges in an earlier age appeared in the next age, and overlaps of both bacterial nodes and edges are significantly higher than would be expected from random (Fisher's exact test for node, P = 0.001 to 0.004; permutation test for edge with 999 permutations, P = 0.001) (Table S5). In contrast, a mean of 27% of fungal nodes and 0.5% of fungal edges in an earlier age appeared in the next chronosequence age, and the numbers of overlapping fungal nodes and edges are not significantly greater than random (Fisher's exact test for node, P = 0.872 to 0.999; permutation test for edge with 999 permutations, P = 0.614 to 0.990) (Table S5).

Among the overlapping bacterial nodes, about 30% belong to the phylum Acidobacteria, which is higher than the relative abundance of Acidobacteria among all OTUs used for network analysis (22%). Bacteroidetes and Planctomycetes were also overrepresented in overlapping nodes compared to their percentages in all OTUs, at 14% and 10% compared to 10% and 6%, respectively. In contrast, taxa that were underrepresented in overlaps included Actinobacteria (5% of overlapping nodes versus 13% of OTUs) and Alphaproteobacteria (11% of overlapping nodes versus 16% of OTUs). Among the 352 overlapping bacterial edges, 223 edges (63%) involve Acidobacteria, and 146 of 352 edges (42%) are between two Acidobacteria nodes. More broadly, about 86% of overlapping edges involve Acidobacteria, Bacteroidetes, or Planctomycetes, even though the combined relative abundance of these 3 taxa is only 38% of the total OTUs used for network analysis. Nodes involved in the overlapping edges have significantly higher degrees and closeness centrality than average (P = 0.001 to 0.006, Welch's two-sample t test), but the relative abundances of OTUs represented by those nodes are not always significantly different from average (P = 0.001 to 0.683, Welch's two-sample t test) (Table S6). The overlapping fungal nodes mainly belong to Agaricomycetes, Sordariomycetes, and Dothideomycetes, with Agaricomycetes being overrepresented (15% of overlapping nodes versus 12% of OTUs). In contrast to bacteria, there were only 3 overlapping fungal edges among all possible pairs of consecutive chronosequence ages.

The bacterial and fungal networks show different patterns of cooccurrence, typically with higher numbers of nodes and edges in bacterial networks. The bacterial networks typically contained a much higher number of interactions dominated by one large connected component that consisted of more than 65% of all OTUs (Fig. 4). In contrast, fungal networks consisted of much fewer interactions scattered across multiple small clusters (i.e., only a few OTUs) that were discrete and densely clustered (higher clustering coefficients) (Fig. 5 and Table 1). It is important to note that the relatively low intra-annual turnover rates of fungi may result in less-detectable cooccurrences. However, the stark differences between the two sets of networks, which were observed consistently even when the networks were reanalyzed with a range of alternative correlation cutoffs (Fig. S3), suggest that there are real and significant differences between these sets of interactions. We also constructed bacterium-fungus cooccurrence networks with the 500 most abundant bacterial OTUs and the 500 most abundant fungal OTUs. In these networks, there are more bacterial nodes and edges than fungal (Fig. S5 and Table S8), which is consistent with our observations in the separate networks. Also, the number of bacterium-bacterium edges increased with increasing chronosequence age, while the number of fungus-fungus edges peaks at 21 years; these findings are also consistent with the trends in the separate networks. The number of bacterium-fungus edges increased with increasing chronosequence age, but they comprised only 7% of all edges on average. Additional characteristics of bacterial and fungal networks, including scale-free indices, small world indices, and modularity, are summarized in Table S7 and Fig. S4.

The richness of both bacteria and fungi peaked at 30 years and was lower in the unmined reference sites. We hypothesize that the short time since disturbance at the younger sites was limiting, such that richness was still increasing up to the 30-year sites. This aspect of ecosystem recovery allows the coexistence of both r- and k-selected species at the older sites, which could explain the peak richness at 30-year site (37, 38). This general trend is in agreement with other observations of microorganisms in aquatic (39) and soil (40, 41) systems.

The β-diversity also followed a clear pattern for both bacteria and fungi, in that the community structure became more similar to the unmined sites as chronosequence age increased. The bacterial community shifts were more gradual and overlapping, while the fungal communities at different ages were more distinct. This agrees with previous observations that plant symbioses with fungi are more common than with bacteria (16) and that plant community composition has strong influences on the fungal community (17–20), which has been shown to extend to bulk soil in the form of common mycelial networks (42, 43). The specificity with vegetation may explain the distinctness of fungal communities at different chronosequence ages, as the plant community is also experiencing taxonomic shifts from being grass dominated to tree dominated during this time frame (44). Furthermore, as dominant decomposers of litter (45), the varied litter properties with different vegetation (46, 47) might also influence fungal community composition. In contrast, most bacteria inhabit small-scale niches in bulk soil and have a less-direct connection with plants than fungi (48). However, it is important to remember that in addition to external environmental factors, the fact that community changes along the chronosequence age were correlated with TMB and SIR suggests that autogenic factors also play important roles in microbial recovery.

The patterns of intra-annual variability among bacterial and fungal communities were also markedly different (Fig. 1 and 2), with bacterial abundance higher in May and July and fungal abundance higher in September and October. This shift was likely related to seasonal dynamics of plant inputs, temperature, and precipitation. For example, fungi are often identified as the primary decomposers of recalcitrant organic matter (45), such as cellulose and lignin, which are increased in the later seasons due to plant litter and may favor higher fungal abundance in the later sampling months. Based on the ANOSIM and NMDS visualization, bacterial community structure shifted significantly among sampling months, while fungal community structures did not. This difference in intra-annual variability may also be related to growth rates and highlights the ecological differences between these two communities. For example, soil fungal biomass turnover times have been estimated to be >100 days, which is about 10-fold slower than those of typical soil bacteria (14, 49). In addition, bacteria are more sensitive to changes in soil moisture (50, 51), which correlated significantly with the intra-annual bacterial community shifts observed in this study (r² = 0.157, P = 0.012). Short-term variability of bacterial communities has been observed previously, but the magnitude of change varies with soil or ecosystem type (52–54). In our study, the β-diversity indicates that the shifts of soil bacterial communities within a 6-month span could be as large as those observed across decades of ecosystem recovery, although the effect of chronosequence age was still observable within each monthly shift (Fig. 2a). In contrast, fungal community structures did not change significantly over 6 months, implying that the majority of the change in soil fungal communities occurs over longer time scales than for bacteria.

Changes in the types of soil microbiota across chronosequence age were apparent at both low and high taxonomic levels, including bacterial phylum and fungal class (Fig. 3c and f), which suggests that significant taxonomic changes were occurring across the entire community and supports the theory that at least some degree of ecological coherence exists among broad phylogenetic groups (55–57). In succession, bacterial communities appeared to shift from copiotrophic groups to oligotrophic groups, with decreased Bacteroidetes and increased Verrucomicrobia (Fig. 3c) (58–61). This shift potentially resulted from the decrease in abundant labile C and N substrates and the increase in nutrient-poor and recalcitrant C compounds, which is supported by a measured increase in mineral soil C:N ratio with increasing chronosequence age at these sites (36). Ecosystem changes are also apparent at lower taxonomic levels. Many of the OTUs that increase significantly over the chronosequence, such as Bradyrhizobium, Chthoniobacteraceae (genus DA101), Mycobacterium, and Rhodoplanes, indicate increased soil organic matter decomposition and nutrient cycling as conditions become more oligotrophic. Interestingly, we also saw an increase in “Candidatus Xiphinematobacter” in our data, suggesting an increase in nematode populations (62) and the potential development of higher trophic levels in the soil ecosystem.

Overall, these composition changes also likely result from broader ecosystem succession and vegetation change, as plant inputs (leaf litter, root biomass, and exudates) vary with vegetation type (46, 47, 63) and have a strong influence on microbial community composition (64, 65). For example, we observed increased Alphaproteobacteria and decreased Actinobacteria (Fig. 3c), which is a change consistent with a shift from grassland to forest soil (66). With regard to fungi, many of the Dothideomycetes are tolerant of extreme environmental conditions, including extreme heat and cold, drought, and UV radiation, and some produce enzymes that help degrade rocks (67, 68). Given that the 5-year sites lacked mature tree coverage that can help buffer rapidly fluctuating conditions and reduce direct UV radiation, Dothideomycetes may have competitive advantages in this relatively harsh environment (Fig. 3f). Among more specific taxa, we saw increases in fungal genera, like Hygrophorus, which includes tree-associated ectomycorrhizae (69, 70), and other common forest saprotrophs, such as Entoloma and Mortierella.

In larger organisms, the complexity of ecological networks has been shown to increase during ecosystem succession, leading to greater food web stability (71). Similarly, we observed increased bacterial cooccurrence network complexity with increased chronosequence age, with more cooccurring OTUs retained and higher connectivity in later ages. Functionally, this development was also accompanied by previously reported increases in microbial activity at these sites (36). Zelezniak et al. (72) found that, across multiple habitat types, species cooccurrence tends to be mainly driven by metabolic dependencies. This supports the connection between network complexity and microbial activity at these sites and illustrates the need to conduct more complex investigations of how interaction networks can be used to gain insight into the relationship between changes in microbial function and diversity and how they relate to ecosystem change.

To further explore changes in cooccurrence networks at different chronosequence ages, we quantified the overlaps of networks at consecutive stages, including the significance of overlapping OTUs and their taxonomic compositions. The significantly higher-than-average numbers of overlaps in both bacterial OTUs and cooccurrences indicate that a portion of the bacterial interaction network persists over relatively long periods of time during the recovery of soil microbiota. In contrast, fungal networks changed more substantially and displayed new interactions in each chronosequence age, with no more overlaps between stages than would be expected from random chance, suggesting minimal to no influence from the previous stage. These network differences are consistent with the β-diversity results, illustrating a gradual shift in bacterial communities between ages but relatively distinct shifts between fungal communities.

In a simulation of cooccurrence networks, Berry and Widder (73) identified keystone species as having higher values for mean degree and closeness centrality. Applying this rationale, we identified bacterial OTUs in overlapping edges as potential keystone species which have degrees and closeness centrality that are significantly higher than average (Table S6). The potential keystone OTUs observed in our study were dominated mostly by Acidobacteria but also by Bacteroidetes and Planctomycetes. It is important to note that these are not simply a result of more abundant taxa having a higher chance to cooccur with other taxa, because the relative abundances of some keystone OTUs identified in our study were not significantly different from average (Table S6). The identification of potentially important OTUs is critical to improving our understanding of the ecological roles played by specific bacterial taxa. For example, although Acidobacteria is one of the most abundant phyla in soils and represented many of the potential keystone OTUs identified in this study, little is known about these organisms, because they have been very difficult to isolate and culture (74, 75).

Overall, these results provide a comprehensive view of how the recovery of the bacterial and fungal components of the soil microbiota can be strikingly different during ecosystem recovery following extreme disturbance. This work highlights numerous ways in which different properties of soil bacteria and fungi can influence their variability during ecosystem recovery within the same environment. Further investigations that focus on the ecological roles of specific taxa, microbial interactions, and potential keystone species are essential to obtaining a better understanding of the mechanisms that shape patterns of microbial diversity and function during ecosystem recovery after disturbances.

The reforested mined sites are located in Wise County, VA, USA, in the central Appalachian Coal Basin (76). The samples for this study were collected in collaboration with work done by Avera et al. (36) to look at ecosystem recovery at the same sites. Details of the sampling sites (climate, aspect, slope, vegetation, etc.) and plot design are described in the “Soil dynamics” section of that publication. Briefly, four reclaimed sites were selected to establish a chronosequence of 5, 11, 21, and 30 years postreforestation, and a nearby unmined site in a natural hardwood forest was used as a reference (Fig. S1). Triplicate study plots were randomly located within each site. Three soil samples were taken from 0 to 10 cm depth in each plot, pooled, and homogenized to more fully capture the diversity at the plot level, while using the replicated plots to characterize variation in the soil microbiota at each site. Soil samples were collected in May, July, September, and October 2013 to estimate the intra-annual variation of microbial communities. All samples were transported to the lab on ice, sieved to 4 mm to remove coarse roots and stones, and stored at −80°C until DNA extraction. As part of their research, Avera et al. (36) measured other ecosystem properties, and some of their results are further analyzed here to aid in an understanding of the factors influencing microbiota recovery. These include soil properties, such as temperature (Temp), moisture, and available NH₄⁺ and NO₃⁻; gas fluxes of N₂O, CO₂, and CH₄; and total microbial biomass (TMB) and substrate-induced respiration (SIR).

Total DNA was extracted from approximately 0.25 g of soil using the Mo Bio PowerSoil DNA isolation kit (Mo Bio Laboratories, Carlsbad, CA, USA). Quantitative PCR (qPCR) was performed to estimate the abundances of bacteria and fungi in soil samples. Bacteria were targeted with the Eub338 (5′-ACT CCT ACG GGA GGC AGC AG-3′) and Eub518 (5′-ATT ACC GCG GCT GCT GG-3′) primer pair, and fungi were targeted with ITS1f (5′-TCC GTA GGT GAA CCT GCG G-3′) and 5.8s (5′-CGC TGC GTT CTT CAT CG-3′) primer pair (77). Each 25-μl reaction mixture contained 12.5 μl of iTaq Universal SYBR green supermix (Bio-Rad, Hercules, CA, USA), 5 μl of PCR water (Mo Bio Laboratories), 1.25 μl of each of the primers (10 μM), and 5 μl of isolated DNA as the template. The conditions were 15 min at 95°C, followed by 40 cycles of 95°C for 1 min, 53°C for 30 s, and 72°C for 1 min. Standards were made from a 10-fold dilution series of plasmids containing cloned target sequences. Standard curves for all qPCR runs had efficiency values ranging from 0.97 to 1.07 and R² values of >0.99.

The V4 region of the bacterial 16S rRNA gene and the fungal ITS1 region were amplified by PCR to characterize bacterial and fungal communities, respectively. The PCR mixture (25 μl) contained 10 μl of 2.5× 5 Prime HotMaster mix (5 Prime, Gaithersburg, MD), 13 μl of PCR water (Mo Bio Laboratories, Carlsbad, CA, USA), 0.5 μl of each of the primers (10 μM), and 1 μl of isolated DNA as the template. The V4 region was amplified with the 515F and 806R primer pair (78), and the ITS1 region was amplified with the ITS1F and ITS2 primer pair (79). Primers were designed for Illumina sequencing with adapters, primer pads, and 2-bp linker sequences, and the reverse primer contained a 12-bp barcode sequence unique to each sample (78, 79). The thermal cycling protocol was 94°C for 5 min, 35 cycles of 94°C for 45 s, 50°C for 45 s, and 72°C for 90 s, followed by a final extension period at 72°C for 10 min. All samples were amplified in triplicate on thermal cyclers (Bio-Rad, Hercules, CA, USA), pooled, and visualized on agarose gels. PCR products were purified using the UltraClean PCR cleanup kit (Mo Bio Laboratories). Amplicon concentrations were quantified using a Qubit 2.0 fluorometer (Invitrogen, USA), according to the manufacturer's protocol, and amplicons from all samples were combined in equimolar ratios. Bacterial and fungal amplicons were sequenced separately on the Illumina MiSeq platform using 250 bp and paired ends. Sequencing reads were quality filtered and processed using the USEARCH pipeline (80). In brief, the forward and reverse reads were merged with at least 150 bp overlap. The merged reads were filtered with a minimum length of 200 bp and expected errors per read lower than 0.5, which are calculated from quality scores. Chimeric sequences were identified with UCHIME and removed (81). The SILVA and UNITE databases were used as references for bacteria and fungi, respectively (82, 83), and taxonomy was assigned using the RDP Classifier (84). The Chao1 index was generated with Quantitative Insights Into Microbial Ecology (QIIME) to estimate the richness of bacterial and fungal communities (85).

The R program (86) was primarily used for statistical analyses and plotting. Significant differences in bacterial and fungal abundance and richness were detected across chronosequence ages or sampling months using ANOVA and Tukey's HSD test. Matrices of pairwise similarities for all samples were calculated independently for bacterial and fungal communities using the Bray-Curtis distance metric. Patterns of similarity among samples were visualized with nonmetric multidimensional scaling (NMDS) using the metaMDS function in the R vegan package (87), and environmental variables (Temp, moisture, and levels of NH₄⁺, NO₃⁻, N₂O, CO₂, and CH₄), TMB, and SIR (36) were added to the ordination with the envfit function in the same package. The ANOSIM function in vegan was used to test the statistical significance of all pairs of samples, with 999 permutations. The multivariate dispersion of all months and chronosequence groups was calculated with function betadisper. Microbial OTUs that showed significant differences in abundance across sampling months or chronosequence ages were selected using the edgeR program, with a false-discovery rate (FDR) of <0.05 and count per million (CPM) of >2⁶ in order to choose highly abundant and significantly variable OTUs for downstream analysis (88). After filtering, 2,480 bacterial OTU and 1,340 fungal OTU remained for further analysis.

Cooccurrence analyses have recently been increasingly used in microbial ecology to better understand community structure, characterize potential intracommunity interactions, and identify potentially shared and contrasting niches (89, 90). The cooccurrence network analysis in this study was performed with the igraph R package (91). The 500 most abundant OTUs (accounting for 61.6% of bacterial sequences and 48.0% of fungal sequences) per chronosequence age were used to build individual networks for each age. The cutoff value of 500 was based upon a similar approach used by Dini-Andreote et al. (23), but we also constructed networks using the most abundant 800 and 1,000 OTU to verify that the interpretation of the general trends of network properties did not change. Spearman rank correlations (SRC) were calculated for each pair of OTUs in the same chronosequence age. The edges in the networks represent statistically significant (FDR < 0.05) correlations, with an absolute SRC value of >0.9. The nodes in the networks represent individual OTUs, with at least one significant correlation with another OTU. The numbers of nodes and edges, average degree, clustering coefficient, and closeness centrality were calculated using the igraph R package. Degree is the number of edges of each node, representing how many other nodes (OTUs) in the network are connected with the given node. The clustering coefficient is a measure of how often a node's neighbors are also connected (92). Closeness centrality measures the inverse distance to all other nodes from one node (93). The overlapping nodes and edges of networks from two consecutive stages were identified, and the significance of the difference between the number of actual overlaps and randomized overlaps were analyzed with Fisher's exact test and permutation methods (999 permutations). The abundance, degree, and closeness centrality of nodes in overlapping edges were compared with the average of the whole network with Welch's two-sample t tests.

The resulting raw reads were deposited in the NCBI BioProject database (accession number PRJNA324696 ).