Dajiuhu (31°24′ – 31°33′ N, 109°56′ – 110°11′ E) is a typical subalpine peatland located in Shennongjia, Hubei province, central China (Supplementary Figure S1A). Peat has developed in a closed intermontane basin (1,730 m, above sea level), with an area around 16 km² and a depth of 2–3 m (Huang et al., 2018). The mean annual temperature and precipitation are 7.2°C and 1,560 mm, respectively. The present vegetation types are mainly Carex spp., S. palustre, Sanguisorba officinalis, and Euphorbia esula accompanying shrubs. More information about the Dajiuhu Peatland was given by Huang et al. (2012) and Li et al. (2016).

In this study, S. palustre and underlying peat samples (0–5 cm) were collected with sterile plastic bags and centrifuge tubes, respectively. Samples were collected in triplicate, in July 2016. Sampling sites were located at the peatland core area are Erhaoba (EHB), Niangniangfen (NNF) and two sites at Yangluchang (YLC) where S. palustre was actively growing (Supplementary Figure S1A). Water table at the four sites was determined in situ at the time of sampling. Visible plant root, litter and debris were removed from S. palustre peat samples. All the samples were transported to the geomicrobiology laboratory at China University of Geosciences (Wuhan) with ice box within 12 h.

To get the microbes from the ectosphere, S. palustre was rinsed three times for 5 min each time with sterile deionized water to remove adherent debris. Then approx. 5 g of S. palustre samples were placed in 20 mL sterile phosphate buffered saline (PBS), and shaken and sonicated to separate the tightly adhering microbes (Bulgarelli et al., 2015; Edwards et al., 2015). The S. palustre plant material was transferred into another tube, and PBS solution was centrifuged at 10,000 g for 5 min and the supernatant was discarded (Supplementary Figure S1B). The remaining PBS fraction was stored at −80°C until DNA extraction.

The retrieved S. palustre designated for endosphere underwent further sonication and shake twice following the above-depicted approach. The S. palustre samples were then placed into 75% ethanol and 20% hydrogen peroxide for 1 min, respectively (Bulgari et al., 2009). The S. palustre samples were rinsed three times for 5 min each time with sterile deionized water and the final sterile water was inoculated on the R2A medium. The S. palustre plant materials were stored at −80°C until DNA extraction.

Additional S. palustre peat samples (approx. 20 g) were dried at 50°C for 24 h and subject to the organic matter and total nitrogen analyses.

Total DNA for each approximately 0.5 g freeze-dried S. palustre peat, ectosphere and endosphere fractions were extracted by PowerSoil DNA Isolation Kit (QIAGEN, Düsseldorf, Germany). The concentration and quality of DNA were checked with NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, United States) and 1% agarose gel stained with ethidium bromide. The V3-V4 regions of bacterial 16S rRNA gene were amplified with the primer set of 347F (5′-CCTACGGRRBGCASCAGKVRVGAAT-3′) and 802R (5′-GGACTACNVGGGTWTCTAATCC-3′) (Ren et al., 2019). PCR reactions contained 25 μL 2× Premix Taq (Takara Biotechnology, Dalian Co., Ltd., China), 1 μL each primer (10 mM) and 3 μL DNA (20 ng/μL) template in a volume of 50 μL. PCR products were then purified and quantified with Quit 2.0 fluorometer (Invitrogen, Carlsbad, CA, United States). DNA libraries were generated using NEBNext Ultra DNA Library Prep Kit for Illumina (New England Biolabs, Ipswich, MA, United States) following manufacturer’s protocols and the quality of libraries was validated by the Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, United States) and the Qubit 2.0 Fluorometer. DNA libraries were sequenced on an Illumina MiSeq platform with PE300 paired-end reads at GENEWIZ, Inc. (Suzhou, China).

