Degradation shaped bacterial and archaeal communities with predictable taxa and their association patterns in Zoige wetland at Tibet plateau

Soil microbes provide important ecosystem services. Zoige Plateau wetland, the largest alpine peat wetland in the world, has suffered from serious degradation in the past 30 years. We studied the composition of the Zoige Plateau alpine wetland soil microbiota and relations among specific taxa using 16S rRNA amplicon sequencing combined with association network analysis. Compared to the pristine swamp soil, taxons DA101, Aeromicrobium, Bradyrhizobium, and Candidatus Nitrososphaera were enriched and several methanogenic Euryarchaeota were depleted in the moderately degraded meadow soil and highly degraded sandy soil. Soil total potassium contents in soils with different degradation levels were significantly different, being the highest in meadow soil and lowest in swamp soil. The association network analysis showed that total potassium positively correlated with specific bacterial and archaeal taxa. Jiangella, Anaerolinea, Desulfobulbus, Geobacter, Flavobacterium, Methanobacterium and Methanosaeta were identified as the keystone genera in the networks. Soil degradation affected soil properties, and caused changes in the bacterial and archaeal community composition and the association patterns of community members. The changes could serve as early warning signals of soil degradation in alpine wetlands.


Results
Soil properties and community composition. Soil properties were different in the three soils associated with the degraded Zoige wetland (Table 1). Soil total potassium content (TK: 6.79 to 23.13 g kg −1 ) was highest in the meadow soil and lowest in swamp soil (Table 1). Principal component analysis (PCA) showed clear separation of the three soils indicating that wetland degradation resulted in distinct microhabitats (ANOSIM: R > 0.9, P < 0.001) (Fig. 1A). Proteobacteria (24.6-31.6%), Acidobacteria (12.3-19.7%), and Soil codes* WC (%) pH SOC (g.kg −1 ) TN (g.kg −1 ) OP (g.kg −1 ) TK (g.kg −1 ) AN (mg.kg −1 ) AP (mg.kg −1 ) AK (mg.kg −1 ) SW11 74.9 ± 0.81a 7.59 ± 0.09b 279.  S1). The relative abundances of Deltaproteobacteria and Betaproteobacteria were higher in the swamp soil than in the other soils. The relative abundance of Alphaproteobacteria was highest in the swamp soil, and that of Gammaproteobacteria in the sandy soil. The relative abundances of ten families were higher than 0.5% of the total abundance (Fig. S2). Compared to the pristine swamp soil, the relative abundance of Nitrososphaeraceae, Sphingomonadaceae, Xanthomonadaceae, and Chitinophagaceae were significantly higher in the degraded soils, while those of Sinobacteraceae and Comamonadaceae were lower (P < 0.05). Moreover, the relative abundances of Myxococcaceae (P < 0.05) and Rhodocyclaceae were lowest in the meadow soil, while those of Syntrophobacteraceae and Rhodospirillaceae were lowest in the swamp soil (Fig. S2).

