Emily C. Farrer1,2*
Dorota L. Porazinska3
Marko J. Spasojevic2,4Andrew J. King5,6
Clifton P. Bueno de Mesquita2,7Samuel A. Sartwell2,7Jane G. Smith2
Caitlin T. White2,7
Steven K. Schmidt7Katharine N. Suding2,7- 1Department of Ecology and Evolutionary Biology, Tulane University, New Orleans, LA, United States
- 2Institute of Arctic and Alpine Research, University of Colorado Boulder, Boulder, CO, United States
- 3Department of Entomology and Nematology, University of Florida, Gainesville, FL, United States
- 4Department of Evolution, Ecology, and Organismal Biology, University of California, Riverside, Riverside, CA, United States
- 5Biosciences Division, Oak Ridge National Laboratory, Oak Ridge, TN, United States
- 6King Ecological Consulting, Knoxville, TN, United States
- 7Department of Ecology and Evolutionary Biology, University of Colorado Boulder, Boulder, CO, United States
While it is well established that microbial composition and diversity shift along environmental gradients, how interactions among microbes change is poorly understood. Here, we tested how community structure and species interactions among diverse groups of soil microbes (bacteria, fungi, non-fungal eukaryotes) change across a fundamental ecological gradient, succession. Our study system is a high-elevation alpine ecosystem that exhibits variability in successional stage due to topography and harsh environmental conditions. We used hierarchical Bayesian joint distribution modeling to remove the influence of environmental covariates on species distributions and generated interaction networks using the residual species-to-species variance-covariance matrix. We hypothesized that as ecological succession proceeds, diversity will increase, species composition will change, and soil microbial networks will become more complex. As expected, we found that diversity of most taxonomic groups increased over succession, and species composition changed considerably. Interestingly, and contrary to our hypothesis, interaction networks became less complex over succession (fewer interactions per taxon). Interactions between photosynthetic microbes and any other organism became less frequent over the gradient, whereas interactions between plants or soil microfauna and any other organism were more abundant in late succession. Results demonstrate that patterns in diversity and composition do not necessarily relate to patterns in network complexity and suggest that network analyses provide new insight into the ecology of highly diverse, microscopic communities.
Introduction
Microbes are important regulators of ecosystem function, and much research has been directed at testing how soil microbial composition and diversity shift along environmental gradients. Work, thus far, has identified many important abiotic and biotic drivers that structure the composition and diversity of microorganisms in the soil, including pH, salinity, climate, soil nutrients, plant diversity, and plant functional traits (Fierer and Jackson, 2006; Lozupone and Knight, 2007; Nemergut et al., 2010; de Vries et al., 2012; Porazinska et al., 2018). However, how interactions among microbial taxa change across the landscape is much less well understood (Nemergut et al., 2013). Previous studies have used microbial co-occurrence patterns to construct networks and make inferences about microbial interactions (Barberán et al., 2012; Widder et al., 2014; Williams et al., 2014; Morriën et al., 2017; de Vries et al., 2018); however, these methods are limited because microbial co-occurrence can reflect both species interactions (e.g., predation, facilitation, decomposition) and shared environmental niches (co-occurrence due to shared preference for abiotic conditions). Recent advances in hierarchical joint species distribution modeling allow us to parse out the effect of environmental variables to better capture interactions per se (Hui and Poisot, 2016; Ovaskainen et al., 2017); yet they have rarely been applied to microbial datasets (Collins et al., 2018). Here we used this technique to test whether patterns in the complexity of bacterial and eukaryotic networks follow patterns in diversity and compositional change across a fundamental gradient, ecological succession.
While the vast majority of research on community patterns along successional gradients has focused on macroorganisms, shifts in microbial communities are beginning to be addressed. Recent work suggests that the composition of microbes changes considerably and predictably during primary succession (Nemergut et al., 2007; Tarlera et al., 2008; Roy-Bolduc et al., 2015; Poosakkannu et al., 2017), with certain groups like Actinobacteria and N fixers prevalent in early successional stages, and other groups such as Acidobacteria and mycorrhizal taxa dominating later stages (Rime et al., 2015; Poosakkannu et al., 2017; Yarwood and Högberg, 2017). Much work also suggests that the taxonomic diversity (Nemergut et al., 2007; Tarlera et al., 2008; Brown and Jumpponen, 2014; Rime et al., 2015; Poosakkannu et al., 2017) and functional diversity (Tscherko et al., 2003) of soil bacterial and fungal communities generally increase during succession. However, some studies find decreases, no change, or variable patterns for some or all taxa (Sigler and Zeyer, 2002; Brown and Jumpponen, 2014; Dini-Andreote et al., 2014; Poosakkannu et al., 2017).
While patterns in diversity and composition of microbial communities during succession are beginning to emerge, virtually nothing is known about how microbial networks shift across that same gradient. The network structure of communities, independent of diversity and composition, has important implications for resilience (Welti and Joern, 2015; Mandakovic et al., 2018), stability (Neutel et al., 2007), and the efficiency of carbon transfer in an ecosystem (D’Alelio et al., 2016; Morriën et al., 2017). Studies of macroorganisms generally find that interaction networks increase in complexity across succession. In pollination biology (Albrecht et al., 2010; Losapio et al., 2015) and work on animal food webs (Wardle et al., 1995; Neutel et al., 2007), network complexity (typically measured as the number of connections per species) increases due to the increase in quantity and diversity of resources over succession. Experiments have also shown that substrate limitation can constrain the complexity of animal food webs (Chen and Wise, 1999). This same reasoning should hold true for microbial food webs and interaction networks over succession (Bardgett et al., 2005), but it has rarely been studied (but see Dini-Andreote et al., 2014). The diversity of carbon substrates increases over succession (Hooper et al., 2000; Nemergut et al., 2007; Milcu and Manning, 2011). More carbon substrates mean a higher diversity of microorganisms to degrade those substrates and more mutualistic or predatory relationships among taxa can arise (Boer et al., 2005). Furthermore, microbial communities often act in consortia to synergistically degrade complex plant-derived compounds, with some microbes utilizing metabolites or taking advantage of breakdown products of extracellular enzymes produced by other taxa (Lynd et al., 2002; Alessi et al., 2017). Thus, due to the increases in carbon quantity and diversity over succession, we expected the complexity of microbial interactions to increase as well.
Here we combined classical community measures of diversity and species composition with measures of putative species interaction networks using hierarchical Bayesian joint species distribution modeling to ask how microbial community structure changes across a high elevation successional gradient. Our system is a high-elevation alpine landscape in which local variation in environmental conditions (e.g., topography, snow depth) results in strong differences in successional stage across space, as measured by differences in a suite of ecological variables including plant cover, plant diversity, soil nutrients, organic matter content, and microbial biomass. We investigated bacterial and eukaryotic communities, including single-celled eukaryotes, fungi, and soil microfauna, as well as plant communities, to capture multitrophic interactions among core soil organisms. Incorporating multiple taxonomic groups is notable as most soil microbial network research has been limited to a single domain of life, primarily bacteria, and thus has not captured higher-order predation and grazing food webs (but see Steele et al., 2011; Morriën et al., 2017). We hypothesized that over the successional gradient, diversity will increase, species composition will change, and soil microbial networks will become more complex.
Materials and Methods
Study Site
Our study site is a south-facing subnival slope in Green Lakes Valley Watershed at Niwot Ridge LTER (40°3′24″N 105°37′30″W, Supplementary Figure S1). It is approximately 2 km2 in area, ranging between 3610 and 3940 m in elevation, and is composed of talus block slopes, late-melting snowbanks covering unvegetated gravel soils, fellfields, and patches of tundra vegetation. Plant cover ranges from 0% in late-melting snowbanks up to ∼75% (170 individuals m–2) in exposed areas. The site is covered with snow from October to June and the deepest snowfields do not melt fully until August or September. This site was previously sampled in 2007 (King et al., 2010, 2012); we report here on a sampling done in 2015.
Sampling locations were based on a grid in which most plots were ∼50 m apart, and in three targeted sampling areas where plots were 5 m apart (King et al., 2010) (Supplementary Figure S1). We sampled 98 plots that were circular with 1 m radius. Soil samples, for microbial community analysis and biogeochemistry, were collected September 8–17, 2015, 1–2 weeks after the main snowbeds melted out. Three subsamples of soil (∼3 cm diameter, 4 cm depth) were taken from each plot and combined. Samples were homogenized in a ziplock bag, and subsamples for sequencing and microfauna extraction were weighed out the same day and stored at −20 and 4°C, respectively, until processing. Samples were then stored overnight at 4°C, and subsampled for gravimetric soil moisture, water holding capacity (WHC), inorganic N, and total C (TC) and N (TN), microbial C and N, and pH.
Vegetation was surveyed during peak biomass, August 19–September 3, 2015. Plant density was measured by counting all plants within the 3.14 m2 circle and ranged from 0 to 169 plants/m2. Clonal plants and cushion plants were counted as clumps, such that counts represent our best estimate of single genets. Moss was also counted as clumps. Lichen was not included in the plant dataset, because the fungi and algae/cyanobacteria that compose it are measured in the microbe datasets. Plant cover was measured using the point intercept method, assessing the presence of vegetation, rock, and bare soil at 40 points within the plot.
Laboratory Molecular Methods
We chose metabarcoding gene markers specific to each microbial group of interest: 16S for bacteria, ITS for fungi, and 18S for eukaryotes. For eukaryotes, we designated two groups: (1) small, often single-celled, eukaryotes (protists and algae), and (2) soil “microfauna” (small invertebrates including nematodes, tardigrades, rotifers, platyhelminthes, small arthropods, and small annelids). All groups except the microfauna were extracted directly from 0.35 g soil. Because this volume is too small to accurately capture abundance of microfauna, they were first extracted from ∼20 g soil subsamples using Whitehead trays (Porazinska et al., 2014). Briefly, each tray (20 cm × 20 cm × 5 cm) was equipped with a supporting sieve lined with a thin tissue. Soil was spread evenly over the tissue, wetted with 150 ml of water, and left to extract at room temperature for 24 h. Water from each tray was passed through a 38 μm mesh sieve, captured microfauna were transferred to 15 ml sterile Falcon tubes, and allowed to settle overnight at 4°C. The volume was reduced to ∼0.5 ml with sterile disposable pipettes and transferred to bead beating tubes from the PowerSoil DNA isolation kit (MOBIO Laboratories Inc., Carlsbad, CA, United States).
All DNA was extracted using a PowerSoil DNA Isolation Kit according to the manufacture’s protocol. Each sample was PCR amplified twice using primers (515F/806R, ITS1-F/ITS2, and 1391f/EukBr), multiplexing barcodes (Goley), and conditions as adapted by the Earth Microbiome Project1 (Amaral-Zettler et al., 2009; Caporaso et al., 2012). Amplicons were purified and normalized with SequalPrep Normalization Kits (Invitrogen Inc., CA, United States), combined into three single pools of 16S, ITS, and 18S amplicon libraries, and sequenced on three lanes using Illumina technology (MiSeq2000, pair-end 2 × 300 bp) at the BioFrontiers Institute, Boulder, CO, United States.
Sequence Analysis
The raw read sequence data were processed using amplicon sequence variants (ASV) methods in QIIME2 version 2018.2 (Callahan et al., 2017). We demultiplexed and trimmed primers and adapters in QIIME2 and denoised the data and joined paired reads using DADA2. For bacteria and ITS reads that showed poor quality at the ends, we truncated reads where median quality score fell below ∼30; however, DADA2 is robust to low-quality sequence because it incorporates read quality information into its error model. We assigned taxonomy using a pre-trained Naive Bayes classifier. The classifier was trained on the Greengenes 13.8 database (DeSantis et al., 2006) for 16S, UNITE 7.2 (Abarenkov et al., 2010) for ITS, and SILVA-ARB release 111 (Pruesse et al., 2007) for 18S. 16S and 18S sequences were aligned using mafft, and trees were built using FastTree within QIIME2.
Biogeochemistry
We measured a number of environmental variables to characterize the successional stage of each plot and to use as explanatory variables in the joint distribution models. Gravimetric soil moisture was measured on 5 g of soil by drying for 48 h at 60°C. WHC was determined by placing 4 g soil in 15 ml bottom-meshed conical tubes, saturating the soil and allowing to drain, weighing the saturated soil, and then drying at 60°C for 15 h and weighing the dried soil.
Microbial biomass N and C were determined using the chloroform fumigation method (Robertson et al., 1999). Briefly, a 5 g subsample of soil was extracted immediately with 0.5 M K2SO4 and another subsample was fumigated with chloroform to kill microbes and then extracted. Total dissolved nitrogen (TDN) and dissolved organic carbon (DOC) were then analyzed using a Shimadzu total organic carbon analyzer equipped with a TDN module (Shimadzu Scientific Instruments, Inc., Columbia, MD, United States), and microbial biomass N and C were calculated as the difference between the fumigated and unfumigated samples. Inorganic N (NH4+, NO3–) was measured on the unfumigated extracts with a Lachat QuickChem 8500 Flow Injection Analyzer (Lachat Instruments, Loveland, CO, United States) and Synergy 2 Multi-Detection Microplate Reader (BioTek Instruments, Inc., Vinooski, VT, United States).
Total C and TN were measured on 4 g air-dried soil, ground manually with mortar and pestle, and analyzed with a Thermo Finnigan Flash EA 1112 Series CHN analyzer (Thermo Fisher Scientific, Inc.). To measure pH, 2 g of soil was suspended in 3 ml of ultrapure water (Honeywell) and shaken for 1 h, and pH was recorded when the reading stabilized using a calibrated Oakton benchtop pH meter (Oakton Instruments, Vernon Hills, IL, United States).
Snow depth in each plot was calculated from annual snow depth surveys conducted in the Green Lakes Valley from 1997 to 2015 (Losleben, 2002; Morse, 2019). Surveys occur in May at peak snowpack along a grid of random points (mean n = 483) spaced ∼50 m apart. We used kriging interpolation to create a continuous raster surface for each year, conferred snow depth to each plot, and then averaged over the 19 years to reflect long-term snow cover conditions. We also calculated the coefficient of variation (CV) in snow depth for each plot (standard deviation/mean) to quantify variability in snow depth. Elevation (m. a. s. l.) was recorded using a handheld Trimble GPS device with an error rate of 3 m.
Statistical Analysis
Our final data set consisted of the 75 plots that had sufficient sequencing depth shared across all organism groups and available data for all biogeochemical analyses. Prior to analyses, all data sets were rarefied to an even sampling depth (bacteria rarefied to 7987, fungi to 1023 small eukaryote to 871, microfauna to 700). We recognize that this sampling depth is relatively shallow; however, our network analysis focuses on abundant species, thus we believe we are capturing relevant community members (see Supplementary Figure S2 for rarefactions). Prior to network construction, microbial taxa were classified as photosynthetic, heterotrophic, chemoautotrophic, or unknown using a number of taxonomic references (Maddison and Schulz, 2007; Rosenberg et al., 2014a, b, c, d, e, f; Horton et al., 2019). All analyses were done in R (R Core Team, 2017).
We classified plots as early, mid, or late succession by running a PCA of environmental and plant variables previously determined to change over successional gradients in alpine areas (Bekku et al., 2004; Schmidt et al., 2008; Porazinska et al., 2018) (TC, TN, NH43+, NO3–, microbial biomass C, microbial biomass N, pH, WHC, soil moisture, snow depth, elevation, plant density, plant Shannon diversity, and plant cover) using the vegan package (Oksanen et al., 2016) (see Supplementary Figure S3 for the PCA and Supplementary Table S1 for pairwise Pearson correlations). We chose to include all variables in the PCA because succession is associated with changes in all variables; however, results are robust to the choice of variables used to define successional stage. Axis 1 of the PCA explained 50.9% of the variation and was used as a proxy for successional stage (Axis 2 only explained 14.8% of the variation so was not used). Plot scores were extracted and divided into three groups of 25 corresponding to early, mid, and late succession (see Supplementary Figure S1 for photos of plots). We recognize that succession in this system is a continual process and splitting the data into three groups with equal number of plots is not ideal. However, it was necessary to split groups evenly (25 plots in each successional stage) so that sample size would not influence the number of microbial taxa being modeled and the statistical power to detect significant relationships. We chose to split data into three groups, because a sample size of n = 25 plots in each successional stage was reasonable based on the complexity of (df needed for) the joint distribution modeling described below.
To compare diversity of bacteria, fungi, small eukaryotes, and microfauna across succession, Faith’s phylogenetic diversity, ASV richness (Chao1), Pielou’s evenness, rarity, and frequency were calculated. Phylogenetic diversity was calculated using the package picante (Kembel et al., 2010) in R for all groups except fungi, because the fungal ITS sequenced region is highly variable and difficult to align (Moore and Frazer, 2007). Rarity was calculated as the proportion of taxa with relative abundances less than 1/S (Camargo, 1992), where S is the mean ASV richness. Frequency was calculated as the number of plots (out of 25) in which each taxon was present in early, mid, and late successional stages. To test the effect of successional stage on all diversity metrics except frequency, ANOVAs were performed using the gls() function in package nlme (Pinheiro et al., 2019) in R accounting for spatial autocorrelation with a spherical autocorrelation structure. To test the effect of successional stage on frequency, a quasi-poisson regression was performed using the glm() function in R (autocorrelation was not included since the observational unit is taxon rather than plot, quasipoisson models were preferred due to slight overdispersion). Differences among successional stages were tested using Tukey post hoc tests with glht() from the multcomp package (Hothorn et al., 2008).
Compositional shifts across the successional gradient were assessed for each taxonomic group using ordination techniques and visualized using relative abundance barplots and krona plots. For ordination, because microbiome datasets are compositional (relative abundance) and thus have a negative correlation bias, the centered log-ratio (clr) transformation was used (Gloor et al., 2017). We used cmultRepl() in package zCompositions to impute zeros (estimate a small non-zero number for the zeros, because you cannot take the log of zero), and then calculated clr. Redundancy analysis (RDA) was performed on the clr transformed data (also known as the Aitchison distance) (Gloor et al., 2017), and permutation tests were done to test the effect of successional stage on community composition using R package vegan (Oksanen et al., 2019). To visualize compositional shifts over succession at a higher taxonomic level, we calculated relative abundance of different groups of bacteria (phylum level), fungi (phylum level), small eukaryotes (phylum or major clade level, much of the taxonomy is not resolved), and microfauna (phylum level), separating photosynthetic, chemoautotrophic, and heterotrophic taxa within the groups. The effect of successional stage on the most abundant groups was tested using linear models including spatial autocorrelation structure (spherical model), and ANOVAs were performed with P-values adjusted to control the false discovery rate (Benjamini and Hochberg, 1995). Krona interactive plots were also created (Ondov et al., 2011) using the psadd package (Pauvert, 2019) in R and are available in FigShare (Farrer, 2019).
Hierarchical Bayesian joint species distribution modeling was performed using the package boral (Hui and Poisot, 2016) in R. We used clr transformed data to reduce negative correlation bias, and we only included ASVs that occurred in at least 12 plots (out of 25) in each successional stage; thus, we suggest that this is a “core” microbial network focused on interactions among frequent taxa. We fit the joint distribution model with three latent variables and included autocorrelation structure (spherical model) to account for non-independence due to spatial arrangement of plots. We use forward selection to select the four environmental explanatory variables that explained the most variation in our models (snow depth, pH, soil moisture, and CV snow depth), with the realization that many environmental variables were correlated with one another. We used actual environmental variables (as opposed to Axis 1 of the PCA which we used to define successional stage), because we wanted more dimensionality to our analysis and we wanted our results to be generalizable to other systems. For early, mid, and late successional models, there was some correlation among the four environmental explanatory variables (seven out of 18 pairs were correlated: two in early, two in mid, and three in late); however, the absolute values of all correlation coefficients were ≤0.62. When models were run including only environmental variables and autocorrelation structure, those four variables explained an average of 27% of the variance in ASV (clr) abundance, suggesting that it is important to remove their effect. We fit the model using MCMC with 40,000 iterations, a burn-in of 10,000, and a thin of 30. Convergence was checked using Geweke diagnostics and trace plots, and model fit was assessed using Dunn–Smyth residuals and a normal quantile plot (see Supplementary Figure S4). We calculated the residual species-to-species correlation matrix, which represents species correlations after accounting for the effects of the environmental variables. We interpret these correlations as representing species interactions; however, they may also represent unmeasured environmental gradients, which is discussed further in the discussion.
We built networks based on correlations whose 95% posterior credible interval did not overlap zero. Networks were visualized using the package igraph (Csardi and Nepusz, 2006). We calculated linkage density as our metric of complexity (average number of edges per node), and we tabulated number of nodes (ASVs), number of edges (connections between the nodes), number of positive and negative correlations, and the number of bacteria, fungi, small eukaryotes, and soil microfauna included in each network. All code can be found on GitHub: https://github.com/ecfarrer/NWT_MovingUphill5.
Randomizations and Network Assessment
Prior studies have shown that species richness can affect the number of network connections (Faust et al., 2015), and our data show that microbial richness increases across succession. We performed four tests to ensure that the observed patterns of network complexity over the successional gradient were real and not due to statistical artifacts.
First, we assessed our false-positive rate in the networks, because microbe–microbe connections in the network may be due to random chance. We randomized our early, mid, and late successional datasets (n = 10 randomizations for each successional stage) using randomizeMatrix() in picante (Kembel et al., 2010) using the method that randomizes abundances within taxa and maintains species occurrence frequency. Then we ran the joint distribution modeling analysis on this randomized data, assessed the number of significant network connections, and converted this to a percent by dividing by the number of network connections in models using our observed data.
Second, we did simulations to test the effect of taxonomic richness on the number of network connections. Our datasets had two properties that we wanted to preserve in our simulated data, (1) our late successional plots had the highest taxonomic richness, but (2) despite high richness, after applying our frequency cutoff (only modeling species that were present in ≥12 plots, see above), the late succession dataset had fewer taxa that were included in joint distribution modeling compared to early and mid succession datasets (in other words, most of the diversity in late successional plots comprised infrequent taxa). Therefore, we simulated count datasets that varied in species richness by drawing from the probabilities (relative abundances) of the bacterial data from our early and late successional plots. We used a Dirichlet–multinomial distribution using the rmultinom() function and the rdirichlet() function from the MCMCpack package (Martin et al., 2011) as in Faust et al. (2015). For each simulation set, we simulated 25 samples with 2000 reads. We performed joint distribution modeling on these simulated data using only latent variables and assessed the number of significant network connections. Then we repeated the simulations 25 times each for “early” and “late” successional stages. This generated a dataset (n = 50) with variability in mean total taxonomic richness and variability in the richness of frequent species (passing our ≥ 12 frequency cutoff). We assessed the effect of species richness and number of taxa used in modeling on number of network connections and network complexity using linear regression.
Third, because our early, mid, and late successional datasets had different numbers of taxa that passed our frequency cutoff (306, 301, and 273, respectively) and were included in distribution modeling, we subsampled our data to make sure that the network complexity patterns were robust to the number of taxa in the analysis. For the early and mid successional plots, we randomly chose 273 taxa out of the 306 and 301 taxa, respectively, and performed distribution modeling and assessed network complexity (n = 10 times for each successional stage).
Fourth, because the PCA (Supplementary Figure S3) showed that the late successional plots comprised a greater range in variability compared to early or mid successional plots, we did subsetting and modeling to ensure that this did not affect our results. We chose a similar range in PCA Axis 1 (approximately 0.35 Axis 1 units), and randomly chose 12 plots from early, mid, and late successional categories within that range (12 plots was the greatest number of late successional plots that could be obtained within that range). We performed distribution modeling and assessed network complexity for these three reduced models.
Results
Diversity and Composition
Across our microbial datasets, phylogenetic diversity of small eukaryotes and microfauna increased over succession (Figure 1). ASV richness of fungi, small eukaryotes, and microfauna increased over succession (Figure 1 and Supplementary Figure S5), and taxonomic evenness of bacteria, small Eukaryotes, and microfauna increased over succession (Figure 1). The proportion of rare (low abundance) taxa was highest in late succession for small Eukaryotes and soil microfauna, but did not change over succession for bacteria and fungi (Supplementary Figure S5). The average frequency (number of plots) in which taxa were present decreased over succession for bacteria, fungi, and small Eukaryotes; however, the opposite trend was present in the soil microfauna (Supplementary Figure S5).
Figure 1. Diversity and evenness of microbial communities across succession: bacteria (A,E), fungi (B,F), small eukaryotes (C,G), and soil microfauna (D,H). For bacteria, small eukaryotes, and soil microfauna, the diversity metric is Faith’s phylogenetic diversity, which is a measure of taxonomic richness. For Fungi, the metric is ASV richness (Chao1), since the ITS region cannot be aligned to give phylogenetic information. Successional stage does not have a significant effect on bacterial diversity (F2,72 = 1.82, P = 0.17), but does affect fungi (F2,72 = 3.41, P = 0.039), small eukaryotes (F2,72 = 21.01, P < 0.001), and soil microfauna (F2,72 = 54.23, P < 0.001). Successional stage affects evenness of bacteria (F2,72 = 5.28, P = 0.007), small Eukaryotes (F2,72 = 3.36, P = 0.040), and microfauna (F2,72 = 10.65, P < 0.001), but not fungi (F2,72 = 1.96, P = 0.15). Results of Tukey post hoc tests for comparing multiple treatments are shown as letters. Values shown are means and standard errors. Note that the y-axis scales are different for each taxonomic group.
Species composition also changed considerably across the gradient. RDA ordinations showed that successional stage explained 8.2–10.2% (P < 0.001) of the variance in ASV composition of bacteria, fungi, small eukaryote, and soil microfauna communities (Table 1 and Supplementary Figure S6). Relative abundance of some photosynthetic bacteria and eukaryotes (Cyanobacteria, Chlorophyta) decreased over succession (Figure 2). Heterotrophic members of the Verrucomicrobia, Planctomycetes, and Stramenopiles increased in abundance. The Ascomycota increased over succession, whereas the Mortierellomycota decreased by about half (Figure 2B). No significant patterns were seen at the phylum level for the soil microfauna (Figure 2D).
Table 1. Permutational multivariate analysis of variance testing the effect of successional stage on community composition for the four microbial groups.
Figure 2. Effect of successional stage on the relative abundance of major taxonomic groups of bacteria (A), fungi (B), small eukaryotes (C), and soil microfauna (D). Note that only abundant groups are shown in the panels. P and H refer to photosynthetic and heterotrophic groups, respectively. The effect of successional stage on each of the groups was tested using separate anovas, with P-values corrected for false discovery rate using Benjamini–Hochberg: †P < 0.10, ∗P < 0.05, ∗∗P < 0.01, ∗∗∗P < 0.001. Solid lines indicate significant or nearly significant relationships and dashed lines indicate non-significant relationships.
Microbial Networks
Microbial networks decreased in complexity across the successional gradient (Figure 3 and Table 2). Complexity (measured as linkage density, the average number of connections per taxon) decreased from 16.7 to 6.9 to 4.7. Both the number of taxa involved in the networks and the number of correlations decreased over succession (Table 2). The majority of network interactions were positive, and the percentage of positive interactions was highest in late succession (early: 77.9%, mid: 82.0%, late: 88.9%). The number of photosynthetic microbes in networks decreased over succession (Table 2). Soil microfauna increased slightly in networks over succession, with one taxon present in early and mid successional networks and two in the late successional network. Plants were present only once in the networks, in late succession (the grass Trisetum spicatum).
Figure 3. Microbial networks across a successional gradient from early (A), mid (B), and late (C) successional stages. Interactions were assessed using the residual species-to-species correlation matrix from a hierarchical joint distribution model, which removes the effect of environmental covariates on species distributions. Displayed here are the correlations whose 95% credible intervals did not overlap zero. Nodes are colored based on taxonomy. See Table 2 for network statistics.
Table 2. Network statistics for microbial networks across a successional gradient (for the network diagrams, see Figure 3).
Multiple tests showed that our networks were robust and that the observed pattern of increasing taxonomic diversity but decreasing network complexity over succession is not due to statistical artifacts. False positive rates in our networks were very low: networks of randomized data had an average of 1.0 (±0.5), 3.1 (±1.1), and 2.9 (± 0.9) significant microbe–microbe connections for early, mid, and late succession networks, respectively, which suggests 0.029, 0.32, and 0.64% of the interactions in our actual networks may be false positives. Performing joint distribution modeling on simulated datasets with varying taxonomic richness showed that neither the total richness nor the richness of frequent species (those passing our ≥ 12 frequency cutoff used in distribution modeling) affected the number of significant network connections (R2 = 0.0088, P = 0.24, R2 = 0.026, P = 0.14, respectively) or network complexity (R2 = 0.021, P = 0.53, R2 = 0.013, P = 0.27, respectively). Furthermore, when we controlled for the number of taxa passing the frequency cutoff by subsetting the early and mid successional datasets so that the same number of species was included in joint distribution modeling as the late successional dataset, the pattern of reduced complexity across the successional gradient held: the total number of network connections was still higher in early and mid successional areas compared to late (2699 ± 48, 774 ± 20 vs. 452) as was the network complexity (number of interactions per taxon, 14.6 ± 0.2, 6.2 ± 0.1, vs. 4.7). Lastly, when we controlled for the range of PCA axis 1 by running models on reduced datasets in which early, mid, and late successional plots all comprised the same spread along PCA axis 1 (0.35 units each), we found that the pattern of decreasing complexity across the gradient held as well: the total number of connections (678, 156, 60) and the network complexity (10.9, 5.4, 2.9) decreased across the gradient in these reduced but standardized datasets.
Discussion
Testing how different metrics of microbial community structure shift across the landscape is an important step in understanding the complex, yet unseen, majority that regulates ecosystem function. Consistent with our hypothesis, we found that microbial diversity generally increased across the successional gradient and species composition shifted in ways that benefited certain groups (heterotrophic Verrucomicrobia, Ascomycota, heterotrophic Stramenopiles) over others (Cyanobacteria, Mortierellomycota, Chlorophyta). However, contrary to our expectations, the complexity (defined as the average number of interactions per taxon) of core microbial networks (of frequent taxa) decreased across the successional gradient. Results highlight a surprising complexity of photosynthetic and heterotrophic microbial interactions in sparsely vegetated soil and suggest that different analytical metrics (diversity vs. microbial interactions) may yield different conclusions about the complexity of microbial communities.
It is commonly found that microbial diversity increases (Nemergut et al., 2007; Cline and Zak, 2015) and microbial composition changes (Nemergut et al., 2007; Dini-Andreote et al., 2014) over succession. We found that the diversity of fungi, small eukaryotes, and soil microfauna (but not bacteria) increases over succession, and we observed marked changes in composition in all taxonomic groups. The lack of (significant) increase in bacterial diversity is somewhat surprising, but previous work has shown that bacteria dominate the active community of unvegetated soils during the summer months at this site (Ley and Schmidt, 2002) likely because they have a much broader range of physiologies (e.g., chemoautotrophy, N-fixation) than do fungi and other eukaryotes (Schmidt et al., 2014). Increases in eukaryotic microbial diversity are likely due, in part, to increases in resource availability (Waldrop et al., 2006), such that a greater number of taxa are able to meet minimum resource requirements as resources increase. Overall, the harsh environmental conditions (low nutrients, low moisture) in early successional areas may limit microbes in the same way they limit plant growth. Also, over successional gradients, the number of complex carbon molecules in the system increases – saccharides, cellulose, phenolics, lignins, tannins – and one of the main controls on microbial diversity is the diversity of carbon substrates (Zhou et al., 2002). The shifting carbon inputs and changing abiotic conditions are also likely driving the compositional shifts observed. Some of the changes in relative abundance observed here agree with findings from other studies, such as an increase in Verrucomicrobia (Nemergut et al., 2007) and Acidobacteria (Nemergut et al., 2007; Yarwood and Högberg, 2017) over succession and an abundance of cyanobacteria (Schmidt et al., 2008; Zumsteg et al., 2012; Yarwood and Högberg, 2017) and algae (Kaštovská et al., 2005) in early successional stages. However, some of our results differ from prior studies that show decreases in Actinobacteria (Zumsteg et al., 2012; Poosakkannu et al., 2017) and shifts from Ascomycota to Basiodiomycota (Zumsteg et al., 2012) with succession in alpine and arctic areas, suggesting that succession may proceed differently at different sites even at a broad taxonomic level.
Despite general increases in diversity, we found that core microbial interaction networks decreased in complexity across the successional gradient. This same pattern was observed along a salt marsh chronosequence, where early successional stages displayed more network complexity (Dini-Andreote et al., 2014). However, other studies show that microbial network complexity increases over primary succession in forests (Yarwood and Högberg, 2017) and secondary succession in old fields (Morriën et al., 2017). Furthermore, studies show that carbon inputs increase network complexity during the short-term development of rhizosphere communities (Shi et al., 2016) and in a CO2 enrichment experiment (Zhou et al., 2011; Tu et al., 2015). All of these studies (with the exception of Morriën et al., 2017) focus only on one taxonomic group (bacteria or fungi), so they do not capture higher order trophic interactions and cascades.
In our system, one explanation for the decrease in network complexity over succession is the shift in primary producers and carbon resources across the gradient, from a dominance of algal/cyanobacterial and recalcitrant aeolian carbon inputs early in succession to a dominance of labile carbon compounds from root exudates as plant density increases. Our results show that early successional soils contain abundant photosynthetic bacteria and eukaryotes including Cyanobacteria, Chorophyta, and Stramenopiles (e.g., Diatomea and Chrysophyceae), and previous work at our study site has measured significant microbial carbon fixation in these barren, plant-free soils (Freeman et al., 2009b). These primary producers could form the base of food chains and interact with heterotrophs that feed on algal exudates, senescent algae, and algal detritus (Cole, 1982). In fact, studies in other systems show that the chemical constituents of algal exudates differ by species (Aluwihare and Repeta, 1999) and this affects the composition of bacterial heterotrophs (Kalscheur et al., 2012). Interactions between heterotrophic organisms and algae and cyanobacteria would result in high network complexity in early succession that would decline in later succession as increasing plant cover shades out photosynthetic microbes. This is borne out in the networks: the number of interactions between (frequent) bacterial heterotrophs and (frequent) photoautotrophic microbes declines over succession from 205 (early) to 55 (mid) to 1 (late). Additionally, the fourth most highly connected microbial taxon (potentially a type of “hub”) in the early successional networks was a photoautotrophic microbe in the Chloroflexi, whereas all highly connected taxa in mid and late successional networks were non-photosynthetic.
Further work at our study site shows that aeolian-deposited plant litter, particularly windblown pollen grains, is another important carbon source in early succession (Freeman et al., 2009a; Naff et al., 2013). Data from our study site indicate that barren talus soils contain a higher proportion of taxa that mineralize complex organic matter, like pollen, compared to vegetated soils (Ley et al., 2001). The breakdown of recalcitrant organic matter can involve numerous microbial taxa in decomposition, as bacteria that degrade complex substrates provide resources for bacteria that feed on metabolic byproducts or take advantage of monomers released by extracellular enzymes (Rakoff-Nahoum et al., 2014; Datta et al., 2016). This could result in complex interaction networks in early succession that would disappear as more labile plant root exudates dominate the system; however, this has not been experimentally tested. In our early succession networks, we found numerous positive interactions (146) between heterotrophic bacteria and members of the Ktedonobacteria (Phylum Chloroflexi), a group that is known to be important in degrading complex polymers, like pollen (Hug et al., 2013; King and King, 2014), and the number of heterotrophic bacteria–Ktedonobacteria interactions decrease over succession (33 in mid succession, four in late succession). Interestingly, the relative abundance of Ktedonobacteria does not change over succession, just their involvement in networks.
A major limitation in this study, and many studies of soil microbial communities, is the limited functional information about microbial taxa. This makes interpretation of network patterns across succession difficult. We took a first step by classifying taxa by mode of nutrition (photosynthetic/chemoautotrophic/heterotrophic) and by looking at patterns in groups for which we have information in our system (e.g., Ktedonobacteria). However, many interactions in our networks remain unexplained because we do not know the function of one or both of the interacting taxa (for 37% of the interactions, we do not even know if interacting organisms are photosynthetic or heterotrophic). As more functional information becomes available, it may open up new hypotheses explaining the observed drop in interaction complexity over succession.
Network Interactions
Although we know very little about many of the taxa in the networks, we can highlight some interactions present in the network as proof of concept that the networks display realistic microbial relationships. In late successional networks, the grass Tristum spicatum was positively associated with a bacterium in the Sinobacteraceae, which has been shown to be a root associate of other plants (Tanaka et al., 2018). Trisetum was also positively associated with the Archaean, Candidatus Nitrososphaera, an ammonia oxidizing prokaryote (the first step in nitrification) (Rosenberg et al., 2014f). Archaeans are increasingly recognized as important plant symbionts especially in alpine areas (Taffner et al., 2018), and evidence suggests that Trisetum tends to prefer nitrate uptake to ammonium (Miller and Bowman, 2002). In late successional communities, we found relationships (positive and negative) between a bacterial-feeding nematode in the Plectidae (Yeates et al., 1993) and 28 different bacteria; both positive and negative relationships could be suggestive of consumption, depending on whether the predator limits prey availability or is limited by prey availability (Halliwell and Macdonald, 1996). Similarly, in early successional networks, we found relationships between a predatory nematode in the Nygolaimidae (Yeates et al., 1993) and 31 different eukaryotes, bacteria, and fungi. We also found positive correlations in mid succession communities between a bacterial feeding Heterolobosean amoeba (Brown et al., 2012; Pánek and Cepicka, 2012) and 15 bacterial taxa. And in early successional communities, we found correlations among two filter-feeding (subclass Scuticociliatia and Hypotrichia) ciliates and 41 different bacteria, which likely indicate feeding relationships (Lee and Capriulo, 1990; Jurgens and Simek, 2000). Overall, despite the lack of functional information for many of the microbial taxa in our networks, the interactions discussed here suggest that the methods used produce interpretable results.
Limitations
While the hierarchical joint distribution modeling yielded realistic microbial relationships, had low false positive rates, and was not highly influenced by species richness like other correlational approaches (Faust et al., 2015), there are limitations. As with any correlation-based analysis, a key limitation in this work is that we cannot definitively conclude that all correlations represent interactions. While we attempted to isolate interactions per se by using hierarchical joint distribution modeling and removing the effect of environmental covariates, the residual correlations could represent shared niches that were not parceled out by the four chosen environmental covariates. For example, in early succession, there were five positive correlations among photosynthetic microbes (bacteria and algae), which might indicate that they share the same high light niche. Similarly, we did not include specific carbon substrates in the joint distribution model, so correlations among heterotrophic microbes may indicate shared niches specializing on a particular substrate (e.g., chitin, cellulose). We suggest that the putative interactions discovered by combining joint distribution modeling with network analysis be taken as hypotheses, which can be tested with future studies.
Our network analysis, by statistical necessity, only focuses on interactions among core, frequent microbial taxa. The analysis does capture a fairly wide range of taxon abundance: for example, for bacteria in early succession, the average relative abundance of taxa in the network ranged two orders of magnitude from 0.00051 to 0.057, and species with higher relative abundance were no more likely to be included in the network (i.e., have interactions with other microbes) vs. not. That said, there were many infrequent taxa that we could not attempt to model: in early succession, there were a total of 8442 taxa (bacteria, fungi, small eukaryotes, soil microfauna, plants) in the plots, but only 306 were frequent enough (present in > 11 plots) to be considered in the joint distribution modeling. It is difficult to predict how interactions with infrequent taxa would change network structure. Work on plants suggests that rare species may be facilitated more than dominant species (Choler et al., 2001); however, predation may have the opposite pattern (rare species may escape predation) (Murdoch and Oaten, 1975). Understanding the functioning of the rare biosphere is one of the major challenges in microbial ecology today, and new technologies and analytical methods are being developed to probe these rare community members (Jia et al., 2018).
Integrating diversity, species composition, and interaction networks gives us a comprehensive assessment of how systems change across succession at the community level; however, our methods do not address how ecosystem properties shift across succession. While there is a large literature on the link between biodiversity and ecosystem function (e.g., Hooper et al., 2005), there is scant information on the impact of network complexity on ecosystem level properties. One previous study showed that microbial network complexity is associated with increased efficiency of carbon uptake (Morriën et al., 2017), and by attempting to isolate interactions per se, our work provides some understanding of energy flow through an ecosystem since many interactions represent trophic relationships. The implications of network complexity for nutrient cycling and ecosystem function are a topic that is ripe for future studies.
Conclusion
Microbes control many key aspects of ecosystem function, and ecologists have well documented the changes in diversity and composition of microbial communities across environmental gradients. However, less is known about how the network structure of microbes shifts across the landscape, and whether patterns in the complexity of species interactions follow patterns in complexity of diversity. Here we find that diversity increases across a successional gradient but network complexity declines. While much work remains to be done regarding the interpretation and implications of ecological network structure, using hierarchical joint distribution modeling to generate interaction networks is an important first step toward elucidating potential species interactions and important microbial players. The more we learn about the organization of microbial communities, the more we can leverage them for the many ecosystem services they provide.
Data Availability Statement
The raw 16S, ITS, and 18S sequence data were deposited in the NCBI Sequence Read Archive (SRA): PRJNA587150.
Author Contributions
EF, DP, MS, AK, KS, and SKS conceived the study. EF, DP, CB, SAS, JS, and CW collected field data and samples. EF, DP, CB, SAS, and CW performed laboratory biogeochemistry and molecular analyses. EF performed bioinformatics and data analysis and wrote the manuscript. All authors edited the manuscript and approved the final version of the manuscript.
Funding
Funding was provided by the National Science Foundation (DEB-1457827) and Niwot Ridge LTER Program (NSF DEB-1637686).
Conflict of Interest
AK was employed by Ecological Consulting (Knoxville, TN, United States).
The remaining 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 Josh Grinath, Lauren Hallett, and Jared Huxley for assistance in the field and lab. We thank Noah Fierer and Jessica Henley for help with molecular work. Five reviewers gave thoughtful comments on the manuscript. Logistical support was provided by the Niwot Ridge LTER Program and CU Mountain Research Station.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmicb.2019.02887/full#supplementary-material
Footnotes
References
Abarenkov, K., Nilsson, R. H., Larsson, K. H., Alexander, I. J., Eberhardt, U., Erland, S., et al. (2010). The UNITE database for molecular identification of fungi – recent updates and future perspectives. N. Phytol. 186, 281–285. doi: 10.1111/j.1469-8137.2009.03160.x
Albrecht, M., Riesen, M., and Schmid, B. (2010). Plant–pollinator network assembly along the chronosequence of a glacier foreland. Oikos 119, 1610–1624. doi: 10.1111/j.1600-0706.2010.18376.x
Alessi, A. M., Bird, S. M., Bennett, J. P., Oates, N. C., Li, Y., Dowle, A. A., et al. (2017). Revealing the insoluble metasecretome of lignocellulose-degrading microbial communities. Sci. Rep. 7:2356. doi: 10.1038/s41598-017-02506-5
Aluwihare, L. I., and Repeta, D. J. (1999). A comparison of the chemical characteristics of oceanic DOM and extracellular DOM produced by marine algae. Mar. Ecol. Prog. Ser. 186, 105–117. doi: 10.3354/meps186105
Amaral-Zettler, L. A., McCliment, E. A., Ducklow, H. W., and Huse, S. M. (2009). A method for studying protistan diversity using massively parallel sequencing of V9 hypervariable regions of small-subunit ribosomal RNA genes. PLoS One 4:e6372. doi: 10.1371/journal.pone.0006372
Barberán, A., Bates, S. T., Casamayor, E. O., and Fierer, N. (2012). Using network analysis to explore co-occurrence patterns in soil microbial communities. ISME J. 6, 343–351. doi: 10.1038/ismej.2011.119
Bardgett, R. D., Bowman, W. D., Kaufmann, R., and Schmidt, S. K. (2005). A temporal approach to linking aboveground and belowground ecology. Trends Ecol. Evol. 20, 634–641. doi: 10.1016/j.tree.2005.08.005
Bekku, Y. S., Nakatsubo, T., Kume, A., and Koizumi, H. (2004). Soil microbial biomass, respiration rate, and temperature dependence on a successional glacier foreland in Ny-Ålesund, Svalbard. Bioone 36, 395–399. doi: 10.1657/1523-0430(2004)036
Benjamini, Y., and Hochberg, Y. (1995). Controlling the false discovery rate: a practical and powerful approach to multiple testing. J. R. Stat. Soc. Ser. B 57, 289–300. doi: 10.1111/j.2517-6161.1995.tb02031.x
Boer, W. D., Folman, L. B., Summerbell, R. C., and Boddy, L. (2005). Living in a fungal world: impact of fungi on soil bacterial niche development. FEMS Microbiol. Rev. 29, 795–811. doi: 10.1016/j.femsre.2004.11.005
Brown, M. W., Silberman, J. D., and Spiegel, F. W. (2012). A contemporary evaluation of the acrasids (Acrasidae, Heterolobosea, Excavata). Eur. J. Protistol. 48, 103–123. doi: 10.1016/j.ejop.2011.10.001
Brown, S. P., and Jumpponen, A. (2014). Contrasting primary successional trajectories of fungi and bacteria in retreating glacier soils. Mol. Ecol. 23, 481–497. doi: 10.1111/mec.12487
Callahan, B. J., McMurdie, P. J., and Holmes, S. P. (2017). Exact sequence variants should replace operational taxonomic units in marker-gene data analysis. ISME J. 11:2639. doi: 10.1038/ismej.2017.119
Camargo, J. A. (1992). Can dominance influence stability in competitive interactions? Oikos 64, 605–609. doi: 10.1002/ecy.1797
Caporaso, J. G., Lauber, C. L., Walters, W. A., Berg-Lyons, D., Huntley, J., Fierer, N., et al. (2012). Ultra-high-throughput microbial community analysis on the Illumina HiSeq and MiSeq platforms. ISME J. 6, 1621–1624. doi: 10.1038/ismej.2012.8
Chen, B., and Wise, D. H. (1999). Bottom-up limitation of predaceous arthropods in a detritus-based terrestrial food web. Ecology 80, 761–772. doi: 10.1890/0012-9658(1999)080
Choler, P., Michalet, R., and Callaway, R. M. (2001). Facilitation and competition on gradients in alpine plant communities. Ecology 82, 3295–3308. doi: 10.1890/0012-9658(2001)082
Cline, L. C., and Zak, D. R. (2015). Soil microbial communities are shaped by plant-driven changes in resource availability during secondary succession. Ecology 96, 3374–3385. doi: 10.1890/15-0184.1
Cole, J. J. (1982). Interactions between bacteria and algae in aquatic ecosystems. Annu. Rev. Ecol. Syst. 13, 291–314. doi: 10.1146/annurev.es.13.110182.001451
Collins, C. G., Stajich, J. E., Weber, S. E., Pombubpa, N., and Diez, J. M. (2018). Shrub range expansion alters diversity and distribution of soil fungal communities across an alpine elevation gradient. Mol. Ecol. 27, 2461–2476. doi: 10.1111/mec.14694
Csardi, G., and Nepusz, T. (2006). The igraph software package for complex network research. InterJ., Comp. Syst. 1695, 1–9.
D’Alelio, D., Libralato, S., Wyatt, T., and Ribera d’Alcalà, M. (2016). Ecological-network models link diversity, structure and function in the plankton food-web. Sci. Rep. 6:21806. doi: 10.1038/srep21806
Datta, M. S., Sliwerska, E., Gore, J., Polz, M. F., and Cordero, O. X. (2016). Microbial interactions lead to rapid micro-scale successions on model marine particles. Nat. Commun. 7:11965. doi: 10.1038/ncomms11965
de Vries, F. T., Griffiths, R. I., Bailey, M., Craig, H., Girlanda, M., Gweon, H. S., et al. (2018). Soil bacterial networks are less stable under drought than fungal networks. Nat. Commun. 9, 3033–3033. doi: 10.1038/s41467-018-05516-7
de Vries, F. T., Manning, P., Tallowin, J. R. B., Mortimer, S. R., Pilgrim, E. S., Harrison, K. A., et al. (2012). Abiotic drivers and plant traits explain landscape-scale patterns in soil microbial communities. Ecol. Lett. 15, 1230–1239. doi: 10.1111/j.1461-0248.2012.01844.x
DeSantis, T. Z., Hugenholtz, P., Larsen, N., Rojas, M., Brodie, E. L., Keller, K., et al. (2006). Greengenes, a chimera-checked 16S rRNA gene database and workbench compatible with ARB. Appl. Environ. Microbiol. 72, 5069–5072. doi: 10.1128/aem.03006-05
Dini-Andreote, F., de Cassia Pereira e Silva, M., Triado-Margarit, X., Casamayor, E. O., van Elsas, J. D., and Salles, J. F. (2014). Dynamics of bacterial community succession in a salt marsh chronosequence: evidences for temporal niche partitioning. ISME J. 8, 1989–2001. doi: 10.1038/ismej.2014.54
Farrer, E. C. (2019). Krona plots from Soil microbial networks shift across a high-elevation successional gradient. Front. Microbiol. doi: 10.6084/m9.figshare.9971552.v1
Faust, K., Lima-Mendez, G., Lerat, J. S., Sathirapongsasuti, J. F., Knight, R., Huttenhower, C., et al. (2015). Cross-biome comparison of microbial association networks. Front. Microbiol. 6:1200. doi: 10.3389/fmicb.2015.01200
Fierer, N., and Jackson, R. B. (2006). The diversity and biogeography of soil bacterial communities. Proc. Natl. Acad. Sci. U.S.A. 103, 626–631. doi: 10.1073/pnas.0507535103
Freeman, K. R., Martin, A. P., Karki, D., Lynch, R. C., Mitter, M. S., Meyer, A. F., et al. (2009a). Evidence that chytrids dominate fungal communities in high-elevation soils. Proc. Natl. Acad. Sci. U.S.A. 106, 18315–18320. doi: 10.1073/pnas.0907303106
Freeman, K. R., Pescador, M. Y., Reed, S. C., Costello, E. K., Robeson, M. S., and Schmidt, S. K. (2009b). Soil CO2 flux and photoautotrophic community composition in high-elevation, ‘barren’ soil. Environ. Microbiol. 11, 674–686. doi: 10.1111/j.1462-2920.2008.01844.x
Gloor, G. B., Macklaim, J. M., Pawlowsky-Glahn, V., and Egozcue, J. J. (2017). Microbiome datasets are compositional: And this is not optional. Front. Microbiol. 8:2224. doi: 10.3389/fmicb.2017.02224
Halliwell, E. C., and Macdonald, D. W. (1996). American mink Mustela vison in the upper thames catchment: relationship with selected prey species and den availability. Biol. Conserv. 76, 51–56. doi: 10.1016/0006-3207(95)00072-0
Hooper, D. U., Bignell, D. E., Brown, V. K., Brussard, L., Dangerfield, J. M., Wall, D. H., et al. (2000). Interactions between aboveground and belowground biodiversity in terrestrial ecosystems: patterns, mechanisms, and FeedbacksWe assess the evidence for correlation between aboveground and belowground diversity and conclude that a variety of mechanisms could lead to positive, negative, or no relationship—depending on the strength and type of interactions among species. Bioscience 50, 1049–1061.
Hooper, D. U., Chapin, F. S., Ewel, J. J., Hector, A., Inchausti, P., Lavorel, S., et al. (2005). Effects of biodiversity on ecosystem functioning: a consensus of current knowledge. Ecol. Monogr. 75, 3–35. doi: 10.1111/brv.12110
Horton, T., Kroh, A., Ahyong, S., Bailly, N., Boyko, C. B., Brandão, S. N., et al. (2019). World Register of Marine Species (WoRMS). Ostend: WoRMS Editorial Board.
Hothorn, T., Bretz, F., and Westfall, P. (2008). Simultaneous inference in general parametric models. Biometr. J. 50, 346–363. doi: 10.1002/bimj.200810425
Hug, L. A., Castelle, C. J., Wrighton, K. C., Thomas, B. C., Sharon, I., Frischkorn, K. R., et al. (2013). Community genomic analyses constrain the distribution of metabolic traits across the Chloroflexi phylum and indicate roles in sediment carbon cycling. Microbiome 1:22. doi: 10.1186/2049-2618-1-22
Hui, F. K. C., and Poisot, T. (2016). BORAL – bayesian ordination and regression analysis of multivariate abundance data in R. Methods Ecol. Evol. 7, 744–750. doi: 10.1111/2041-210x.12514
Jia, X., Dini-Andreote, F., and Falcão Salles, J. (2018). Community assembly processes of the microbial rare biosphere. Trends Microbiol. 26, 738–747. doi: 10.1016/j.tim.2018.02.011
Jurgens, K., and Simek, K. (2000). Functional response and particle size selection of Halteria cf. grandinella, a common freshwater oligotrichous ciliate. Aqua. Microb. Ecol. 22, 57–68. doi: 10.3354/ame022057
Kalscheur, K. N., Rojas, M., Peterson, C. G., Kelly, J. J., and Gray, K. A. (2012). Algal exudates and stream organic matter influence the structure and function of denitrifying bacterial communities. Microb. Ecol. 64, 881–892. doi: 10.1007/s00248-012-0091-1
Kaštovská, K., Elster, J., Stibal, M., and Šantrùèková, H. (2005). Microbial assemblages in soil microbial succession after glacial retreat in Svalbard (High Arctic). Microb. Ecol. 50:396. doi: 10.1007/s00248-005-0246-4
Kembel, S. W., Cowan, P. D., Helmus, M. R., Cornwell, W. K., Morlon, H., Ackerly, D. D., et al. (2010). Picante: R tools for integrating phylogenies and ecology. Bioinformatics 26, 1463–1464. doi: 10.1093/bioinformatics/btq166
King, A. J., Farrer, E. C., Suding, K. N., and Schmidt, S. K. (2012). Cooccurrence patterns of plants and soil bacteria in the high-alpine subnival zone track environmental harshness. Front. Microbiol. 3:347. doi: 10.3389/fmicb.2012.00347
King, A. J., Freeman, K. R., McCormick, K. F., Lynch, R. C., Lozupone, C., Knight, R., et al. (2010). Biogeography and habitat modelling of high-alpine bacteria. Nat. Commun. 1:53. doi: 10.1038/ncomms1055
King, C. E., and King, G. M. (2014). Description of Thermogemmatispora carboxidivorans sp. nov., a carbon-monoxide-oxidizing member of the class Ktedonobacteria isolated from a geothermally heated biofilm, and analysis of carbon monoxide oxidation by members of the class Ktedonobacteria. Int. J. Syst. Evol. Microbiol. 64, 1244–1251. doi: 10.1099/ijs.0.059675-0
Lee, J. J., and Capriulo, G. M. (1990). “The Ecology of Marine Protozoa: an Overview,” in Ecology of Marine Protozoa, ed. G. M. Capriulo (New York, NY: Oxford University Press), 1–15
Ley, R. E., Lipson, D. A., and Schmidt, S. K. (2001). Microbial biomass levels in barren and vegetated high altitude talus soils. Soil Sci. Soc. Am. J. 65, 111–117.
Ley, R. E., and Schmidt, S. K. (2002). Fungal and bacterial responses to phenolic compounds and amino acids in high altitude barren soils. Soil Biol. Biochem. 34, 989–995. doi: 10.1016/s0038-0717(02)00032-9
Losapio, G., Jordán, F., Caccianiga, M., and Gobbi, M. (2015). Structure-dynamic relationship of plant–insect networks along a primary succession gradient on a glacier foreland. Ecol. Modell. 314, 73–79. doi: 10.1016/j.ecolmodel.2015.07.014
Losleben, M. (2002). Snow Survey for Green Lakes Valley, Niwot Ridge LTER, Colorado. Colorado: Environmental Data Initiative.
Lozupone, C. A., and Knight, R. (2007). Global patterns in bacterial diversity. Proc. Natl. Acad. Sci. U.S.A. 104, 11436–11440. doi: 10.1073/pnas.0611525104
Lynd, L. R., Weimer, P. J., van Zyl, W. H., and Pretorius, I. S. (2002). Microbial Cellulose Utilization: fundamentals and Biotechnology. Microbiol. Mol. Biol. Rev. 66, 506–577. doi: 10.1128/mmbr.66.3.506-577.2002
Maddison, D. R., and Schulz, K. S. (2007). The Tree of Life Web Project. Available at: http://tolweb.org (accessed November 7, 2019).
Mandakovic, D., Rojas, C., Maldonado, J., Latorre, M., Travisany, D., Delage, E., et al. (2018). Structure and co-occurrence patterns in microbial communities under acute environmental stress reveal ecological factors fostering resilience. Sci. Rep. 8:5875. doi: 10.1038/s41598-018-23931-0
Martin, A. D., Quinn, K. M., and Park, J. H. (2011). MCMCpack: Markov Chain Monte Carlo in R. J. Stat. Softw. 42. doi: 10.18637/jss.v042.i09
Milcu, A., and Manning, P. (2011). All size classes of soil fauna and litter quality control the acceleration of litter decay in its home environment. Oikos 120, 1366–1370. doi: 10.1111/j.1600-0706.2010.19418.x
Miller, A. E., and Bowman, W. D. (2002). Variation in Nitrogen-15 natural abundance and nitrogen uptake traits among co-occurring alpine species: do species partition by nitrogen form? Oecologia 130, 609–616. doi: 10.1007/s00442-001-0838-8
Morriën, E., Hannula, S. E., Snoek, L. B., Helmsing, N. R., Zweers, H., de Hollander, M., et al. (2017). Soil networks become more connected and take up more carbon as nature restoration progresses. Nat. Commun. 8:14349. doi: 10.1038/ncomms14349
Morse, J. (2019). Annual Snow Survey, Green Lakes Valley, Niwot Ridge, Colorado, 2013-Ongoing. Environmental Data Initiative. Available at: https://doi.org/10.6073/pasta/50e635a1ac210481b1ef9a5df3a3df60 (accessed December 3, 2019).
Murdoch, W. W., and Oaten, A. (1975). “Predation and Population Stability,” in Advances in Ecological Research, ed. A. MacFadyen (Cambridge, MA: Academic Press), 1–131. doi: 10.1016/s0065-2504(08)60288-3
Naff, C. S., Darcy, J. L., and Schmidt, S. K. (2013). Phylogeny and biogeography of an uncultured clade of snow chytrids. Environ. Microbiol. 15, 2672–2680. doi: 10.1111/1462-2920.12116
Nemergut, D. R., Anderson, S. P., Cleveland, C. C., Martin, A. P., Miller, A. E., Seimon, A., et al. (2007). Microbial community succession in an unvegetated, recently deglaciated soil. Microb. Ecol. 53, 110–122. doi: 10.1007/s00248-006-9144-7
Nemergut, D. R., Cleveland, C. C., Wieder, W. R., Washenberger, C. L., and Townsend, A. R. (2010). Plot-scale manipulations of organic matter inputs to soils correlate with shifts in microbial community composition in a lowland tropical rain forest. Soil Biol. Biochem. 42, 2153–2160. doi: 10.1016/j.soilbio.2010.08.011
Nemergut, D. R., Schmidt, S. K., Fukami, T., Bilinski, T. M., Stanish, L. F., Knelman, J. E., et al. (2013). Patterns and processes of microbial community assembly. Microbiol. Mol. Biol. Rev. 77, 342–356. doi: 10.1128/MMBR.00051-12
Neutel, A. M., Heesterbeek, J. A. P., van de Koppel, J., Hoenderboom, G., Vos, A., Kaldeway, C., et al. (2007). Reconciling complexity with stability in naturally assembling food webs. Nature 449, 599–602. doi: 10.1038/nature06154
Oksanen, J., Blanchet, F. G., Friendly, M., Kindt, R., Legendre, P., McGlinn, D., et al. (2016). vegan: Community Ecology Package. R Package Version 2.4–1.
Oksanen, J., Blanchet, F. G., Friendly, M., Kindt, R., Legendre, P., McGlinn, D., et al. (2019). vegan: Community Ecology Package. R Package Version 2.5–5.
Ondov, B. D., Bergman, N. H., and Phillippy, A. M. (2011). Interactive metagenomic visualization in a Web browser. BMC Bioinformatics 12:385. doi: 10.1186/1471-2105-12-385
Ovaskainen, O., Tikhonov, G., Norberg, A., Guillaume Blanchet, F., Duan, L., Dunson, D., et al. (2017). How to make more out of community data? a conceptual framework and its implementation as models and software. Ecol. Lett. 20, 561–576. doi: 10.1111/ele.12757
Pánek, T., and Cepicka, I. (2012). “Diversity of heterolobosea,” in Genetic Diversity in Microorganisms, ed. M. Caliskan (Rijeka: Intech Open), 3–26.
Pauvert, C. (2019). psadd: Additions to Phyloseq Package for Microbiome Analysis, R Package Version 0.1.2.
Pinheiro, J., Bates, D., DebRoy, S., Sarkar, D., and Team, R. C. (2019). Nlme: Linear and Nonlinear Mixed Effects Models, R Package Version 3.1–141.
Poosakkannu, A., Nissinen, R., Männistö, M., and Kytöviita, M. M. (2017). Microbial community composition but not diversity changes along succession in arctic sand dunes. Environ. Microbiol. 19, 698–709. doi: 10.1111/1462-2920.13599
Porazinska, D. L., Farrer, E. C., Spasojevic, M. J., Sartwell, S. A., Smith, J. G., White, C. T., et al. (2018). Plant diversity and density predict belowground diversity and function in an early successional alpine ecosystem. Ecology 99, 1942–1952. doi: 10.1002/ecy.2420
Porazinska, D. L., Morgan, M. J., Gaspar, J. M., Court, L. N., Hardy, C. M., and Hodda, M. (2014). Discrimination of plant-parasitic nematodes from complex soil communities using ecometagenetics. Phytopathology 104, 749–761. doi: 10.1094/PHYTO-08-13-0236-R
Pruesse, E., Quast, C., Knittel, K., Fuchs, B. M., Ludwig, W., Peplies, J., et al. (2007). SILVA: a comprehensive online resource for quality checked and aligned ribosomal RNA sequence data compatible with ARB. Nucleic Acids Res. 35, 7188–7196. doi: 10.1093/nar/gkm864
R Core Team (2017). R: A Language and Environment for Statistical Computing. Vienna: R Foundation for Statistical Computing.
Rakoff-Nahoum, S., Coyne, M. J., and Comstock, L. E. (2014). An ecological network of polysaccharide utilization among human intestinal symbionts. Curr. Biol. 24, 40–49. doi: 10.1016/j.cub.2013.10.077
Rime, T., Hartmann, M., Brunner, I., Widmer, F., Zeyer, J., and Frey, B. (2015). Vertical distribution of the soil microbiota along a successional gradient in a glacier forefield. Mol. Ecol. 24, 1091–1108. doi: 10.1111/mec.13051
Robertson, G. P., Coleman, D. C., Bledsoe, C. S., and Sollins, P. (1999). Standard Soil Methods for Long-Term Ecological Research. New York, NY: Oxford University Press.
Rosenberg, E., DeLong, E. F., Lory, S., Stackebrandt, E., et al. (2014a). The Prokaryotes: Actinobacteria, 4th Edn. Berlin: Springer-Verlag.
Rosenberg, E., DeLong, E. F., Lory, S., Stackebrandt, E., et al. (2014b). The Prokaryotes: Alphaproteobacteria and Betaproteobacteria, 4th Edn, Berlin: Springer-Verlag.
Rosenberg, E., DeLong, E. F., Lory, S., Stackebrandt, E., et al. (2014c). The Prokaryotes: Deltaproteobacteria and Epsilonproteobacteria, 4th Edn, Berlin: Springer-Verlag.
Rosenberg, E., DeLong, E. F., Lory, S., Stackebrandt, E., et al. (2014d). The Prokaryotes: Firmicutes and Tenericutes, 4th Edn, Berlin: Springer-Verlag.
Rosenberg, E., DeLong, E. F., Lory, S., Stackebrandt, E., et al. (2014e). The Prokaryotes: Gammaproteobacteria, 4th Edn, Berlin: Springer-Verlag.
Rosenberg, E., DeLong, E. F., Lory, S., Stackebrandt, E., et al. (2014f). The Prokaryotes: Other Major Lineages of Bacteria and The Archaea, 4th Edn, Berlin: Springer-Verlag.
Roy-Bolduc, A., Bell, T. H., Boudreau, S., and Hijri, M. (2015). Comprehensive sampling of an isolated dune system demonstrates clear patterns in soil fungal communities across a successional gradient. Environ. Microbiol. Rep. 7, 839–848. doi: 10.1111/1758-2229.12312
Schmidt, S. K., Nemergut, D. R., Darcy, J. L., and Lynch, R. (2014). Do bacterial and fungal communities assemble differently during primary succession? Mol. Ecol. 23, 254–258. doi: 10.1111/mec.12589
Schmidt, S. K., Reed, S. C., Nemergut, D. R., Cleveland, C. C., Weintraub, M. N., Hill, A. W., et al. (2008). The earliest stages of ecosystem succession in high-elevation (5000 metres above sea level), recently deglaciated soils. Proc. R. Soc. B 275, 2793–2802. doi: 10.1098/rspb.2008.0808
Shi, S., Nuccio, E. E., Shi, Z. J., He, Z., Zhou, J., and Firestone, M. K. (2016). The interconnected rhizosphere: high network complexity dominates rhizosphere assemblages. Ecol. Lett. 19, 926–936. doi: 10.1111/ele.12630
Sigler, W. V., and Zeyer, J. (2002). Microbial diversity and activity along the forefields of two receding glaciers. Microb. Ecol. 43, 397–407. doi: 10.1007/s00248-001-0045-5
Steele, J. A., Countway, P. D., Xia, L., Vigil, P. D., Beman, J. M., Kim, D. Y., et al. (2011). Marine bacterial, archaeal and protistan association networks reveal ecological linkages. ISME J. 5, 1414–1425. doi: 10.1038/ismej.2011.24
Taffner, J., Erlacher, A., Bragina, A., Berg, C., Moissl-Eichinger, C., and Berg, G. (2018). What is the role of Archaea in plants? new insights from the vegetation of alpine bogs. mSphere 3:e122-18. doi: 10.1128/mSphere.00122-18
Tanaka, Y., Tamaki, H., Tanaka, K., Tozawa, E., Matsuzawa, H., Toyama, T., et al. (2018). “Duckweed-Microbe Co-Cultivation Method” for Isolating a Wide Variety of Microbes Including Taxonomically Novel Microbes. Microb. Environ. 33, 402–406. doi: 10.1264/jsme2.ME18067
Tarlera, S., Jangid, K., Ivester, A. H., Whitman, W. B., and Williams, M. A. (2008). Microbial community succession and bacterial diversity in soils during 77 000 years of ecosystem development. FEMS Microbiol. Ecol. 64, 129–140. doi: 10.1111/j.1574-6941.2008.00444.x
Tscherko, D., Rustemeier, J., Richter, A., Wanek, W., and Kandeler, E. (2003). Functional diversity of the soil microflora in primary succession across two glacier forelands in the Central Alps. Eur. J. Soil Sci. 54, 685–696. doi: 10.1046/j.1351-0754.2003.0570.x
Tu, Q., Yuan, M., He, Z., Deng, Y., Xue, K., Wu, L., et al. (2015). Fungal communities respond to long-term CO2 elevation by community reassembly. Appl. Environ. Microbiol. 81, 2445–2454. doi: 10.1128/AEM.04040-14
Waldrop, M. P., Zak, D. R., Blackwood, C. B., Curtis, C. D., and Tilman, D. (2006). Resource availability controls fungal diversity across a plant diversity gradient. Ecol. Lett. 9, 1127–1135. doi: 10.1111/j.1461-0248.2006.00965.x
Wardle, D. A., Yeates, G. W., Watson, R. N., and Nicholson, K. S. (1995). Development of the decomposer food-web, trophic relationships, and ecosystem properties during a three-year primary succession in sawdust. Oikos 73, 155–166.
Welti, E. A. R., and Joern, A. (2015). Structure of trophic and mutualistic networks across broad environmental gradients. Ecol. Evol. 5, 326–334. doi: 10.1002/ece3.1371
Widder, S., Besemer, K., Singer, G. A., Ceola, S., Bertuzzo, E., Quince, C., et al. (2014). Fluvial network organization imprints on microbial co-occurrence networks. Proc. Natl. Acad. Sci. U.S.A. 111, 12799–12804. doi: 10.1073/pnas.1411723111
Williams, R. J., Howe, A., and Hofmockel, K. (2014). Demonstrating microbial co-occurrence pattern analyses within and between ecosystems. Front. Microbiol. 5:358. doi: 10.3389/fmicb.2014.00358
Yarwood, S. A., and Högberg, M. N. (2017). Soil bacteria and archaea change rapidly in the first century of Fennoscandian boreal forest development. Soil Biol. Biochem. 114, 160–167. doi: 10.1016/j.soilbio.2017.07.017
Yeates, G. W., Bongers, T., de Goede, R. G. M., Freckman, D. W., and Georgieva, S. S. (1993). Feeding habits in soil nematode families and genera - an outline for soil ecologists. J. Nematol. 25, 315–331.
Zhou, J., Xia, B., Treves, D. S., Wu, L. Y., Marsh, T. L., O’Neill, R. V., et al. (2002). Spatial and resource factors influencing high microbial diversity in soil. Appl. Environ. Microbiol. 68, 326–334. doi: 10.1128/aem.68.1.326-334.2002
Zhou, J. Z., Deng, Y., Luo, F., He, Z. L., and Yang, Y. F. (2011). Phylogenetic molecular ecological network of soil microbial communities in response to elevated CO2. mBio 2:e122-11. doi: 10.1128/mBio.00122-11
Keywords: bacteria, diversity, fungi, interaction network, joint distribution model, 16S, 18S, ITS
Citation: Farrer EC, Porazinska DL, Spasojevic MJ, King AJ, Bueno de Mesquita CP, Sartwell SA, Smith JG, White CT, Schmidt SK and Suding KN (2019) Soil Microbial Networks Shift Across a High-Elevation Successional Gradient. Front. Microbiol. 10:2887. doi: 10.3389/fmicb.2019.02887
Received: 19 August 2019; Accepted: 29 November 2019;
Published: 18 December 2019.
Edited by:
Manuel Delgado Baquerizo, University of Alicante, SpainReviewed by:
Jeff R. Powell, Western Sydney University, AustraliaHugo Saiz, Rey Juan Carlos University, Spain
Copyright © 2019 Farrer, Porazinska, Spasojevic, King, Bueno de Mesquita, Sartwell, Smith, White, Schmidt and Suding. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Emily C. Farrer, ZWZhcnJlckB0dWxhbmUuZWR1