Raw sequences were processed in Quantitative Insight Into Microbial Ecology (QIIME v. 1.9.1) platform (Caporaso et al., 2010). The sequences were assigned to each sample based on their unique barcode, followed by the removal of primers and barcodes. The forward and reverse sequences were paired-end assembly with overlaps more than 20 bp. The high quality reads were determined according to without any ambiguous base “N”, length >200 bp and base quality score >20. The clean sequences were clustered into operational taxonomic units (OTUs) by UCLUST (Edgar, 2010) based on 3% sequence dissimilarity and singletons (OTU with only one sequence in all samples) were eliminated. The mitochondria and chloroplast were removed and the chimera sequences were identified by UCHIME algorithm (Edgar et al., 2011). Represent sequences for each OTU were classified using Ribosomal Database Project classifier at a threshold of 0.8 (Cole et al., 2009) and annotated in the SILVA 119 database. All the samples were resampled to 29,577 sequences per sample by mothur (Schloss et al., 2009).

Raw 16S rRNA gene data was available at the NCBI Sequence Read Archive database¹ with the accession number of PRJNA512496.

We set relative abundance thresholds as 0.01% for rare taxa (RT) and 1% for abundant taxa (AT) and classified all OTUs (archaea and unclassified OTUs excluded) into six subcategories (Dai et al., 2016; Xue et al., 2018) to compare the roles of taxa with different abundance in the whole bacterial communities (WBC): always abundant taxa (AAT) with a relative abundance >1% in all samples; conditionally abundant taxa (CAT) with a relative abundance ≥1% in some samples and always ≥0.01% in all samples; always rare taxa (ART) with a relative abundance <0.01% in all samples; conditionally rare taxa (CRT) with a relative abundance <0.01% in some samples but never abundant (≥1%) in any samples; moderate taxa (MT) with a relative abundance between 0.01% and 1% in all samples; conditionally abundant and rare taxa (CART) with a relative abundance varying from rare (<0.01%) to abundant (≥1%). The detailed depictions of OTU category were listed in Supplementary Table S1.

The contribution of each subcategory to community variations was determined by similarity percentage (SIMPER) analysis in PAST 3.26². The alpha diversity indices (Observed species, Chao1, ACE, Simpson, Pielou, and Shannon), NMDS, PERMANOVA, and ANOSIM analyses were performed using “vegan” package in R 3.5.1³. Functional profiles were predicted using Tax4Fun package (Aßhauer et al., 2015) and principal component analysis (PCA) was employed to reveal compartment differences of potential functional communities based on Euclidean distance of the relative abundance of overall Kyoto Encyclopedia of Genes and Genomes (KEGG) Orthologs (KOs) with “ggord” package in R 3.5.1. Differences in alpha diversity and potential function were tested with one-way analysis of variance (ANOVA) in SPSS 18⁴. Beta diversity analysis was performed based on Bray-Curtis distance of the relative abundance of OTUs in non-metric multidimensional scaling (NMDS). Analysis of similarity (ANOSIM) was used to test the differences among groups. Monte–Carlo tests of redundancy analysis determined the relationships between environmental factors and bacterial communities of S. palustre peat in Canoco 5⁵.

The OTUs with relative abundance above 0.01% were selected, and the correlation matrix between two OTUs was calculated in “psych” R package. Correlation coefficient (r > 0.8 or r < −0.8) and significance (FDR-corrected P < 0.01) were integrated into network analysis (Jiao et al., 2016). Topology attributes of the network (numbers of node and edge, average degree, network diameter, graph density, modularity, clustering coefficient and average of path length) were calculated using “igraph” package in R 3.5.1. Network visualization was conducted with Gephi 0.9.2⁶. Meanwhile, 1,000 Erdös-Réyni random networks were constructed to compare with the real network (Erdös and Réyni, 1960). Nodes with high degree (>100) and low betweenness centrality (<5,000) were considered as the keystone species in co-occurrence networks (Ma et al., 2016).

In total, 2,413,649 high quality sequences and 3,144 operational taxonomic units (OTUs) were obtained. After the removal of archaea and unclassified OTUs, 2,848 OTUs were subject to further analysis. Among these, 13 OTUs representing 9.01% of all sequences were abundant taxa (AT) and present in all samples, whereas 2,706 OTUs occupying 39.2% of all sequences belonged to rare taxa (RT). Conditionally abundant and rare taxa (CART) included 126 OTUs representing 51.13% of all sequences (Supplementary Table S1).