Differentially abundant taxa.
Enriched taxa (e-taxa) and depleted taxa (d-taxa), representing taxa with relative abundances twice higher or lower (P = 0.05) in the swamp soil than in the meadow and sandy soils, were identified by differential abundance analysis. Using taxa counts from swamp soil as a control and an FDR-adjusted P < 0.1, eleven and 89 differentially abundant taxa were detected at the phylum and genus levels, respectively ( Fig. 2). At the phylum level, Euryarchaeota, WWE1, OP8 and GOUTA4 were enriched in the swamp soil when compared with the meadow soil (Figs 2A and S3A). Compared with the sandy soil, WWE1, Euryarchaeota, OP8, GOUTA4 and Bacteroidetes were enriched, and TM7 and AD3 were depleted in the swamp soil (Figs 2B and S3B). At the genus level, additional 35 e-taxa and 16 d-taxa were detected in the swamp soil when compared with the meadow soil (Fig. 2D). The top five e-taxa were DA101, Crocinitomix, Aeromicrobium, Bradyrhizobium, and Candidatus Nitrososphaera (phylum Crenarchaeota), and the top five d-taxa were Syntrophus, Methanobacterium, Methanosarcina, Chloronema, and Methanosaeta (Dataset S1). When compared with the sandy soil, 31 e-taxa and 6 d-taxa were detected in the swamp soil (Fig. 2E), among which DA101, Aeromicrobium, Sporichthya, JG37-AG-70, and Zoogloea were the top five e-taxa and Chryseobacterium, Syntrophus, Chloronema, Oryzihumus, and Methanosarcina the top five d-taxa (Dataset S2). When comparing sandy soil with the meadow soil, only one d-taxa was detected at the genus level (Fig. 2C,F).
The Mantel test showed significant correlations between differentially abundant phyla and multiple edaphic factors. Soil TK correlated negatively with most of the differentially abundant phyla (P = 0.001) (Fig. S3C,D).
Community structure, variation, and determinants. Principal coordinate analysis (PCoA) with weighted and unweighted UniFrac distances clearly separated the swamp communities from the communities in degraded soils (Fig. 3A,B). Weighted UniFrac takes the abundance of taxa into account whereas unweighted UniFrac is based on the presence and absence of species, making it more sensitive to rare taxa (Fig. 3B). The clear separation between the swamp soil and the degraded soils in unweighted Unifrac-based PCoA indicates that wetland degradation affects especially the rare taxa. In the distance-based redundancy analysis (dbRDA), soil TK was the strongest driver of microbial community structure among the constrained nine edaphic factors (Fig. 3C,D). In addition, Mantel test showed that soil TK correlated positively and significantly with soil microbial community composition (Fig. S4A,B). These results further suggest that the total potassium content of the soils was the key factor in shaping the microbial community composition in the wetland soils.
Network associations. Based on differential taxa abundance analysis, association networks were generated to characterize bacterial and archaeal communities in the soils (Fig. 4). The topological features of the network are presented in Table S2. In total, 1893 edges connected the 229 differentially abundant taxa in the network. Many of the differentially abundant taxa correlated with each other positively (Fig. 4), indicating that these taxa likely respond similarly to changes associated with wetland degradation. Actinobacteria, Bacteroidetes, Chloroflexi, Euryarchaeota and Proteobacteria showed higher degree centrality in the network, suggesting that these taxa are strongly associated with the other members of the community (Fig. 5). Taxa belonging to Actinobacteria, Bacteroidetes, Chloroflexi, Proteobacteria and Planctomycetes had higher betweenness centrality compared to the other taxa, indicating that these taxa likely influence the interactions among taxa in the community (Fig. 6).
Based on greedy modularity optimization, three potential functional modules (M1, M2, and M3) were detected according to association strength among the differentially abundant taxa (Fig. 4). Module M1 included Jiangella (phylum Actinobacteria), Anaerolinea (phylum Chloroflexi), Desulfobulbus (Proteobacteria), Flavobacterium (Bacteroidetes), and Methanobacterium and Methanosaeta (phylum Euryarchaeota). The abundances of these taxa were significantly higher in the swamp soil than in the degraded soils (Fig. 6, Table 2). With the exception of DA101 in M2 and Acidovorax in M3, most of the differentially abundant taxa in modules M2 and M3 were not identified at genus level. The relative abundances of almost all taxa in each module were distinct in different soils (Fig. 6). The variation in abundance and taxonomic composition may provide clues on the functional modules that potentially rapidly and reversibly respond to environmental changes induced by wetland degradation through species sorting (plastic adjustment). Although the microhabitats in the meadow and sandy soils were distinct, the composition of the functional modules was similar, especially in modules M1 and M2 (Fig. 6).
PCoA result showed significant (p = 0.001) correlations between network modules and soil properties, with soil TK ranking highest ( Fig. S5A-C). This was supported by Mantel test, which showed high correlations between soil TK and these three potential functional modules, further indicating that soil TK drives structuring in these functional modules (Fig. S5D).

