Moira Hough1*
Amelia McClure2Benjamin Bolduc3Ellen Dorrepaal4Scott Saleska1
Vanja Klepac-Ceraj2*
Virginia Rich3*- 1Department of Ecology and Evolutionary Biology, University of Arizona, Tucson, AZ, United States
- 2Department of Biological Sciences, Wellesley College, Wellesley, MA, United States
- 3Department of Microbiology, College of Arts and Sciences, The Ohio State University, Columbus, OH, United States
- 4Climate Impacts Research Centre, Umeå University, Umeå, Sweden
Plant-associated microbiomes are structured by environmental conditions and plant associates, both of which are being altered by climate change. The future structure of plant microbiomes will depend on the, largely unknown, relative importance of each. This uncertainty is particularly relevant for arctic peatlands, which are undergoing large shifts in plant communities and soil microbiomes as permafrost thaws, and are potentially appreciable sources of climate change feedbacks due to their soil carbon (C) storage. We characterized phyllosphere and rhizosphere microbiomes of six plant species, and bulk peat, across a permafrost thaw progression (from intact permafrost, to partially- and fully-thawed stages) via 16S rRNA gene amplicon sequencing. We tested the hypothesis that the relative influence of biotic versus environmental filtering (the role of plant species versus thaw-defined habitat) in structuring microbial communities would differ among phyllosphere, rhizosphere, and bulk peat. Using both abundance- and phylogenetic-based approaches, we found that phyllosphere microbial composition was more strongly explained by plant associate, with little influence of habitat, whereas in the rhizosphere, plant and habitat had similar influence. Network-based community analyses showed that keystone taxa exhibited similar patterns with stronger responses to drivers. However, plant associates appeared to have a larger influence on organisms belonging to families associated with methane-cycling than the bulk community. Putative methanogens were more strongly influenced by plant than habitat in the rhizosphere, and in the phyllosphere putative methanotrophs were more strongly influenced by plant than was the community at large. We conclude that biotic effects can be stronger than environmental filtering, but their relative importance varies among microbial groups. For most microbes in this system, biotic filtering was stronger aboveground than belowground. However, for putative methane-cyclers, plant associations have a stronger influence on community composition than environment despite major hydrological changes with thaw. This suggests that plant successional dynamics may be as important as hydrological changes in determining microbial relevance to C-cycling climate feedbacks. By partitioning the degree that plant versus environmental filtering drives microbiome composition and function we can improve our ability to predict the consequences of warming for C-cycling in other arctic areas undergoing similar permafrost thaw transitions.
Introduction
As ecosystems warm and undergo transitions due to climate change, we are seeing significant shifts in microbial community ecology and function with important implications for carbon (C) storage (Zhou et al., 2012). This is especially true in arctic peatlands which are warming fast and undergoing major ecosystem transitions as permafrost thaws (McCalley et al., 2014; Mondav et al., 2017; Monteux et al., 2018). Microbial communities are strongly impacted by the environmental changes associated with permafrost thaw, particularly hydrology, available C-substrates, and plant communities. It has been well established that different stages of permafrost-thaw harbor substantially different bulk peat microbial communities, and these changes have major implications for C-cycling, particularly methane (Tveit et al., 2013; McCalley et al., 2014; Mondav et al., 2017; Monteux et al., 2018; Woodcroft et al., 2018a). However, plant-associated microbial communities remain under-studied. In particular, it is unclear whether plant microbiome changes are determined more by environmental changes as permafrost thaws, or by the biotic effects of plant community changes especially in the rhizosphere (Tkacz et al., 2015) and phyllosphere (Redford et al., 2010). Abiotic effects are often thought of as operating at a macro-scale as a filter that selects a potential species pool, upon which biotic effects act at the micro-scale (Aguilar et al., 2017; Cadotte and Tucker, 2017; Thakur and Wright, 2017). However, in the case of a permafrost-thaw front, both effects take place on a similar scale, with partial overlap of plant species (and therefore biotic effects) across a strong but geographically small environmental gradient. This raises the question of whether microbial community composition is more strongly determined by plant associate (i.e., biotic filtering) or thaw stage (i.e., environmental filtering), and whether plant selection on the microbial community has implications for C-cycling. This is especially important because permafrost is estimated to contain more than a third of the C in the top 1–3 m of the Earth’s soil and its thaw is expected to lead to C releases to the atmosphere in the range of 37–174 Pg by 2,100 which would increase climate warming by 0.13–0.27°C (Tarnocai et al., 2009; Hugelius et al., 2014; Schuur et al., 2015). Microbes serve as the primary decomposers and gatekeepers determining how much C will remain in the ecosystem versus being released to the atmosphere as CO2 and CH4. Indeed, in some systems changes in microbial community structure and function have been shown to mitigate C-release predicted due to warming alone (Zhou et al., 2012). Therefore accurate prediction of the magnitude of C releases from the arctic requires understanding how microbial communities will be impacted by permafrost thaw.
While thawing permafrost peat microbial communities have been well studied (Tveit et al., 2013; McCalley et al., 2014; Mondav et al., 2017; Monteux et al., 2018; Woodcroft et al., 2018b), the plant-associated microbial communities in these systems have not been well characterized, nor is it clear to what extent plant associations are responsible for the changes in bulk peat microbial community. In most systems, rhizosphere and phyllosphere microbial communities seem to draw at least some members from the communities present in the surrounding environment but they form distinct communities based on factors such as the species, genotype, and health of the host-plant as well as the compartment in question (Delmotte et al., 2009; Redford et al., 2010; Fonseca et al., 2016; Nuccio et al., 2016). For instance, the phyllosphere is likely to be enriched in alpha-proteobacteria, particularly those adapted to stress and with either one-C central metabolism (e.g., methanotrophy and methylotrophy) such as Methylobacterium, or a diverse range of substrate use such as Sphingomonas (Delmotte et al., 2009; Knief et al., 2012). Since both metabolic types have been found on the same phyllosphere samples, this divergence in metabolic strategy may indicate that some organisms capitalize on plant-associated substrates whereas others remain generalists in a resource-limited environment. Controls on rhizosphere community structure are complex, including the structure of the roots, the density, size of their filaments, and on fungal associations of the host plant (Parniske, 2008). These rhizosphere-associated communities are often less diverse than the microbial communities in the surrounding soil matrix but more diverse than the phyllosphere (Berg and Smalla, 2009; Knief et al., 2012). Rhizosphere communities are often enriched in alpha- and beta-proteobacteria, and depleted in Actinobacteria compared to soil microbial communities, indicating plant influence on community structure and composition (Knief et al., 2012; Nuccio et al., 2016).
Plant selection on leaf and root microbiomes is potentially important in mediating microbial controls on decomposition (Strickland et al., 2009) and ultimately C flux to the atmosphere. The first communities that have the opportunity to decompose plant material are the microbes colonizing the plant itself (Tláskal et al., 2016). Plants often have specialized phyllosphere communities which accompany senescing leaf tissue into the soil (Redford et al., 2010; Nissinen et al., 2012; Bragina et al., 2015; Lebeis, 2015; Tláskal et al., 2016; Baldrian, 2017). Members of the phyllosphere have been shown to participate in the initial decomposition of fresh organic material before being replaced by microbes from the bulk soil community, and leaf litter composition can influence decomposition dynamics such as fungal to bacterial ratios (Kögel-Knabner, 2002; Thiessen et al., 2013; Tláskal et al., 2016; Baldrian, 2017). Additionally, plant root exudates are known to be a source of C to soil and to determine rhizosphere microbial community composition (Tkacz et al., 2015). Microbial response to differences in C quality in root exudates (such as through shifts in activity or bacterial to fungal ratios) can influence decomposition dynamics such as by stimulating decomposition or causing additional C release from bulk peat through priming effects (Cheng et al., 2014; Pegoraro et al., 2019). Permafrost thaw has been shown to dramatically alter root growth patterns including quantity and quality of litter inputs which can be expected to influence microbial community composition (Blume-Werry et al., 2019). However, research characterizing permafrost and bulk peat microbial communities in arctic systems (Mackelprang et al., 2011; Hultman et al., 2015; Mondav et al., 2017; Woodcroft et al., 2018a), has yet to clearly describe the role of plants in influencing the microbial community composition, decomposition, and C cycling.
Since permafrost thaw results in changes to both environmental and plant filtering effects on the microbial community, we investigated their relative importance in structuring rhizosphere versus phyllosphere microbial communities, at a well-studied permafrost thaw gradient in Stordalen Mire, Sweden. In this area, as permafrost thaws, palsa peat mounds (which rise above the surrounding wet areas due to a supportive core of ice) collapse to form waterlogged bog areas which initially remain hydrologically disconnected from the later thaw stages but eventually sink and form fully inundated, hydrologically connected, nutrient-rich fens. As these three soil environments differ more strongly than the air, we hypothesized that biotic filtering would be more important than environmental filtering in the phyllosphere, but the reverse would be true in the rhizosphere. Microbial communities are highly diverse and often include taxa which may be more important to community structure than expected based on abundance measures (keystone species) as well as many taxa that are of little importance to community structure. Therefore, we consider two metrics by which microbial community differences may be measured across sites, compartments (rhizosphere, phyllosphere, or peat), and plant associates: (1) differences in relative abundance and (2) differences in importance to network structure (or keystoneness). Furthermore, to investigate whether these community composition differences could be important to the ecosystem’s C-cycling, we investigated functional groups important to CH4 cycling (as a particularly important C-cycling function which is also phylogenetically constrained). This approach allows us to identify the degree to which plant and environmental filtering influence C cycling through changes in microbial functional groups as well as overall microbial community composition, structure, and putative function.
Materials and Methods
Site and Sample Collection
Samples were collected at three time points through the 2015 growing season: early season (20 and 21 June), peak growing season (20 and 23 July), and late growing season (2 September). Stordalen Mire is located 10 km east of Abisko, at 68° 219 N, 18° 499 E, and is 363 m above sea level (with ecologically relevant microtopography across the Mire spanning several meters’ elevation). The site is managed by the Abisko Scientific Research Station, the University of Stockholm, and the Integrated Carbon Observation System. Within the mire, there are three main habitats spanning a permafrost thaw gradient (Figure 1): no thaw (palsa type), initial partial thaw (bog type), and complete thaw (fen type). Palsas consist of raised, permafrost underlain areas characterized by low shrubs (e.g., Betula nana, Empetrum nigrum, Andromeda polifolia, Vaccinium spp.), forbs (e.g., Rubus chamaemorus), graminoids (e.g., Eriophorum vaginatum), lichens, and drier mosses (e.g., Sphagnum fuscum). Bogs are wetter low-lying areas often still underlain by permafrost lenses characterized by more hydric species of sphagnum (e.g., Sphagnum balticum) and small sedges (e.g., E. vaginatum and Carex rotundata). Fens are formed after complete permafrost thaw and collapse. They are characterized by standing water, and hydric sphagnum and sedge species (e.g., Eriphorum angustifolium, C. rotundata, and Carex rostrata). As of 2010, 49% of the area within Stordalen was made up of intact palsa, 37% was made up of partially thawed bog, 12% was made up of fully thawed fen (Bäckstrand et al., 2010). From 1970 to 2000, palsa extent shrank while bog sites expanded by 3% and fen sites expanding by 54% (Bäckstrand et al., 2010).
Figure 1. Permafrost thaw progression showing directions of possible microbial dispersal in the context of changes in plant community composition and hydrology as ecosystems transition from: I. permafrost-underlain palsa, to II. partially thawed bog, to III. fully thawed and inundated fen. The underlying thaw progression figure has been adopted from Johansson et al. (2006) and modified substantially.
Three areas were selected where all three habitats could be found within close (20 m2) proximity. Samples of leaf, root, and associated peat for seven abundant, habitat-typical plant species were selected for phyllosphere and rhizosphere characterization across the mire’s three habitats. Species representative of palsa areas were: R. chamaemorus, A. polifolia, E. vaginatum, and Sphagnum sp. Species representative of bog areas were: A. polifolia, E. vaginatum, C. rotundata, and Sphagnum sp. Species representative of fen areas were: E. angustifolium, and C. rotundata. Note that palsa and bog Sphagnum spp. differed, and that several species were representative in more than one habitat (and were sampled separately in each habitat in case of physiological differences). Fewer species were selected to represent the fen habitat based on its much lower plant diversity (Johansson et al., 2006). Replicate plants were sampled in two ways: spatially distributed and spatially proximal. For the former, one individual of each habitat-representative species was sampled at each of the three locations (generating three spatially distributed replicates). For the latter, at one location, two additional sets of plant samples were collected for each habitat-representative species (generating a total of three spatially proximal replicates). In addition, a single peat sample was taken at each location, in each habitat (for three spatially distributed replicates per habitat).
All material was collected using forceps sterilized with ethyl alcohol and rinsed twice with water. For each plant, phyllosphere samples were collected first, consisting of approximately 2 g of living Sphagnum spp. tissue, all leaf tissue from A. polifolia and R. chamaemorus plants, and 1 cm2 of sedge plants (Eriophorum spp. and C. rotundata). When all phyllosphere sample material had been collected, a serrated edge knife was used to cut around the base of the plant. The plant was then pulled out with associated peat still attached. “Bulk peat” samples were collected as 15 mL of peat from 2 to 3 cm deep, from as close to the root as possible without touching the root. Rhizosphere samples were collected by snipping small pieces of the root fibers into a 1.5 mL Eppendorf tube. Both the primary root and lateral roots were included and the root tip was excluded. A total of 4 cm of main and lateral roots were sampled from A. polifolia and R. chamaemorus plants, and one segment of the main root node and 2–4 fine root hairs were sampled for each sedge rhizosphere sample. While Sphagnum spp. plants do not have roots, 2 g samples were taken from the base of the stalk where the tissue was brown but still attached to the living stalk. After collection, the samples were stored in a cooler, and then transferred to a −80°C freezer within 4 h of collection. The samples were shipped on dry ice to Wellesley College and immediately stored at −80°C until further processing.
DNA Extraction and Sequencing of 16S rRNA Gene
For each sample a standardized quantity of 50 mg of Sphagnum spp. photosynthetic and non-photosynthetic tissue, 75 mg of angiosperm leaf tissue, 100 mg of angiosperm root material, or 50 mg of peat was used for DNA extraction. In all cases bulk samples were used, meaning that results include both endo- and epiphytic microbial communities. Total DNA from each sample was extracted using the PowerSoil® DNA Isolation Kit (MoBio, Carlsbad, CA, United States) according to the manufacturer’s instructions and eluted in 50 μL of elution buffer. After extraction, DNA was quantified using NanoDrop (Thermo Scientific, Inc., Wilmington, DE, United States) and stored at −20°C.
Sequencing at Argonne National Laboratories followed the protocols of Caporaso et al. (2012). PCR amplification of the V4 region of 16S rRNA gene made use of the forward primer 515F (5′-GTGCCAGCMGCCGCGGTAA-3′) and the reverse primer 806R (5′-GGACTACHVGGGTWTCTAAT-3′). PCR amplicons were quantified by PicoGreen (Invitrogen, Carlsbad, CA, United States) using a plate reader, and amplicons were pooled in equal concentrations into a single tube. This pool was cleaned up using UltraClean PCR Clean-Up Kit (MoBio, Carlsbad, CA, United States) and quantified using the Qubit (Invitrogen, Carlsbad, CA, United States). The pooled samples were sequenced on the Illumina MiSeq platform (Illumina, San Diego, CA, United States) according to the sequencing procedures described in Caporaso et al. (2012). Sequences are available in the National Center for Biotechnology Information Sequence Read Archive, Study Accession: PRJNA599435.
Bioinformatic Analysis
An in-house workflow was used to process operational taxonomic unit (OTU) data, which was based primarily on Quantitative Insights into Microbial Ecology (Caporaso et al., 2010) and SparCC. Detailed information about the SparCC workflow can be found in Friedman and Alm (2012). Briefly, paired-end sequences were first joined, then split into libraries according to their index. Sequences were then filtered according to quality using USEARCH fastq_filter (-fastq_maxee 0.5), dereplicated and screened against phiX174. Surviving sequences were clustered into OTUs at 97% average nucleotide identity, then taxonomy assigned using the Greengenes database (DeSantis et al., 2006) version 13_8 using the Ribosomal Data Project classifier (McDonald et al., 2012). Finally, OTUs matching chloroplast and mitochondrial sequences were removed from the analysis. Unifrac analysis was performed in QIIME (Lozupone and Knight, 2005). OTU tables were pruned to address specific experimental questions, then LEfSe biomarker analysis was completed (Segata et al., 2011).
Network analysis began with filtering of OTU tables to exclude any OTUs with fewer than 10 counts and represented in fewer than three total samples across all habitat, tissue or plant types. To generate networks, OTU tables were split according to habitat, tissue type, and plant type, and mean OTU abundances for each sample type were calculated based on filtered tables of counts as described above. Compositionality-robust correlations of each OTU table (15 total) was generated using SparCC using the median of 20 iterations. Data were resampled and bootstrapped 100 times (20 iterations) to generate correlation p-values. Correlation matrices were then processed using python packages Numpy (Walt et al., 2011), Scipy (Virtanen et al., 2020), and Pandas (McKinney, 2010). Each correlation matrix has its correlations filtered according to its p-value (<0.05), and networks were built using the filtered correlations as an adjacency matrix and rendered using python-igraph (Csardi and Nepusz, 2006) with visual refinement using Adobe Illustrator CC such that OTUs served as “nodes” or vertices, and absolute correlation values as the “edges” between OTUs.
We interpreted OTU PageRank score, which quantifies connectivity, as a metric of ecological importance. We considered the most important OTUs to be those in the 99th percentile for each network based on PageRank scores, and compared interpretation based on PageRank to the additional network metrics of degree, transitivity, betweenness, and closeness. Demonstration of the strong correlation of PageRank with the “consequentiality” of nodes in diverse settings: for identification of key metabolic genes in pathways of co-occurrence gene networks (Browne et al., 2016; Wang et al., 2017, extinction cascades in food webs (Allesina and Pascual, 2009), and complex phylogenetic signals in communities of co-occurring bacteria (Estrada-Pena et al., 2018). This metric has been successfully applied to extinction food webs (Allesina and Pascual, 2009). The nature of the metric (which quantifies the connections of a node to others, with weighting for those neighbors’ respective connectivity too), and these other studies, supported its use as an indicator of keystone organisms. In order to relate OTU importance and abundance (which are on different scales), we used the rank of each within every sample, rather than its absolute value.
Microbial functional groups of interest were methanogens and methanotrophs. While functional assignation of organisms identified using 16S sequencing cannot be definitive, methane cycling traits are generally highly conserved (Martiny et al., 2015) meaning that organisms belonging to taxonomic groups known to produce or consume methane are likely to have similar lifestyles. Methanogens were identified based on classification into one of the following orders known to be dominated by methanogens and active at our site: Methanosarcinales, Methanobacteriales, Methanomicrobiales (Mondav et al., 2014; Woodcroft et al., 2018b; Evans et al., 2019). Methanotrophs were identified based on classification into one of the following groups known to be predominantly methanotrophic and active at our site: order Methylacidiphilales, families Beijerinckiaceae, Hyphomicrobiaceae, Methylocystaceae, Crenotrichaceae, and Methylococcaceae (Tveit et al., 2013; Knief, 2015; Singleton et al., 2018; Smith and Wrighton, 2019).
Statistical Analysis
NMDS and variance partitioning were performed on both Unifrac and Bray–Curtis distances of OTU tables for all samples. Significance tests were performed with redundancy analysis for Unifrac data and distance-based redundancy analysis for Bray–Curtis distances using the vegan package in R (Oksanen et al., 2019). Hierarchical clustering (based on Jaccard and Bray–Curtis distances) was performed on both abundance and importance scores of OTUs in networks. All statistical analyses were performed in R studio running R version 3.6.0 (R Studio Team, 2018; R Core Team, 2019).
Results
A total of 350 samples were sequenced and rarefied to include 1,000 high-quality bacterial or archaeal 16S rRNA gene sequences per sample. Alpha diversity varied between rarefaction analysis indicating that full sampling depth was achieved only for samples with lower levels of diversity, and not with samples from the fully thawed fen peat (Supplementary Figure S1). Out of all sequences, 99.6% were assigned to Bacteria and 0.4% to Archaea – the latter is consistent with previous work from this site (Mondav et al., 2017) and reflects the known biases of the V4 primer set against Archaea (Parada et al., 2016). Across all samples, there were a total of 373 unique OTUs (1,436 in the unrarefied dataset) which were assigned to 170 genera and 31 known phyla (170 genera and 32 phyla in the unrarefied dataset). Diversity within a sample ranged from 124 to 361 OTUs (181 to 609 in the unrarefied dataset). The four most dominant phyla when considering the entire dataset were Proteobacteria (44%), Acidobacteria (22.9%), Bacteroidetes (8.9%), and Verrucomicrobia (7.8%). This profile matched the general profiles found in previous investigation into bulk peat communities at this site by Mondav et al. (2017) but with stronger and more consistent dominance of Proteobacteria and lower relative abundance of Acidobacteria across all habitat types, particularly in the phyllosphere. This is generally in line with the findings of Mondav et al. (2017) despite major differences in methodology.
Diversity and Composition of Bacterial Taxa Across Sample Types
A total of 66 OTUs were ubiquitous across all sample types (habitat, compartment, and plant) after removal of chloroplast sequences. These belonged to 10 phyla: Acidobacteria, Actinobacteria, Armatimonadetes, Bacteroidetes, Chlamydiae, Cyanobacteria, Planctomycetes, Proteobacteria, Verrucomicrobia, and WPS-2. Bacterial composition at the phylum level showed the same major phyla represented across the permafrost thaw gradient with variation in relative abundance by habitat type and compartment (Figure 2). Based on LEfSe analysis, we identified phyla that were enriched in abundance in one habitat as compared to the others within each compartment (LDA > 2). The palsa was enriched in Acidobacteria and Planctomycetes in all three compartments, as well as Actinobacteria and Chlamydiae in the rhizosphere and peat, and WPS-2 and FCPU426 in the rhizosphere only. The bog was enriched in WS4 in all compartments, TM6 in the phyllosphere and rhizosphere, Armatimonadetes and Verrucomicrobia in the rhizosphere and peat, OD1 in the phyllosphere only, FBP and Proteobacteria in the rhizosphere only, and WPS-2 in the peat only. The fen contained the greatest number of enriched phyla with Bacteroidetes, Chlorobi, Chloroflexi, Fibrobacteres, Firmicutes, Nitrospirae, Spirochaetes, Euryarchaeota, and Parvarchaeota enriched in all three compartments. The fen rhizosphere and peat were also enriched in Crenarchaeota, Elusimicrobia, OP3, and WWE1, with the rhizosphere enriched for OD1, TM7, AD3, Cyanobacteria, and Fusobacteria, and the peat enriched in Gemmatimonadetes (Supplementary Figure S2).
Figure 2. Changes in relative abundance of bacterial phyla across habitats (palsa = brown, bog = green, and fen = blue) and compartments (phyllosphere, rhizosphere, and bulk peat). Circle size corresponds to percent relative abundance.
NMDS based on Bray–Curtis distances between OTUs showed separation of centroids by habitat with only partial separation by compartment and plant type and very little separation by season (Figure 3). PERMANOVA with plant and habitat nested within compartment indicated significant differences (p < 0.001) between communities across all factors (compartment, plant, habitat, and month) as well as significant interaction terms (p < 0.001). Variance partitioning indicated that when considered individually, habitat and plant accounted for the majority of the variation with compartment playing a smaller role and sample month having almost no impact (Table 1). However, the most variance was explained when habitat, plant associate, and compartment were all considered (adjusted r2 = 0.27, p < 0.001). Repeating the analyses using phylogenetically weighted Unifrac values to take phylogenetic similarity into account provided similar results but with an increase in the unique amounts of variance explained, particularly by compartment (adjusted r2 = 0.39, p < 0.001, Table 1).
Figure 3. NMDS plots of Bray–Curtis distances between samples from each (A) habitat, (B) compartment, (C) plant associate, and (D) month. Outer ellipses include 95% of points for a habitat, inner ellipses are 95% confidence around the centroid. Stress = 0.17.
Table 1. Variance explained by each factor based on Bray–Curtis and Unifrac distances of total dataset and Bray–Curtis distances of only OTUs included in networks.
Considering each compartment separately, NMDS plots indicated some separation by habitat type and plant associate, especially for phyllosphere and rhizosphere samples (Figure 4). Variance partitioning based on Bray–Curtis distance between communities showed different patterns for each compartment. In the rhizosphere, plant associate explained a similar to the amount of variance as habitat, but in the phyllosphere plant associate explained substantially more variance than habitat (Table 2). Plant associate explained a similar amount of variance in both rhizosphere and phyllosphere. However, habitat explained substantially more variance in the rhizosphere and the bulk peat than the phyllosphere. The total amount of variance explained by plant and habitat together was much higher in the rhizosphere and bulk peat than in the phyllosphere (Table 2). Much of the explained variance overlapped between both habitat and plant type, likely because plant species distribution was not even across all habitats.
Figure 4. NMDS plots of Bray–Curtis distances separated by compartment and color-coded by (A) habitat type and (B) plant associate. Outer ellipses include 95% of points for a group, inner ellipses are 95% confidence around the centroid. Stress = 0.17.
Table 2. Variance partitioning by compartment based on Bray–Curtis distance between communities considering each compartment (phyllosphere, rhizosphere, and bulk peat) separately.
Methane-Cycling Functional Groups
Putative methanotrophs were identified as bacteria belonging to groups which have previously been shown to have active methanotrophic representatives at this site based on metagenomic and metatranscriptomic data (Tveit et al., 2013; Knief, 2015; Singleton et al., 2018; Smith and Wrighton, 2019). These were order Methylacidiphilales and families: Beijerinckiaceae, Hyphomicrobiaceae, Methylocystaceae, Crenotrichaceae, and Methylococcaceae. Methanotrophs, particularly members of Methylocystaceae, were found across all habitats, compartments, and plants though in varying abundances (Figure 5A and Supplementary Table S3). Members of Methylacidiphilales were common in palsa and bog habitats but rare in the fen, whereas members of Methylococcales were common in the fen but rare in the bog and not seen in the palsa. Variance partitioning on putative methanotrophs showed similar patterns to the full dataset except that habitat explained less and plant associate more of the variance in phyllosphere putative methanotroph distribution than in the full dataset (Table 3).
Figure 5. Relative abundances of (A) Methanotrophs and (B) Methanogens, grouped by habitat and compartment. Colors indicate plant associates.
Table 3. Variance partitioning of putative methanotrophs by compartment based on Bray–Curtis distance between communities considering each compartment (phyllosphere, rhizosphere, and bulk peat) separately.
Putative methanogens were defined as Euryarchaeota in orders Methanosarcinales, Methanomicrobiales, Methanococcales, Methanobacteriales, Methanopyrales, Methanocellales, Methanomassiliicoccales, and “Candidatus Methanophagales” (Evans et al., 2019). Methanogens were found across all habitats and compartments as well as all plants except R. chamaemorus. Members of the family Methanobacteriaceae reached the highest abundances, especially in the fen, and were the most widespread (Figure 5B and Supplementary Table S4). Variance partitioning based on relative abundances did not match the pattern seen in the overall dataset. None of the phyllosphere variance was well explained. In the rhizosphere, plant associate explained substantially more variance than habitat, though overlap was high (Table 4).
Table 4. Variance partitioning of putative methanogens by compartment based on Bray–Curtis distance between communities considering each compartment (phyllosphere, rhizosphere, and bulk peat) separately.
Network Analysis
Because networks were based on the un-rarefied dataset, they contained greater diversity with a total of 1,436 OTUs. We avoided rarefication in order to preserve as many members of the community as possible; focusing on those whose importance was greater than their relative abundance. Rarefying data masks this signal. Since PageRank scores community members more highly when their surrounding community members are highly scored, rarefying even a single members’ abundance would diminish the importance of those to which it is connected. The mean diversity per network type was 361 with a minimum of 181 and a maximum of 609. Variance partitioning on the OTUs included in the networks recapitulated the results seen in the entire dataset with slight increases in the amount of variance explained by plant and habitat, especially when using only keystone OTUs (Table 1). Network structure varied substantially by plant associate, retaining some similarity across habitat types for a given plant species (Figure 6). Sphagnum, in particular, showed a characteristic “donut” structure regardless of habitat – despite these being different species in the different habitats. This structure reflects a lack of a central hub of OTUs, with instead dispersed interconnected hubs. This may indicate the presence of multiple sub-communities connected by a few interacting species, which could be caused by distinct niches. A similar structure was found in a few other networks such as E. vaginatum in the bog, but was not consistent across habitats. Similarities in network structure between tissue types of a given plant species are not meaningful since the shape of the network structure visualization was held constant for each plant species in each habitat so only node size and connectivity varied across tissue types.
Figure 6. Network diagrams of the microbial community associated with each sample type (separated by habitat, compartment, and plant associate). Blue lines indicate positive correlations and red lines indicate negative correlations. Colors indicate phyla.
To further test whether microbial communities were more similar when sampled from the same habitat versus the same plant associate, we used hierarchical clustering using Bray–Curtis distances based on both relative abundance and network importance score (PageRank or keystoneness). The phyllosphere microbial communities mostly clustered by habitat rather than plant type with the exception of Sphagnum which clustered on its own based on importance to network but with the bog plants based on relative abundance (Figures 7A,B). Rhizosphere microbial communities clustered primarily by habitat when importance to network was considered, however, when relative abundance was considered sedges from the bog clustered with the fen sedges rather than with the other bog plants (Figures 7C,D). As would be expected, bulk peat microbial communities clustered by habitat based on both relative abundance and importance to network (Figures 7E,F). Hierarchical clustering based on Jaccard dissimilarity which considers only presence/absence not relative abundance showed no differences in the structure of groupings but resulted in increased dissimilarity between groups (Supplementary Figure S3).
Figure 7. Hierarchical clustering of Bray–Curtis dissimilarity between community networks weighting organisms by (A,C,E) relative abundance for a given habitat/plant and (B,D,F) relative importance for a given habitat/plant. Main colors correspond to habitat with brown = palsa, green = bog, blue = fen. Accent colors represent plant species with orange = A. polifolia, brown = R. chamaemorus, olive = E. vaginatum, dark green = Sphagnum, teal = C. rotundata, dark blue = E. angustifolium.
Keystone Organisms
Keystone organisms were defined as OTUs in the 99th percentile of importance for each network based on each organism’s PageRank score. PageRank has not typically been used as a metric for identification of keystone organisms in microbial communities. However, it has been shown to be strongly correlated with extinction cascades in food webs (Allesina and Pascual, 2009), identification of key metabolic genes in pathways of co-occurrence gene networks (Browne et al., 2016; Wang et al., 2017), and used to disentangle phylogenetic signals in communities of co-occurring bacteria (Estrada-Pena et al., 2018) thus indicating that it may indeed be a strong indicator of keystone organisms. PageRank score within a network was strongly correlated with betweenness centrality (Supplementary Figures S4, S5), which has previously been used as a metric for identifying keystone organisms (Martín González et al., 2010). Keystone OTUs also scored highly for closeness centrality, betweenness centrality, degree, and transitivity (Supplementary Figure S4). A total of 75 unique OTUs were identified as keystone organisms across the 27 different networks (Supplementary Table S1). This represented an average of four keystone taxa per network with a minimum of two and a maximum of seven (Supplementary Table S2). OTUs were generally only keystone for one or a few networks and none were keystones across all networks, habitats, compartments, or plants. Considering relative abundances of keystone OTUs across all samples indicated that they ranged from generalists (found across habitats, plants, and compartments) to specialists at any of the three levels (by habitat, plant type, and compartment) (Figure 8). Generalists included members of families Sinobacteraceae (OTU 12), Acetobacteraceae (OTU 137), and Acidobacteriaceae (OTU 274). Specialists included members of an unknown family in Chloroflexi (OTU 264) and Nakamurellaceae (OTU 737) found only in the fen, a member of an unknown family of Acidobacteria (OTU 90) found only in palsa, a member of Neisseriaceae (OTU 5321) found only in the rhizosphere, and members of Methylocystaceae (OTU 167) and Cytophagaceae (OTU 371) found only in phyllosphere. Variance partitioning based on relative abundances of the keystone organisms increased the total amount of variation explained to 33% but did not change the overall patterns in variance (Table 1). Six keystone OTUs belonged to functional groups of interest, one methanogen: OTU 36 (family Methanobacteriaceae) in the bog peat near a Sphagnum sp.; and five putative methanotrophs: OTU 22 (family Methylocystaceae) associated with A. polifolia rhizosphere in the bog and with peat near C. rotundata in the fen, OTU 40 (family Bradyrhizobiaceae) associated with the rhizosphere of R. chamaemorus in the palsa and E. angustifolium in the fen, OTU 72 (order Methylacidiphilales) associated with the rhizosphere of A. polifolia in the palsa and C. rotundata in the bog, OTU 167 (family Methylocystaceae) associated with the phyllosphere of A. polifolia and E. vaginatum in the palsa, OTU 171 (order Methylococcales) associated with the rhizosphere of A. polifolia, E. vaginatum, and the basal area of Sphagnum sp. in the bog, and OTU 320 (family Crenotrichaceae) associated with the bulk peat near E. angustifolium in the fen (Supplementary Table S1).
Figure 8. Comparison of (A) relative importance and (B) relative abundance of organisms in the 99th percentile of importance to network structure across habitat, compartments, and plant associates. Stars indicate methanotrophs and diamonds indicate methanogens.
Discussion
We investigated whether biotic filtering by plant associate could outweigh the importance of environmental filtering by habitat on microbial community composition across a permafrost thaw gradient. We found that the metric used for assessing differences between communities (importance to network versus relative abundance) has an impact on which taxa appear relevant. Under both approaches, we found evidence that in the rhizosphere, the environment first filters microbial community (especially through moisture), and then biotic interactions with plants further shape composition. In the phyllosphere, environmental conditions are less variable spatially so the effects of environmental filtering are less than the effects biotic influence on microbial community – though both are fairly small. However, we also found that the role of plant associate and habitat in shaping microbial communities may be very different for different microbial groups. Putative methanogens and methanotrophs appeared to be particularly strongly influenced by plant-associations, indicating that changes in plant community composition could have a significant influence on CH4 production from thawing permafrost ecosystems.
Below, we consider the following questions raised by our observations: What are the characteristics of the plant-associated microbiome as compared to the peat microbiome observed in other studies? Why did this study find relatively low influence of biotic filtering by plant associate on microbial community composition compared to other research? Does focusing analysis on keystone organisms that determine network structure provide clearer indication of the effects of filtering? Does plant community transition impact CH4 cycling in a way that contributes to permafrost thaw feedbacks?
Characteristics of the Plant-Influenced Microbiome
Microbial communities of Stordalen mire have been well studied but previous work has focused on bulk peat samples especially at depth (Mondav et al., 2017; Monteux et al., 2018; Woodcroft et al., 2018a). This represents the first known examination of plant-associated microbial communities at this site. Our most dominant phyla matched those seen in previous work in which Proteobacteria and Acidobacteria were the most common with Bacteroidetes and Actinobacteria also common (Mondav et al., 2017). This profile matched the general profiles found in previous investigation into bulk peat communities at this site by Mondav et al. (2017) but with stronger and more consistent dominance of Proteobacteria and lower relative abundance of Acidobacteria across all habitat types, particularly in the phyllosphere. This is consistent with previous work that has shown enrichment of alpha- and beta-proteobacteria in plant-associated communities – particularly the rhizosphere – as compared to bulk soil (Knief et al., 2012; Nuccio et al., 2016). While Mondav et al. (2017) found high relative abundance of Euryarchaeota in fen samples, we saw very low numbers with the highest abundance in fen bulk peat. We found four phyla in much higher abundance than Mondav et al. (2017) (in which they were absent or comprised less than 1% of the community): Firmicutes, Cyanobacteria, Armatimonadetes, and Spirochaetes. In contrast, we saw lower than 1% abundance of the following phyla that were reported at high abundance by Mondav et al. (2017): Microgenomates, Woesearchaeota, Thaumarchaeota, and Bathyarchaeota. These differences were strongest in the plant-associated microbial communities. Thus while our overall community composition was similar to that found in examinations of the deeper soil profile, our results are consistent with the idea that there are plant-mediated differences in microbial communities in the surface soil.
Our study also found substantially lower diversity than Mondav et al. (2017), who found a mean OTU richness of 721 with a range of 309–1,226. This may be because they report total diversity in pooled samples as opposed to our singletons which may have allowed resulted in the discovery of more different taxa. Additionally, our study used amplification of the V4 region and Illumina MySeq while Mondav et al. (2017) amplified the V6–V8 region and used 454 sequencing. Different regions of the 16S gene are known to provide different degrees of taxonomic resolution (Bukin et al., 2019) as well as biases in taxonomic composition (Sinclair et al., 2015; Fischer et al., 2016). Finally, our sampling approach for bulk peat was substantially different. Because we were sampling in the rooting zone close to plants, it is possible that even our bulk peat samples were more strongly plant influenced than those of Mondav et al. (2017), which were targeted at the peat profile. Thus both sequencing differences and plant influence may account for some of the differences in diversity and community composition between our bulk peat samples and those of Mondav et al. (2017).
Biotic Versus Environmental Filtering
Variance partitioning of relative abundances of the microbial community provided support for our hypothesis that habitat is less important for structuring phyllosphere microbial communities. In the phyllosphere, plant associate explained more variation than habitat – particularly when overlapping factors were excluded – indicating that biotic filtering was stronger than environmental filtering (Table 2). Additionally, both plant associate and habitat explained more variance in community composition in the rhizosphere than the phyllosphere. This is consistent with there being higher dispersal rates and less filtering (biotic or abiotic) on phyllosphere communities. Indeed, habitat and plant associate together only explained 20% of the variation in phyllosphere community composition (as compared to 35% in the rhizosphere and 33% in the bulk peat). This is consistent with previous findings that the phyllosphere tends to be more variable than the rhizosphere (Lebeis, 2015). We hypothesize that phyllosphere microbial communities are either more strongly determined by stochastic processes or simply more strongly impacted by variables not associated with plant-associate or thaw stage such as humidity or radiation. Testing these hypotheses was beyond the scope of the current study but would be a useful step for improving our understanding of phyllosphere microbial communities.
In contrast, in the rhizosphere plant associate and habitat explained similar amounts of variance in community composition (Table 2), indicating that both were important but neither strongly overrode the other. This was not entirely in line with our hypothesis which proposed that environmental filtering by habitat would override biotic effects in the rhizosphere. Indeed, biotic effects explained slightly more variance than habitat (8% vs 7%). This is particularly striking given the large differences in habitat type, ranging from dry palsa to fully inundated fen. This difference was large enough that our plant species did not span all three habitats, and some were found in only a single habitat. This last is a major cause of the high degree of overlap in the variance explained by both plant and habitat. While it has long been known that plants strongly influence their rhizospheres (Berg and Smalla, 2009; Knief et al., 2012; Ferrera-Rodríguez et al., 2013; Philippot et al., 2013; Lebeis, 2015; Lee et al., 2015), this influence is not always sufficient to override environmental differences (Kumar et al., 2016, 2017; Nuccio et al., 2016). Thus, it is striking that even with such a large difference between habitat types, plant associate still explained as much or slightly more variation than habitat.
One hypothesis to explain why previous studies have found that environmental factors were more important than species associate – in contrast to our finding and despite the strong influence of plants on rhizosphere communities – is that differences between plant functional groups are a key factor. Nuccio et al. (2016) found that rhizosphere community composition was most strongly correlated with moisture and then with region but not with species. However, this study was conducted within grasslands and focusing on a single plant genus (Avena). While our study included a strong moisture gradient, it also included a number of different plant functional types (sphagnum mosses, sedges, forbs, and ericaceous shrubs). It may be that the larger differences between plant functional types explain why plant associate had such a strong influence on rhizosphere community across our large environmental gradient. Indeed, research has shown substantial differences between plant species in quantity and quality of root litter inputs in thawing permafrost ecosystems which could influence microbial community composition (Blume-Werry et al., 2019). Cluster analysis based on relative abundances showed that sedge rhizosphere communities from the wet bog and fen habitats grouped together – though the sedge-associated community in the dry palsa grouped with its habitat (Figure 8). This is consistent with the finding of Nuccio et al. (2016) that moisture is a key factor filtering rhizosphere communities. Similarly, Sphagnum phyllosphere communities from palsa and bog habitats were more similar to each other than to any other communities, particularly when clustered by organisms’ network importance (Figure 8). Sphagnum mosses are known to cultivate a specialized endophytic microbial community (Bragina et al., 2012, 2014; Bragina et al., 2015; Holland-Moritz et al., 2018), however, this community has also been known to demonstrate high specificity to host and environmental factors. This means it is somewhat surprising that the phyllosphere communities from the different Sphagnum species and very different hydrologic conditions of the palsa and bog were so similar. Overall, these findings are consistent with the idea that environment first filters microbial community (especially through moisture), then biotic interactions with plants further shape composition.
Keystone Organisms
The high diversity of microbial communities provides a challenge in considering how to weight each taxonomic group in terms of its importance to community composition. Typically, community composition is compared using relative abundance as a weighting, which essentially assumes that the higher abundance of an organism, the larger is its importance to the community composition. While this may be true in many cases, when considering the impacts of community changes on ecosystem function often low-abundance organisms can be highly influential (Paine, 1969; Shade et al., 2014; Banerjee et al., 2018). Network analysis provides a mechanism for identifying these keystone organisms in microbial communities. Previous studies have used high betweenness and centrality scores to identify keystone taxa within microbial communities, however, there is not yet a consensus over which network metrics best measure keystoneness (Martín González et al., 2010; Berry and Widder, 2014; Vick-Majors et al., 2014; Banerjee et al., 2016a, b). An alternate approach has been applied by Allesina and Pascual (2009) who found that PageRank scores were strongly associated with extinction cascades in food webs. This provides a strong biological rationale for use of this metric. Additionally, we found that PageRank was strongly correlated with the other metrics of betweenness and centrality (Supplementary Figure S4). Therefore, we used PageRank scores as an indication of an organism’s importance to its network and defined keystone organisms as those in the 99th percentile of importance (PageRank score) for each community (network).
Environmental filtering appeared to influence keystone organisms more strongly than the overall microbial community. Habitat explained an increased portion of the variance in microbial community across all compartments when only keystone organisms were considered (Table 2). Similarly, clustering communities based on keystoneness (PageRank score) changed the grouping of rhizosphere communities such that sedges no longer clustered together, indicating that habitat was a more critical factor in determining keystone taxa in this case (Figure 7). On the other hand, in the phyllosphere, sphagnum-associated communities clustered farther apart from others when only keystone organisms were considered, suggesting that these organisms were more strongly influenced by their plant associate. It is not surprising that some organisms will be more strongly influenced by habitat and others by plant-associate, but these findings suggest that it is important to think carefully about which metrics we use to measure changes in community composition since they may alter our results.
We found some evidence that keystone organisms can be important to C-cycling as well as community assembly processes. Five keystone OTUs belonged to functional groups of interest, one methanogen and four putative methanotrophs. The methanogen was a keystone in bulk peat and all but one of the methanotrophs were keystones in rhizosphere and/or phyllosphere communities (Supplementary Table S1). Since the effect of methanotrophs can be large enough to remove up to 50–97% of CH4 from an ecosystem (Segarra et al., 2015), these keystone organisms may indeed play an important role in ecosystem processes.
Potential Impacts on CH4-Cycling
An important change taking place as permafrost thaws is conversion to increasingly hydric habitats with increasing methane (CH4) production (Johansson et al., 2006). This change has the potential to have a major impact on climate change since CH4 has a warming potential roughly 25 times that of CO2 (IPCC, 2013). In order to understand the implications of biotic and abiotic filtering of microbial communities on C-cycling in these systems we separately investigated organisms associated with CH4 cycling. Use of 16S data to infer functionality carries with it caveats since any inferences are based solely on taxonomic associations with organisms known to perform particular functions. However, complex traits such as methanogenesis and methanotrophy tend to be more highly conserved than other microbial traits, with methanogenesis in particular conserved across entire families or even orders (Martiny et al., 2015). Methanotrophy is also generally conserved within families, though there are examples where it has been lost by particular species (Tamas et al., 2014). While activity of an organism or any particular gene cannot be inferred from presence alone, previous research at this site has used metagenomic and metatranscriptomic data to show the presence and activity of the families identified here as putative methane-cycling organisms (Mondav et al., 2014; Singleton et al., 2018; Woodcroft et al., 2018b).
Putative methanotrophs identified in this dataset were members of groups previously identified as encoding genes for methanotrophy at this site (Singleton et al., 2018): Alpha-proteobacterial families Methylocystaceae, Beijerinckiaceae, Hyphomicrobiaceae, and Methylocystaceae; Gammaproteobacterial families Methylococcaceae and Crenotrichaceae; and the Verrucomicrobia order Methylacidi philales (Smith and Wrighton, 2019). We found putative methanotrophs across all habitats, plants, and compartments. This is perhaps unsurprising in an environment with high CH4 production. The most ubiquitous were members of Methylocystaceae, and this family also reached the highest relative abundances, which were in the phyllosphere – particularly when associated with A. polifolia and sedges (Eriophorum spp. and C. rotundata). In the palsa and fen, Methylacidiphilales and Crenotrichaceae, respectively, were also very common. Variation in methanotroph communities was driven by both habitat and plant associate, with the latter more important in the phyllosphere and the former more important in the rhizosphere and peat (Table 3). This suggests that methanotrophs in the rhizosphere are more influenced by soil parameters such as moisture and oxygenation, whereas methanotrophs in the phyllosphere are shaped by plant characteristics such as structure or biochemistry. If so, future investigations of functional gene abundance and activity may identify associations between plant functional types and microbial functional types. For instance, sedges have aerenchymous roots that may vent methane to the phyllosphere (Iversen et al., 2015) thereby shaping the microbial community.
Putative methanogens at this site were all members of genuses previously identified as methanogens at this site based on functional genes (Woodcroft et al., 2018a) belonging to Methanobacteriaceae (Methanobacterium), Methanoregulaceae (Candidatus Methanoregula), and Methanosarcinaceae (Methanosarcina). These organisms were distributed across all habitats and compartments, though they were rare in the phyllosphere and most abundant in the bulk soil (Figure 5b). Although methanogenesis is an anaerobic process, methanogens and even methane production have been previously observed in bulk aerobic environments such as the phyllosphere (Keppler et al., 2006; Taffner et al., 2019), and it is possible they inhabited anaerobic pockets of leaf tissue. Plant associate explained much more variation in relative abundances of putative methanogens in the rhizosphere than it did for any other group of organisms (Table 4). Tighter interaction between plant roots and putative methanogens than other groups could be driven by differences in substrate availability due to root exudates or differences in electron acceptor availability as a result of gas transport by roots. This strong influence of plant associate on explaining putative methanogen abundance has important implications for climate change since CH4 production from thawing permafrost ecosystems is a major potential feedback to global warming.
Conclusion
We began by asking to what degree plant community composition versus environmental shifts are important to structuring microbial communities under climate change. We additionally asked whether the relative contribution of each differs when considering only keystone organisms or particular functional groups versus the whole microbiome. Such questions are important because microbes are responsible for driving many important processes, such as C-cycling, which are changing rapidly with climate. If we can better understand which forces structure microbial communities or key functional groups, we can better predict the implications of changing conditions for these processes. Our results showed that such prediction is difficult, since much of the variation in microbial community composition remains unexplained by either biotic (plant) or environmental filtering. This may be an argument for investigating the role of more specific environmental variables such as moisture or pH in driving microbial community composition, but it should be noted that when investigating more specific drivers it becomes more likely that their effects will vary amongst microbial groups. We did show that both types of filtering play a role, with their relative importance differing between rhizosphere and phyllosphere and for particular functional groups but not for keystone organisms. Putative methane-cycling microbes were especially strongly impacted by plant associates beyond the effects of ecosystem transition, suggesting that plant community composition may be particularly important to shifts in the CH4 cycle. Further research should confirm the functionality and quantify activity levels of these putative CH4-cycling organisms when associated with different plant species. Such information may prove particularly useful to remote-sensing approaches, which could potentially use maps of plant community composition to infer information about CH4-cycling microbes. However, such approaches should consider microbial functional groups individually, rather than entire microbiomes, since the importance of plant versus environment seems to vary by group.
Data Availability Statement
The datasets generated for this study can be found in the IsoGenie Database: https://isogenie-db.asc.ohio-state.edu/, National Center for Biotechnology Information Sequence Read Archive, and Study Accession: PRJNA599435.
Author Contributions
AM, MH, VK-C, and VR designed the study in consultation with ED. AM and MH conducted fieldwork. AM and VK-C performed the lab work. AM, VK-C, and MH performed data processing, analyses, and generated figures. BB performed network analyses, including generation of associated figures and text. AM, MH, VK-C, VR, and SS secured funding for this work. All authors participated in manuscript preparation.
Funding
This work was funded by: a Wellesley College Daniel’s Fellowship to AM, a National Science Foundation Graduate Research Fellowship to MH, Wellesley College startup funds for VK-C, and the United States Department of Energy Office of Biological and Environmental Research Genomic Science Program grants DE-SC0010580 and DE-SC0016440.
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Acknowledgments
We would like to thank Abisko Scientific Research Station for critical support in fieldwork. Thanks also to the IsoGenie Investigators, Saleska Lab, Dr. Jana U’Ren, and Yueh-Fen Li for instrumental feedback on the analyses and writing.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmicb.2020.00796/full#supplementary-material
References
Aguilar, T. C., Rillig, M. C., and Ballhausen, M. B. (2017). Environmental filtering is a relic. A response to cadotte and tucker. Trends Ecol. Evol. 32, 882–884.
Allesina, S., and Pascual, M. (2009). Googling food webs: can an eigenvector measure species’ importance for Coextinctions? PLoS Comput. Biol. 5:e1000494. doi: 10.1371/journal.pcbi.1000494
Bäckstrand, K., Crill, P. M., Jackowicz-Korczyñski, M., Mastepanov, M., Christensen, T. R., and Bastviken, D. (2010). Annual carbon gas budget for a subarctic Peatland, northern Sweden. Biogeosciences 7, 95–108.
Baldrian, P. (2017). Forest microbiome: diversity, complexity and dynamics. FEMS Microbiol. Rev. 41, 109–130. doi: 10.1093/femsre/fuw040
Banerjee, S., Kirkby, C. A., Schmutter, D., Bissett, A., Kirkegaard, J. A., and Richardson, A. E. (2016a). Network analysis reveals functional redundancy and keystone taxa amongst bacterial and fungal communities during organic matter decomposition in an arable soil. Soil Biol. Biochem. 97, 188–198.
Banerjee, S., Mark, B., Carlyle, C. N., Bissett, A., Richardson, A. E., Siddique, T., et al. (2016b). Determinants of bacterial communities in Canadian Agroforestry systems. Environ. Microbiol. 18, 1805–1816. doi: 10.1111/1462-2920.12986
Banerjee, S., Schlaeppi, K., Marcel, G., and van der Heijden, A. (2018). Keystone taxa as drivers of microbiome structure and functioning. Nat. Rev. Microbiol. 16, 567–576. doi: 10.1038/s41579-018-0024-1
Berg, G., and Smalla, K. (2009). Plant species and soil type cooperatively shape the structure and function of microbial communities in the rhizosphere. FEMS Microbiol. Ecol. 68, 1–13. doi: 10.1111/j.1574-6941.2009.00654.x
Berry, D., and Widder, S. (2014). Deciphering microbial interactions and detecting keystone species with co-occurrence networks. Front. Microbiol. 5:219. doi: 10.3389/fmicb.2014.00219
Blume-Werry, G., Milbau, A., Teuber, L. M., Johansson, M., and Dorrepaal, E. (2019). Dwelling in the deep - strongly increased root growth and rooting depth enhance plant interactions with thawing permafrost soil. New Phytol. 223, 1328–1339. doi: 10.1111/nph.15903
Bragina, A., Berg, C., and Berg, G. (2015). The core microbiome bonds the alpine bog vegetation to a transkingdom metacommunity. Mol. Ecol. 24, 4795–4807. doi: 10.1111/mec.13342
Bragina, A., Berg, C., Cardinale, M., Shcherbakov, A., Chebotar, V., and Berg, G. (2012). Sphagnum mosses harbour highly specific bacterial diversity during their whole lifecycle. ISME J. 6, 802–813. doi: 10.1038/ismej.2011.151
Bragina, A., Oberauner-Wappis, L., Zachow, C., Halwachs, B., Thallinger, G., Müller, H., et al. (2014). The sphagnum microbiome supports bog ecosystem functioning under extreme conditions. Mol. Ecol. 23, 4498–4510. doi: 10.1111/mec.12885
Browne, F., Wang, H., Zheng, H., Roehe, R., Dewhurst, R. J., and Walsh, P. (2016). “A network analysis of methane and feed conversion genes in the rumen microbial community,” in Proceedings of the 2016 IEEE International Conference on Bioinformatics and Biomedicine (BIBM), Shenzhen, 1477–1484.
Bukin, Y. S., Galachyants, I. V., Morozov, S. V., Bukin, A. S., Zakharenko, P., and Zemskaya, T. I. (2019). The effect of 16S rRNA region choice on bacterial community metabarcoding results. Sci. Data 6:190007. doi: 10.1038/sdata.2019.7
Cadotte, M. W., and Tucker, C. M. (2017). Should environmental filtering be abandoned? Trends Ecol. Evol. 32, 429–437. doi: 10.1016/j.tree.2017.03.004
Caporaso, J. G., Kuczynski, J., Stombaugh, J., Bittinger, K., Bushman, F. D., Costello, E. K., et al. (2010). QIIME allows analysis of high-throughput community sequencing data. Nat. Methods 7, 335–336.
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
Cheng, W., Parton, W. J., Gonzalez-Meler, M. A., Phillips, R., Asao, S., Gordon, G., et al. (2014). Synthesis and modeling perspectives of rhizosphere priming. New Phytol. 201, 31–44. doi: 10.1111/nph.12440
Csardi, G., and Nepusz, T. (2006). The igraph software package for complex network research. InterJournal 1695, 1–9.
Delmotte, N., Knief, C., Chaffron, S., Innerebner, G., Roschitzki, B., Schlapbach, R., et al. (2009). Community proteogenomics reveals insights into the physiology of phyllosphere bacteria. Proc. Natl. Acad. Sci. U.S.A. 106, 16428–16433. doi: 10.1073/pnas.0905240106
DeSantis, T. Z., Hugenholtz, N., Larsen, M., Rojas, E. L., Brodie, K., Keller, T., et al. (2006). Greengenes, a Chimera-checked 16S rRNA gene database and workbench compatible with ARB. Appl. Environ. Microbiol. 72, 5069–5072.
Estrada-Pena, A., Cabezas-Cruz, A., Pollet, T., Vayssier-Taussat, M., and Cosson, J.-F. (2018). High throughput sequencing and network analysis disentangle the microbial communities of ticks and hosts within and between ecosystems. Front. Cell. Infect. Microbiol. 8:236. doi: 10.3389/fcimb.2018.00236
Evans, P. N., Boyd, J. A., Leu, A., Woodcroft, B. J., Donovan, H., Hugenholtz, P., et al. (2019). An evolving view of methane metabolism in the Archaea. Nat. Rev. Microbiol. 17, 219–232. doi: 10.1038/s41579-018-0136-7
Ferrera-Rodríguez, O., Greer, C. W., Juck, L. D., Consaul, E., Martínez-Romero, L., and Whyte, L. G. (2013). Hydrocarbon-degrading potential of microbial communities from Arctic plants. J. Appl. Microbiol. 114, 71–83. doi: 10.1111/jam.12020
Fischer, M. A., Güllert, S., Sven, C., Neulinger, S., Wolfgang, R., Streit, T., et al. (2016). Evaluation of 16S rRNA gene primer pairs for monitoring microbial community structures showed high reproducibility within and low comparability between datasets generated with multiple Archaeal and bacterial primer Pairs. Front. Microbiol. 7:1297. doi: 10.3389/fmicb.2016.01297
Fonseca, G., Citlali, D., Coleman, D., Etzel, G., Visel, A., Susannah, G., et al. (2016). The cacti microbiome: interplay between habitat-filtering and host-specificity. Front. Microbiol. 7:150. doi: 10.3389/fmicb.2016.00150
Friedman, J., and Alm, E. J. (2012). Inferring correlation networks from genomic survey data. PLoS Comput. Biol. 8:e1002687. doi: 10.1371/journal.pcbi.1002687
Holland-Moritz, H., Stuart, J., Lewis, L. R., Miller, S., Mack, M. C., Mcdaniel, S. F., et al. (2018). Novel bacterial lineages associated with boreal moss species. Environ. Microbiol. 20, 2625–2638.
Hugelius, G., Strauss, S., Zubrzycki, J. W., Harden, E. A. G., Schuur, C. L., Ping, L., et al. (2014). Estimated stocks of circumpolar permafrost carbon with quantified uncertainty ranges and identified data gaps. Biogeosciences 11, 6573–6593.
Hultman, J., Waldrop, M. P., Mackelprang, R., David, M. M., McFarland, J., Blazewicz, S. J., et al. (2015). Multi-omics of permafrost, active layer and thermokarst bog soil microbiomes. Nature 521, 208–212. doi: 10.1038/nature14238
IPCC (2013). “Climate change 2013: the physical science basis,” in Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change, eds T. F. Stocker, D. Qin, G.-K. Plattner, M. Tignor, S. K. Allen, J. Boschung, et al. (Cambridge: Cambridge University Press).
Iversen, C. M., Sloan, V. L., Sullivan, P. F., Euskirchen, E. S., Mcguire, A. D., Norby, R. J., et al. (2015). The unseen iceberg: plant roots in Arctic Tundra. New Phytol. 205, 34–58. doi: 10.1111/nph.13003
Oksanen, J., Guillaume Blanchet, F., Friendly, M., Kindt, R., Legendre, P., McGlinn, D., et al. (2019). vegan: Community Ecology Package, R package version 2.5-5. Available online at: https://CRAN.R-project.org/package=vegan
Johansson, T., Malmer, N., Crill, P. M., Friborg, T., Jonas, H., Mikhail, M., et al. (2006). Decadal vegetation changes in a northern Peatland, greenhouse gas fluxes and net radiative forcing. Glob. Chang. Biol. 12, 2352–2369.
Keppler, F., Hamilton, J. T. G., Brass, M., and Röckmann, T. (2006). Methane emissions from terrestrial plants under aerobic conditions. Nature 439, 187–191.
Knief, C. (2015). Diversity and habitat preferences of cultivated and uncultivated aerobic methanotrophic bacteria evaluated based on pmoA as molecular marker. Front. Microbiol. 6:1346. doi: 10.3389/fmicb.2015.01346
Knief, C., Delmotte, N., Chaffron, S., Stark, M., Innerebner, G., Wassmann, R., et al. (2012). Metaproteogenomic analysis of microbial communities in the phyllosphere and rhizosphere of rice. ISME J. 6, 1378–1390. doi: 10.1038/ismej.2011.192
Kögel-Knabner, I. (2002). The macromolecular organic composition of plant and microbial residues as inputs to soil organic matter. Soil Biol. Biochem. 34, 139–162.
Kumar, M., Brader, G., Sessitsch, A., Mäki, A., Van Elsas, J., and Nissinen, R. (2017). Plants assemble species specific bacterial communities from common core taxa in three Arcto-alpine climate zones. Front. Microbiol. 8:12. doi: 10.3389/fmicb.2017.00012
Kumar, M., Männistö, M. K., Van Elsas, J., and Nissinen, R. (2016). Plants impact structure and function of bacterial communities in arctic soils. Plant Soil 399, 319–332.
Lebeis, S. L. (2015). Greater than the sum of their parts: characterizing plant microbiomes at the community-level. Curr. Opin. Plant Biol. 24, 82–86. doi: 10.1016/j.pbi.2015.02.004
Lee, H. J., Jeong, S. E., Kim, P., Madsen, E. L., and Jeon, C. (2015). High resolution depth distribution of bacteria, archaea, methanotrophs, and methanogens in the bulk and rhizosphere soils of a flooded rice paddy. Front. Microbiol. 6:639. doi: 10.3389/fmicb.2015.00639
Lozupone, C., and Knight, R. (2005). UniFrac: a new phylogenetic method for comparing microbial communities. Appl. Environ. Microbiol. 71, 8228–8235.
Mackelprang, R., Waldrop, M. P., DeAngelis, K. M., David, M. M., Chavarria, K. L., Blazewicz, S. J., et al. (2011). Metagenomic analysis of a permafrost microbial community reveals a rapid response to thaw. Nature 480, 368–371. doi: 10.1038/nature10576
Martín González, A. M., Dalsgaard, B., and Olesen, J. M. (2010). Centrality measures and the importance of generalist species in pollination networks. Ecol. Complexity 7, 36–43.
Martiny, A., Jennifer, B. H., Stuart, E., Jones, J. T., Lennon, M., and Adam, C. (2015). Microbiomes in light of traits : a phylogenetic perspective. Science 350, 649–657. doi: 10.1126/science.aac9323
McCalley, C. K., Woodcroft, B. J., Hodgkins, S. B., Wehr, R., Kim, E., Mondav, R., et al. (2014). Methane dynamics regulated by microbial community response to permafrost thaw. Nature 514, 478–481. doi: 10.1038/nature13798
McDonald, D., Price, M. N., Goodrich, J., Nawrocki, E. P., DeSantis, T. Z., Probst, A., et al. (2012). An improved greengenes taxonomy with explicit ranks for ecological and evolutionary analyses of bacteria and archaea. ISME J. 6, 610–618. doi: 10.1038/ismej.2011.139
McKinney, W. G. (2010). “Data structures for statistical computing in python,” in Proceedings of the 9th Python in Science Conference, Vol. 445, Austin, TX, 51–56.
Mondav, R., McCalley, C. K., Hodgkins, S. B., Frolking, S., Scott, R., Saleska, V., et al. (2017). Microbial network, phylogenetic diversity and community membership in the active layer across a permafrost Thaw Gradient. Environ. Microbiol. 19, 3201–3218. doi: 10.1111/1462-2920.13809
Mondav, R., Woodcroft, B. J., Kim, E., McCalley, C. K., Hodgkins, S. B., Patrick, M., et al. (2014). Discovery of a novel methanogen prevalent in thawing permafrost. Nat. Commun. 5:3212. doi: 10.1038/ncomms4212
Monteux, S., Weedon, J. T., Blume-Werry, G., Gavazov, K., Vincent, E. J., Johansson, M., et al. (2018). Long-term in situ permafrost thaw effects on bacterial communities and potential aerobic respiration. ISME J. 13, 2140–2142. doi: 10.1038/s41396-018-0176-z
Nissinen, R. M., Männistö, M. K., and Van Elsas, J. (2012). Endophytic bacterial communities in three arctic plants from low arctic fell tundra are cold-adapted and host-plant specific. FEMS Microbiol. Ecol. 82, 510–522. doi: 10.1111/j.1574-6941.2012.01464.x
Nuccio, E. E., Anderson-Furgeson, J., Estera, K. Y., Pett-Ridge, J., De Valpine, P., Eoin, L., et al. (2016). Climate and edaphic controllers influence rhizosphere community assembly for a wild annual grass. Ecology 97, 1307–1318.
Parada, A. E., Needham, D. M., and Fuhrman, J. (2016). Every base matters: assessing small subunit rRNA primers for marine microbiomes with mock communities, time series and global field samples. Environ. Microbiol. 18, 1403–1414. doi: 10.1111/1462-2920.13023
Parniske, M. (2008). Arbuscular Mycorrhiza: the mother of plant root endosymbioses. Nat. Rev. Microbiol. 6, 763–775. doi: 10.1038/nrmicro1987
Pegoraro, E., Mauritz, M., Bracho, R., Ebert, C., Dijkstra, P., Hungate, B. A., et al. (2019). Glucose addition increases the magnitude and decreases the age of soil respired carbon in a long-term permafrost incubation study. Soil Biol. Biochem. 129, 201–211.
Philippot, L., Raaijmakers, J. M., Lemanceau, P., and van der Putten, W. H. (2013). Going back to the Roots: the microbial ecology of the rhizosphere. Nat. Rev. Microbiol. 11, 789–799. doi: 10.1038/nrmicro3109
R Core Team (2019). R: A Language and Environment for Statistical Computing. Vienna: R Foundation for Statistical Computing.
Redford, A. J., Bowers, R. M., Knight, R., Linhart, Y., and Fierer, N. (2010). The ecology of the phyllosphere: geographic and phylogenetic variability in the distribution of bacteria on tree leaves. Environ. Microbiol. 12, 2885–2893. doi: 10.1111/j.1462-2920.2010.02258.x
Schuur, E. A. G., McGuire, C., Schadel, G., Grosse, J. W., Harden, D. J., Hayes, G., et al. (2015). Climate change and the permafrost carbon feedback. Nature 520, 171–179. doi: 10.1038/nature14338
Segarra, K. E., Schubotz, F., Samarkin, V., Yoshinaga, M. Y., Hinrichs, K.-U., and Joye, S. B. (2015). High rates of anaerobic methane oxidation in freshwater wetlands reduce potential atmospheric methane emissions. Nat. Commun. 6:7477. doi: 10.1038/ncomms8477
Segata, N., Izard, J., Waldron, L., Gevers, D., Miropolsky, L., Garrett, W. S., et al. (2011). Metagenomic biomarker discovery and explanation. Genome Biol. 12:R60. doi: 10.1186/gb-2011-12-6-r60
Shade, A., Jones, S. E., Caporaso, J. G., Handelsman, J., Knight, R., Fierer, N., et al. (2014). Conditionally rare taxa disproportionately contribute to temporal changes in microbial diversity. mBio 5:e01371-14. doi: 10.1128/mBio.01371-14
Sinclair, L., Osman, O. A., Bertilsson, S., and Eiler, A. (2015). Microbial community composition and diversity via 16S rRNA gene amplicons: evaluating the illumina platform. PLoS One 10:e0116955. doi: 10.1371/journal.pone.0116955
Singleton, C. M., McCalley, C. K., Woodcroft, B. J., Boyd, J. A., Evans, P. N., Hodgkins, S. B., et al. (2018). Methanotrophy across a natural permafrost thaw environment. ISME J. 12, 2544–2558. doi: 10.1038/s41396-018-0065-5
Smith, G. J., and Wrighton, K. C. (2019). Metagenomic approaches unearth methanotroph phylogenetic and metabolic diversity. Curr. Issues Mol. Biol. 33, 57–84. doi: 10.21775/cimb.033.057
Strickland, M. S., Lauber, C., Fierer, N., and Bradford, M. A. (2009). Testing the functional significance of microbial community composition. Ecology 90, 441–451.
Taffner, J., Cernava, T., Erlacher, A., and Berg, G. (2019). Novel insights into plant-associated archaea and their functioning in arugula (Eruca Sativa Mill.). J. Advert. Res. 19, 39–48. doi: 10.1016/j.jare.2019.04.008
Tamas, I., Smirnova, A. V., He, Z., and Dunfield, P. F. (2014). The (d)evolution of methanotrophy in the Beijerinckiaceae–a comparative genomics analysis. ISME J. 8, 369–382. doi: 10.1038/ismej.2013.145
Tarnocai, C., Canadell, J. G., Schuur, E. G., Kuhry, P., Mazhitova, G., Zimov, S., et al. (2009). Soil organic carbon pools in the northern circumpolar permafrost region. Glob. Biogeochem. Cycles 23, 1–11.
Thakur, M. P., and Wright, A. J. (2017). Environmental filtering, niche construction, and trait variability: the missing discussion. Trends Ecol. Evol. 32, 884–886.
Thiessen, S., Gleixner, G., Wutzler, T., and Reichstein, M. (2013). Both priming and temperature sensitivity of soil organic matter decomposition depend on microbial biomass E an incubation study. Soil Biol. Biochem. 57, 739–748.
Tkacz, A., Cheema, J., Chandra, G., Grant, A., and Poole, P. S. (2015). Stability and succession of the rhizosphere microbiota depends upon plant type and soil composition. ISME J. 9, 2349–2359. doi: 10.1038/ismej.2015.41
Tláskal, V., Vořís̃ková, J., and Baldrian, P. (2016). Bacterial succession on decomposing leaf litter exhibits a specific occurrence pattern of cellulolytic taxa and potential decomposers of fungal mycelia. FEMS Microbiol. Ecol. 92, 1–10. doi: 10.1093/femsec/fiw177
Tveit, A., Schwacke, R., Svenning, M. M., and Urich, T. (2013). Organic carbon transformations in High-Arctic peat soils: key functions and microorganisms. ISME J. 7, 299–311. doi: 10.1038/ismej.2012.99
Vick-Majors, T. J., Priscu, J. C., and Amaral-Zettler, L. A. (2014). Modular community structure suggests metabolic plasticity during the transition to polar night in ice-covered Antarctic lakes. ISME J. 8, 778–789. doi: 10.1038/ismej.2013.190
Virtanen, P., Gommers, R., Oliphant, T. E., Haberland, M., Reddy, T., Cournapeau, D., et al. (2020). SciPy 1.0: fundamental algorithms for scientific computing in python. Nat. Methods 17, 261–272. doi: 10.1038/s41592-019-0686-2
Walt, S., Colbert, S. C., and Varoquaux, G. (2011). The NumPy array: a structure for efficient numerical computation. Comput. Sci. Eng. 13, 22–30.
Wang, H., Zheng, H., Browne, F., Roehe, R., Dewhurst, R. J., Engel, F., et al. (2017). Integrated metagenomic analysis of the rumen microbiome of cattle reveals key biological mechanisms associated with methane traits. Methods 124, 108–119. doi: 10.1016/j.ymeth.2017.05.029
Woodcroft, B. J., Singleton, C. M., Boyd, J. A., Evans, P. N., Hoelzle, R. D., Timothy, O., et al. (2018a). Genome-centric metagenomic insights into microbial carbon processing across a permafrost thaw gradient. Nature 560, 49–54.
Woodcroft, B. J., Singleton, C. M., Boyd, J. A., Evans, P. N., Emerson, J. B., Zayed, A. A. F., et al. (2018b). Genome-centric view of carbon processing in thawing permafrost. Nature 560, 49–54. doi: 10.1038/s41586-018-0338-1
Keywords: microbial community assembly, permafrost thaw, plant–microbial interactions, keystone species, climate feedbacks, IsoGenie Project, Stordalen Mire
Citation: Hough M, McClure A, Bolduc B, Dorrepaal E, Saleska S, Klepac-Ceraj V and Rich V (2020) Biotic and Environmental Drivers of Plant Microbiomes Across a Permafrost Thaw Gradient. Front. Microbiol. 11:796. doi: 10.3389/fmicb.2020.00796
Received: 15 November 2019; Accepted: 03 April 2020;
Published: 15 May 2020.
Edited by:
Giorgio Gambino, Institute for Sustainable Plant Protection (IPSP-CNR), ItalyReviewed by:
Andrew Ogram, University of Florida, United StatesChristopher Ryan Penton, Arizona State University, United States
Copyright © 2020 Hough, McClure, Bolduc, Dorrepaal, Saleska, Klepac-Ceraj and Rich. 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: Moira Hough, houghm@email.arizona.edu; Vanja Klepac-Ceraj, vklepacc@wellesley.edu; Virginia Rich, rich.270@osu.edu; virginia.isabel.rich@gmail.com