Significant differences were observed in alpha diversities of whole bacterial communities (WBC) and RT among compartments (Supplementary Table S2). Observed species, Simpson, Pielou, and Shannon indices of WBC and CART significantly decreased from Sphagnum palustre peat (SP) to S. palustre endosphere (Endo), with those of S. palustre ectosphere (Ecto) in-between (Figure 1). In contrast, there were no significant differences in alpha diversity across sampling sites (Supplementary Table S2). Venn diagrams indicated over 50% OTUs in WBC were shared in three compartments, and the distributions of shared and specific OTUs in RT were similar to those of WBC (Figure 2B).

The dissimilarities in bacterial communities across compartments were visualized by NMDS analysis based on Bray-Curtis matrix. The communities of the same compartment clustered together (Figure 2A). The pairwise comparison of between group distances (Figure 2C), and the results of ANOSIM (Table 1) indicated that bacterial communities were significantly different among compartments, rather than among sampling sites.

As for WBC, Acidobacteria and Alphaproteobacteria dominated SP, whereas the predominant phyla in Ecto were Cyanobacteria, Alphaproteobacteria, and Gammaproteobacteria, followed by Acidobacteria. In Endo, Cyanobacteria and Alphaproteobacteria dominated the communities, followed by Acidobacteria (Figure 3A). Alphaproteobacteria was commonly found across all compartments and Acidobacteria presented frequently in SP, whereas Cyanobacteria was abundant in Ecto and Endo.

In total, 13 OTUs were classified into AT, which were affiliated to Actinobacteria, Alphaproteobacteria and Gammaproteobacteria phyla (Figure 3B), and accounted for 0.46% of the OTU richness (Supplementary Table S3). As many as 2,706 OTUs were identified as RT, which were affiliated with 30 phyla (Figure 3C) and represented 95.01% of the OTU richness (Supplementary Table S3). Furthermore, 126 OTUs were recognized as CART belonging to 13 phyla (Figure 3D), and represented 4.42% of the OTU richness (Supplementary Table S3).

A co-occurrence network was constructed based on correlation relationships among OTU relative abundances. The resulting network contained 674 nodes linked by 8,455 edges, where positive correlations were commonly observed compared to negative ones (Figure 4C). The scale-free degree distribution of the real network (power-law: R² = 0.945, P < 0.0001, Supplementary Figure S2) indicated a scale-free and non-random network structure. Meanwhile the network had small-world features and modular structure as indicated by higher values of the modularity index (0.762) and clustering coefficient (0.521) compared with those of identical size Erdös-Réyni random networks (Erdös and Réyni, 1960; Table 2).

The top 10 phyla in the network accounted for 93% of all nodes (Figure 4A). The whole network was grouped into 4 major modules (Figure 4B). Module I included 36.35% of the nodes in the entire network, which mainly fell into Proteobacteria. Module II and module III were mainly comprised of Proteobacteria and Acidobacteria occupying 41.99% of the nodes in the whole network, whereas module IV mainly consisted of Proteobacteria and Bacteroidetes (Supplementary Figure S3A). Ternary plots showed the specificity of different modules to a particular compartment (Supplementary Figures S3B–E). For example, most OTUs in module I showed high abundances in Ecto and Endo, whereas the majority of OTUs in modules II, III, and IV showed high proportions in SP.

Node-level topological features of microbial groups with different abundances (RT, CART, AT, and MT) indicated that the degree, betweenness and eigenvector centrality values of CART were significantly higher than that of RT (Supplementary Figures S4A,B,D). The closeness centrality values of AT were the highest among 4 microbial groups (Supplementary Figure S4C). Based on high degree (>100) and low betweenness centrality (<5,000) in the co-occurrence networks, 5 OTUs were identified as keystone species including Acidobacteria (1 OTU), Betaproteobacteria (1 OTU), Alphaproteobacteria (1 OTU), Ignavibacteriae (1 OTU) and Spirochaetae (1 OTU), which belonged to RT and CART (Supplementary Table S4).