Discussion
We studied shifts in bacterial and archaeal community composition in three soil types associated with wetland degradation: pristine swamp soil, moderately degraded meadow soil and highly degraded sandy soil. As in previous studies based on denaturing gradient gel electrophoresis 16,20 , the bacterial and archaeal community compositions in the three soils differed significantly, but all were dominated by the phyla Proteobacteria, Acidobacteria and Chloroflexi. Proteobacteria are prevalent in various ecosystems and involved in many biogeochemical processes 21 . As in Röskeet al. 22 , Proteobacterial communities were mainly composed of Deltaproteobacteria and Betaproteobacteria, both classified as copiotrophs stimulated by increased nutrients in the environment 15 , which could partially explain why these two groups were more abundant in the swamp soil with higher organic matter content. Acidobacteria include many oligotrophs 23,24 , and their relative abundance was higher in the degraded soils that contained less organic matter, thus possibly favoring these oligotrophs. Interestingly, Chloroflexi that are known to survive on plant residues and grow even under nutrient limitation 25,26 were more abundant in the sandy soil which had less nutrients than the pristine swamp soil.
Understanding the environmental variables that shape soil microbial community structuring is one of the key goals of microbial ecology. Soil water content was previously found to be the major factor connected with Zoige wetland degradation, with significant impacts on the decrease of fungal community abundance and prokaryotic  20,27 . With the lowering of water table, some of the wetland plant species that are essential in swamp soil formation are being replaced by meadow vegetation 28 . Different vegetation types may release different types of root exudates: the exudates are significant sources of carbon, nitrogen and energy of the soil microbiota 29 . Different plant species also have diverse nutritional needs, promoting the development of particular soil microbiome depending on the growth substrate 30 . In our study, water content, organic carbon content, and the total and available nitrogen and phosphorus contents all decreased along with the wetland degradation, making them  valuable indicators for evaluating wetland degradation. Soil total potassium content was found to be an important factor shaping the bacterial and archaeal communities. Similarly, Pereira et al. 31 found that the bacterial diversity and community composition in the soil were highly affected by the total potassium content. Potassium is easily leached because of its mobility in soils, and thus directly affected by soil degradation.
The microbial community structure in the pristine swamp soil differed from those in the meadow and sandy soils. Consistent with our hypothesis, more taxa were enriched in the swamp soil, indicating that wetland degradation decreased the abundance of these bacterial and archaeal taxa. Methanogenic Euryarchaeota is mainly distributed in anoxic environments, for example in wetlands, paddy soils and marine sediments, and is involved in greenhouse gas production 32 . In this study, as expected, the abundances of Euryarchaeota were significantly higher in the pristine soils than in the degraded soils. Previous studies indicated that Zoige wetland was one of the hotspots of methane emissions in the Qinghai-Tibetan plateau, harboring a unique methanogen community 33 . In contrast, there were hardly any methanogens in the meadow soil in Zoige wetland 34 , probably due to the decreased soil water table and content that had lead to less anaerobic conditions. Although we did not measure  the oxygen content in the soils, soil water content in the degraded soils had decreased significantly and water table was lower. A change in water table and content could directly affect the oxic-anoxic interfaces, increase the thickness of the boreal peatland aerobic layer, and ultimately change the microbial community composition and abundances 35 . Moreover, low water tables affect soil organic carbon (SOC) by changing the vegetation cover and SOC decomposition rate which is positively correlated with the methanogenic archaea community 36 . At the genus level, the enrichment of Methanobacterium, Methanosarcina and Methanosaeta in the swamp soil indicate their potential roles as 'keystone species' in methanogenesis in this ecosystem.
Syntrophy, defined as a long-term stable relation of organisms which can be either beneficial or unbeneficial, plays a key role in many anaerobic ecosystems, even in systems where metabolism is sustained by the reduction of electron acceptors such as sulphate or nitrogen oxides 37 . Species related to Syntrophobacteraceae generally constitute consortia with methanogenic Archaea to syntrophically convert propionate to methane, which is a typical metabolic process during anaerobic degradation. It is likely that the relative abundance of Syntrophobacteraceae affiliated species in this study would similarly vary with the Methanogenic Euryarchaeota in response to wetland degradation. However, phylogenetic analysis of the community showed that contrary to the decreased abundance of the methanogenic Euryarchaeota, the abundances of members of the Syntrophobacteraceae actually increased. Syntrophobacteraceae are extremely diverse in their nutritional transformation capabilities, with some members being able to use a wide variety of electron donors and carbon sources syntrophically or in pure cultures, and in the presence or absence of sulfate as electron acceptor 38 . This metabolic peculiarity confers these microorganisms with adaptability to changing environmental conditions, and may explain the increase of Syntrophobacter-related species in the degraded meadow and sandy microcosms. DA101 can occupy various ecological niches and perform oligotrophic life strategy in soil 39 . Aeromicrobium are lignocellulose decomposers and play important roles in organic matter turnover in soil 40 . Bradyrhizobium are beneficial for plant growth due to their nitrogen-fixing ability. Candidatus Nitrososphaera thrives in organic-rich habitats 41 . The enrichment of these taxa in the degraded meadow and sandy soils suggest their subtle differences in habitat adaption, which also contributes differently to the carbon and nitrogen cycling in this specific ecosystem.
Network analysis provides understanding of potential interactions such as cooperation, competition and niche partitioning in bacterial and archaeal communities, and may identify keystone populations 42 . A node with high betweenness centrality value is in a core location in the network, and may affect other interactions in the community 43 . Keystone species with maximum betweenness centrality values are likely to serve as gatekeepers in ecosystem functions 44,45 . In our study, most archaeal keystone species in the networks were affiliated with Euryarchaeota that participate in wetland carbon cycling as discussed above, while most bacterial keystone species were from Actinobacteria, Proteobacteria and Bacteroidetes. Of all the keystone species detected in our study, Desulfobulbus are the only ones involved in the sulfur biogeochemical cycle, particularly in the reduction of sulfate to sulfide 46 . Geobacter are able to oxidize recalcitrant organic carbon 47 , and Flavobacterium are stimulated in organic matter rich environments 48 . The organic matter content in the swamp soil was higher than in the degraded soils, which may explain the significantly higher abundance of Flavobacterium in the swamp soil. Taken altogether, the keystone species may be pivotal determinants of the succession of bacterial and archaeal taxa during the wetland degradation.
In summary, wetland degradation changed soil properties, which further affected specific bacterial and archaeal taxa involved in nutrient cycling resulting in distinct association patterns. Proteobacteria, Acidobacteria, and Chloroflexi were the most abundant phyla. Soil total potassium strongly affected the shaping of bacterial and archaeal communities. Interpretation of the interaction networks in soil bacterial and archaeal community may open new avenues for understanding the ecological roles of bacterial and archaeal taxa in the degradation of the Zoige alpine wetland soil. The observed changes may serve as early warning signals of soil degradation in alpine wetlands.