Monte–Carlo tests indicated that three physicochemical parameters, i.e., water table, total nitrogen and organic matter contents were significantly correlated with the variations of WBC, RT, CART, and AT in SP samples (Table 3). The WBC was significantly affected by water table and total nitrogen content (P < 0.05), which were also identified as the significant determinants (P < 0.05) in RT and CART. In contrast, organic matter content rather than water table and total nitrogen significantly (P < 0.05) shaped the structure of AT.

Similar to the structural variations of communities, the potential functions of communities showed an obvious separation among WBC, RT, and CART (Supplementary Figures S5A,C,D) compared to that in AT across compartments (Supplementary Figure S5B) as indicated by principal component analysis (PCA) based on Euclidean distance. The dissimilarity comparisons in each two compartments and sampling sites were further quantified using PERMANOVA (Supplementary Table S5), which showed significant (α = 0.05) differences in functional potentials across compartments rather than those among sampling sites.

Among WBC, RT, and CART, the relative abundance of potential nitrogen-fixing functional groups was the highest in SP, followed by those in Endo, and the lowest was in Ecto (Figure 5A). The relative abundance of predicted methane-oxidizing functional groups significantly decreased from Endo to SP, with that of Ecto in-between (Figure 5B). Together, specific potential functional groups such as N₂-fixers and CH₄-oxidizers in RT and CART can reflect consistent variation patterns with those in WBC across compartments, whereas AT did not.

Our results validated that compartment effects drove the shift of bacterial communities in S. palustre as indicated by alpha and beta diversities. Alpha diversity decreased from SP to Ecto and was the lowest in Endo (Figure 1). This result agrees with the view of plant-associated microbial colonization, which shows a reduction in the species diversity from rhizosphere to rhizoplane to endosphere (Edwards et al., 2015; Beckers et al., 2017). The changes of alpha diversity from SP to Ecto to Endo may closely relate to the specific tissues of Sphagnum comprised of chlorocytes and hyalocytes. Microbial colonization in Sphagnum is strictly selected by chlorocytes (Opelt et al., 2007b; Bragina et al., 2012), and reduced the bacterial diversity in Sphagnum compartments (Figure 1). The significant increase of community dissimilarities (beta diversity) between compartments from Ecto vs. Endo to SP vs. Ecto and further to SP vs. Endo (Figure 2C and Table 1) also confirmed plant compartment effects on bacterial communities. Each compartment is a unique and heterogeneous microenvironment (Coleman-Derr et al., 2016; Stacie et al., 2016), which selects their specific microbial communities (Beckers et al., 2017). The plant ectosphere is usually oligotrophic with frequent fluctuations in temperature, precipitation, moisture and radiation, which is challenging for microbial survival (Müller et al., 2016). Therefore microorganisms living on the surface of plants have to deal with hard conditions, which resulted in a dispersive distribution of bacterial communities in Ecto (Figure 2A). In contrast, plant cells can act as oases for microbial survival in peatlands. Usually, microorganisms with specific traits, e.g., chemotaxis, motility, and quorum sensing are recruited by the host-dependent selection or passive transport, and thus lead to the formation of a set of specific plant-associated microbial communities (Compant et al., 2010; Kardol and Wardle, 2010; Bulgarelli et al., 2012; Philippot et al., 2013; Müller et al., 2016). Unlike vascular plants, rootless Sphagnum mosses intercept and hold water and nutrients efficiently through non-photosynthetic cells and adjacent photosynthetic cells (Weston et al., 2015; Kostka et al., 2016), which can facilitate the capture of their endosymbiotic microorganisms.

Our results for the first time indicated that microbial taxa with various abundances accounted for different variations of microbial communities in a subalpine peatland. The richness of RT was more than 60 times higher than that of AT (Figure 1), suggesting their important contributions to microbial diversity (Lynch and Neufeld, 2015; Zhan and Macisaac, 2015). Variation patterns of microbial taxa with different abundances in alpha and beta diversities among S. palustre compartments (Figures 1, 2) supported the hypothesis 1 that rare species mainly contributed to bacterial community variations compared to abundant taxa.

The results of co-occurrence network analysis demonstrated keystone species were crucial to the whole bacterial communities and might exert large effects on other community components. Co-occurrence network can provide profound and unique insights into microbial interactions and community assembly (Barberán et al., 2012). Our results indicated that more positive interactions were observed between RT and CART (Figure 4C) than negative ones, which suggested a possible cooperation among them (Ju et al., 2014). Keystone species are characterized by high degree and low betweenness centrality based on scale-free feature (Barabási, 2009; Berry and Widder, 2014), and recognized as the initial components and affect other community components in networks. In our study, all the keystone species identified were affiliated to RT and CART, which confirmed that rare species might play a critical role in community co-occurrence network. Among keystone species, Candidatus Nitrotoga are cold-adapted nitrite-oxidizing bacteria, which involved in nitrogen cycling (Alawi et al., 2007). Spirochaeta are anaerobic bacteria, and some members are capable of cellulose-derivative degradation (glucose and cellobiose) and methane production in peat soils (Schmidt et al., 2015; Juottonen et al., 2017). These keystone species may benefit Sphagnum development under extreme conditions and mediate the recycle of macroelements (C, N) (Bragina et al., 2014) and sustain ecosystem network (Layeghifard et al., 2017). Therefore, rare species may be important in supporting peatland ecosystem functioning.

To be noted, rare species widely spread in all modules, which may affect the functioning of host-associated environments (Jousset et al., 2017). For example, nitrogen-fixing and methane-oxidizing microbes play essential roles in nitrogen supply for Sphagnum (Larmola et al., 2014; Vile et al., 2014) and decomposition of organic matter (Chen and Murrell, 2010). In the predicted functional communities, the variations in relative abundance of genes encoding nitrogenase and methane monooxygenase in RT and CART were consistent with those in WBC among compartments (Figure 5). Therefore, we surmised that rare species were important in the co-occurrence network and specific functions of bacterial communities, and the hypothesis 2 was supported. Although network analysis can give us some insight of microbial interactions, it only shows a statistical correlation among microbial groups and does not directly prove microbial interactions. Combination of microscopy-based or co-culture experiments with co-occurrence analysis will provide more direct evidences of microbial interactions.

Environmental conditions strongly impact the structure of microbial communities (Bahram et al., 2018). The effects of water table and total nitrogen content on bacterial communities in SP samples (Table 3) matches well with previous studies in different types of peatlands (Mishra et al., 2014; Urbanová and Bárta, 2016; Zhong et al., 2017). Interestingly we further found that RT and CART also responded to water table and total nitrogen content. In contrast, AT was sensitive to organic matter content. This indicated that microbial taxa with various abundances responded differently to environmental factors. It has been demonstrated that the content and form of organic matter affect the bioavailability and the quality of organic carbon, and further affect the priority of microbial metabolism (Stres et al., 2008). Our results further demonstrated that organic matter content mainly affect AT in WBC. Therefore, the hypothesis 3 that rare species and abundant taxa responding differently to environmental factors was supported.

Microbial groups with different abundances showed distinct environmental responses, which may result from the environmental filtering or dispersal limitation (Chen et al., 2017; Wu et al., 2017; Xue et al., 2018). Our results (Table 3) indicated RT and CART were more sensitive to environmental filtering than AT did. Usually RT and CART showed weak competition capability and low growth rate (Campbell et al., 2011; Logares et al., 2015), whereas AT can compete for more substrates and well adapt to an ecosystem via active growth and proliferation in a wide niche (Pedrós-Alió, 2006; Jousset et al., 2017; Xue et al., 2018). This suggests that the composition of microbial communities is strongly affected by environmentally induced species sorting (Roberto et al., 2018). Therefore, taxa with various abundances in bacterial communities from SP were impacted by unequal environmental filtering.