Materials and Methods
Study area and soil sampling. The Zoige Plateau is in the northeast of Qinghai-Tibetan Plateau (101°30′E-103°30′E, 32°20′N-34°00′N) at an average altitude of 3500 m (Table S1). The mean annual precipitation in this area is 654 mm 49 . Eighty-five percent of precipitation is delivered from April to September. The mean annual air temperature is 1.1 °C 33 , with the average temperatures ranging from −10.1 °C in January to 10.7 °C in July. Swamp soil (SW), meadow soil (MD) and sandy soil (SD) in Zoige were sampled in August 2014. The swamp was covered with water with an annual average water table of 5.87 cm above ground. The vegetation was dominated by hydrophytes and sedge hydro-mesophytes (Carexmuliensis, Carexlasiocarpa and Carexmeyeriana). The meadow surface was in a humid state and the vegetation was dominated by mesophytes and hydro-mesophytes (Kobresiatibetica), and the water table was approximately 4.26 cm below ground. The sandy soil surface was continuously dry and had little or only few Psammophytes as vegetation cover 50 , and the water table was at least 21 cm below ground. The soil texture was peat soil in SW, sandy loam in soil MD, and sandy soil in SD. The parent materials of swamp and meadow soils were homogeneous silt clay and Triassic slate residues, and sandstones and siltstone, respectively 51 . In the sandy soil, the top layer (0-30 cm) was mostly sand and rest was composed of aeolian parent material 52 . Based on the FAO soil taxonomy system 53 , the SW soil is classified as Histosol, the MD soil as Haplic Calcisol, while the SD was Arenosol.
For each soil type, six sampling plots of at least 50 × 50 m 2 with a minimum distance of 100 m from each other were selected. At each sampling plot, three soil cores from the top 20 cm below the litter layer were collected using a 6.5 cm diameter steel push corer, and mixed into one approximately 1.5 kg homogenized composite sample. Visible plant roots and residues were removed before mixing. The corer was sterilized using 70% ethanol between samplings. Samples were stored in plastic bags on ice and divided into two portions. Sub-samples for physicochemical analyses and DNA extraction were stored at 4 °C and at −80 °C, respectively. Soil physicochemical properties. Soil pH and soil organic carbon (SOC) were determined using a soil-to-water ratio of 1:5 and dichromate oxidization method, respectively 54 . Soil gravimetric water content (WC) was determined by oven-drying at 105 °C for 48 h 54 . Soil total nitrogen (TN) and available nitrogen (AN), total phosphorus (TP) and available phosphorus (AP), total potassium (TK) and available potassium (AK) were determined by Kjeldahl digestion, alkaline hydrolysis diffusion method, and molybdenum blue method and flame photometry, respectively 54 .
DNA extraction, amplification and sequencing. DNA was extracted from 0.75-0.91 g fresh weight soils (corresponding to 0.50 g dry weight) using the FastDNA Spin Kit for Soil (MP Biomedicals, Solon, OH, USA) following the manufacturer's instructions. The quality of DNA was examined using the Nano-200 spectrophotometer (Aosheng, Hangzhou, China). DNA samples were stored at −20 °C.
PCR amplification was carried out using the primers 515F (5′-GTGCCAGCMGCCGCGGTAA-3′) and 806R (5′-GGACTACVSGGGTATCTAAT-3′) targeting the V4 hyper variable region in the bacterial and archaeal 16S rRNA gene 55 . The reverse primer was combined with the adapter and barcode sequences for multiplexing. Amplification was done in 50 μl-reaction volumes containing 3 U of TaKaRa Ex Taq HS (TaKARA Shuzo Co., Shiga, Japan), 5 mM dNTP Mixture (TaKARA), 2.0 mM MgCl 2 , 5 μl of 10× Ex Taq Buffer (TaKARA), 0.5 mM of each primer and 2.5 ng of soil DNA. The PCR procedure included an initial denaturation at 94 °C for 4 min, 30 cycles of 15S at 94 °C, 15S at 55 °C and 30S at 72 °C, and a final extension at 72 °C for 10 min.
PCR products were purified using PCR Clean-up Purification Kit (MP Biomedicals) and quantified using Qubit 2.0 fluorimeter (Invitrogen, Carlsbad, CA, USA). Purified amplicons were pooled in equimolar concentrations and sequenced using MiSeq Reagent Kit V2 for Illumina MiSeq. The raw sequence data were submitted to NCBI Sequence Read Archive (https://www.ncbi.nlm.nih.gov/sra/) with accession number SRS2127453.

Bioinformatics analysis.
Sequence reads were quality-filtered in QIIME (v1.9.1) 56 . Sequences with phred-quality score over 20 and longer than 300 bps were kept for downstream analyses while those could not be assembled were discarded. Chimeric sequences were removed using UCHIME against a reference alignment 55 . Sequences were clustered into operational taxonomic units (OTUs) at 97% using the UPARSE pipeline 57 . The RDP classifier was used to assign OTU representative sequences at 70% threshold 56 . For inter-group comparison, OTUs were randomly resampled at a depth of 12822 sequences. Singletons were eliminated during resampling, results in 35,014 OTUs. The phylogenetic tree was generated based on the representative sequences and used for UniFrac distance calculations. Statistical analysis. We applied principal component analysis (PCA) to assess differences in properties between soils. Analysis of similarity (ANOSIM) was used to test significant differences among abundances of taxa. The relative abundances of taxa at phylum and family levels were analyzed using R 58 . Two-way ANOVA was performed to test differences in the abundances of taxa at family level at P < 0.05. Principal coordinate analysis (PCoA) based on weighted and unweighted UniFrac distances were carried out in the R package vegan 59,60 . Distance based redundancy analysis (dbRDA) was performed to investigate factors driving community variation. Significance of factors was calculated by the permu test function in vegan over the dbRDA model using a maximum of 999 permutations. The goodness of fit for each edaphic factor was estimated by applying the envfit function in vegan (999 permutations). Mantel test was used to determine correlations between soil properties and microbial community composition.
The R package DESeq2 was used to analyze differential log 2 -fold transformed abundances at phylum and genus levels 61 . Differential abundance analysis was performed by fitting a generalized linear model with a negative binomial distribution to the normalized value for each type of taxa, and testing for differential abundance using a Wald test 61 . P-values were adjusted for multiple testing following the procedure of Benjamini and Hochberg 62 , and a false discovery rate (FDR) of 10% was selected to denote statistical significance 63 . Taxa with significant changes in abundance were defined as those that had differential abundance >1.0 and FDR-adjusted P < 0.1. Correlations between taxa with significant changes in abundance and edaphic features were determined using Spearman's rank correlation.