High‐throughput analysis of anammox bacteria in wetland and dryland soils along the altitudinal gradient in Qinghai–Tibet Plateau

Abstract This study investigated the diversity, community composition, and abundance of anaerobic ammonium oxidation (anammox) bacteria along the altitudinal gradient in Qinghai–Tibet Plateau. Two types of soil samples (wetland and dryland soils, n = 123) were collected from 641 m to 5,033 m altitudes. Polymerase chain reaction (PCR) screening showed that anammox were not widespread, and were only detected in 9 sampling sites of the 50 sites tested by amplifying the 16S rRNA genes. Then, only samples collected from Linzhi (2,715 m), Rikaze (4,030 m), and Naqu (5,011 m), which were positive for the presence of anammox, were further processed to explore the biogeography of anammox bacteria in Qinghai–Tibet Plateau. Results of high‐throughput sequencing targeting the hydrazine synthesis β‐subunit (hzsB) gene revealed the presence of three known anammox genera (Candidatus Brocadia, Candidatus Jettenia, and Candidatus Kuenenia) in both soil types. Their diversity, community composition, and abundance did not show significant variation with altitude at large scale. However, it was the small‐scale environmental heterogeneities between wetland and dryland soils that determined their biogeographical distribution. Specifically, the dryland soils had higher diversity of anammox bacteria than the wetland soils, but their abundance patterns varied. The community composition of anammox bacteria were found to be influenced by soil nitrate content.

physiological characterization of some representatives and led to the description of five genera, that is, Candidatus Kuenenia (Schmid et al., 2000), Candidatus Brocadia (Strous et al., 1999), Candidatus Scalindua (Damsté et al., 2003), Candidatus Anammoxoglobus (Kartal et al., 2007), Candidatus Jettenia (Quan et al., 2008), and Candidatus Anammoximicrobium (Khramenkov et al., 2013). Most knowledge on the diversity of these taxa were mainly based on environmental surveys of amplified 16S rRNA gene, and to date, has been the most widely used approach in detecting these uncultivable taxonomic groups (Schmid et al., 2005). However, the highly conserved nature of 16S rRNA limits its efficacy and sensitivity in examining the true diversity of anammox bacterial communities and their distribution in various environments (Hirsch, Long, & Song, 2011). Several studies have also shown that the commonly used 16S rRNA gene primers could not reveal the extent of diversity of these anammox bacteria Harhangi et al., 2012;Lipsewers et al., 2014).
Among these, the hzsB gene has been regarded as a good biomarker to determine the diversity and abundance of anammox bacteria (Wang et al., 2012;Zhu et al., 2013Zhu et al., , 2015, and has been widely used to study the occurrence and composition of anammox bacteria in various ecosystems (Bai, Chen, He, Shen, & Zhang, 2015;Xi, Ren, Zhang, & Fang, 2016;Yang et al., 2015;Zhu et al., 2015). In addition, the fast-paced developments in high-throughput sequencing make it feasible to obtain a more comprehensive understanding, especially in the field of microbial community ecology (Adetutu et al., 2012). However, the use of high-throughput sequencing approach to profile the anammox community based on hszB gene has not yet been carried out.
Anammox bacteria are ubiquitously distributed at redox transition zones in marine Kuypers et al., 2003), freshwater (Schubert et al., 2006;Zhang et al., 2007), and terrestrial ecosystems (Humbert et al., 2010;Shen et al., 2013). Their distribution, community composition, and abundance in these various ecosystems are affected by numerous environmental factors, such as sediment organic C/organic N (OrgC/OrgN), nitrite concentration, and sediment median grain size (Dang et al., 2010). In Pearl Estuary, for example, patterns of distribution of anammox bacteria were associated with salinity, temperature, and pH of the overlying water mass (Fu et al., 2015). Most studies conducted on the community of these taxa, however, were based on aquatic ecosystem, and to date, only few studies are available based on terrestrial ecosystems including the alpine ecosystems (Long, Heitman, Tobias, Philips, & Song, 2013;Xi et al., 2016).
These alpine ecosystems are quite unique as they receive low reactive nitrogen input, and understanding the role of the anammox bacteria in such conditions would have profound implications on the general understanding of nitrogen cycle. In particular, there is still a large gap existing in the study on anammox bacterial distribution in representative environments (wetland and dryland) of the alpine ecosystems.
Hence, the goals of this study are: (1) to explore the occurrence of anammox bacteria in alpine ecosystems and (2) to study their diversity, community composition, and abundance, and the crucial environmental factors that influence their biogeography. To achieve these, we collected samples from the Qinghai-Tibet Plateau, which is considered as the Earth's third polar region. Two representative habitats, a wetland and a dryland, situated along an altitudinal gradient were selected as the sampling sites. Community composition and diversity of anammox bacteria were determined by Illumina sequencing analysis targeting the hzsB gene. We then applied multivariate statistical analyses to link the community patterns with potential environmental factors influencing their biogeography in Qinghai-Tibet Plateau. Our findings contribute to the further understanding of the anammox process and the biogeography of anammox bacteria.

| Site description and soil sampling
The Qinghai-Tibet Plateau (26°00′-39°47′ N, 73°19′-104°47′ E) ranges from the southern Himalaya Mountains to northern Kunlun Mountains, Altun Mountains, and Qilian Mountains, bounded by the Pamir and Karakoram Mountains in the west, and Qinling Mountains and Loess Plateau in the east. It is the highest plateau on earth and known as "the roof of the world" with an average altitude of more than 4,000 m, covering 2.5 million·km 2 , nearly a quarter of the total land area of China. It is surrounded and traversed by several snowcapped mountain ranges and is the origin of many major rivers. It is characterized by strong solar radiation (3,000-6,000 MJ/m 2 ), but its temperature decreases with altitude, with an annual average temperature below 5°C, particularly −15°C to −2°C in winter and 8°C to 18°C in summer. The mean annual precipitation in the Qinghai-Tibet Plateau is 360 mm and varies significantly from the southeast to the northwest. Diverse soil types including river and lake wetlands, grassland and permafrost developed in this region along the altitude, and anthropogenic activities have rare impact on those above 5,000 m. All climatological information was obtained from China Meteorological Data Service Center (http://data.cma.cn/).  Table S1) were selected. Wetland and dryland are two typical ecosystems, and the characters of the soils in each ecosystem are also very different. Wetland ecosystems include marsh, swamp, mires, and aquatic ecosystems like river, lake, and stream (Keddy, 2010). But dryland is a relatively arid region with forest, grassland, and farmland as the most common types of land use. So we collected both two types of soils at each sampling region to understand the biogeography of anammox bacteria better. At each sampling site, after the surface dead plants were removed, soil samples were collected at 0-10 cm depth by using soil drill with three parallel samples. The soil samples were placed in sterile plastic bags, sealed, and transported to the laboratory with ice packs and protected from light. At each site, three parallel samples were mixed to form one composite sample. Each mixed sample was divided into two parts: one part was passed through a 2-mm sieve for physicochemical properties analysis, and the other part was frozen at −80°C for molecular analysis.

| Analytical procedures for soil environmental variables
The soil chemical and physical properties were measured according to Bao (2000). Briefly, soil ammonia nitrogen, nitrite nitrogen, and nitrate nitrogen concentrations were measured by the SEAL Auto Analyzer 3 HR (Seal Analytical, UK) after extraction with KCl (1:5 soil/2 mol/L KCl solution) with detection limits of 0.015 mg·kg −1 , 0.015 mg·kg −1 , and 0.03 mg·kg −1 , respectively. The soil pH was measured by the DELTA 320 pH Analyzer (Mettler Toledo, USA) in a suspension of 1:5 soil:water. The soil moisture contents were analyzed by oven drying 2 g of fresh soil at 108°C until a constant weight was achieved, and then the samples were placed into a microwave muffle furnace at 550°C for 5 hr to determine the total organic matter (TOM) with detection limit of 0.02 g·kg −1 . The total nitrogen (TN), total carbon (TC), and total sulfur (TS) concentrations were determined using a Vario EL III Analyzer (Elementar Analyses System GmbH, Germany) with detection limits of 0.05 mg·kg −1 , 0.2 mg·kg −1 , and 0.25 mg·kg −1 . Triplicates were run for quality assurance/quality control (QA/QC) for all above measurements.

| DNA extraction, polymerase chain reaction, and Illumina sequencing
About 0.33 g of freeze-dried soil of each sample was used for DNA extraction by following the manufacturer's protocol of FastDNA SPIN Kit for Soil (MP Biomedical, Solon, OH, USA). The extracted DNA was checked on 1% agarose gel and its concentration was determined using a Nanodrop ® ND-2000 ultraviolet-visible spectrophotometer (Thermo Fisher Scientific, Schwerte, Germany). The ratio of the absorbance at 260 nm and 280 nm were all about 1.8, which indicated that DNA with a good quality was obtained.

| Quantitative real-time PCR
SYBR Green I based real-time PCR assays were carried out in a volume of 20 μl, containing 10 μl of SYBR ® Premix Ex Taq™ (TAKARA, Dalian, China), 0.4 μl of each primer (10 pmol/μl), and 2 μl of 10-fold diluted DNA template, and was topped up with ddH 2 O to a total volume of 20 μl. Amplification and detection were carried out using an ABI Prism 7300 Sequence Detection System (Applied Biosystems, USA) with the primer sets and thermal profiles described as above.
Three no-template controls (NTCs) were run for each quantitative PCR assay. Ten-fold serial dilutions of a known copy number of the plasmid DNA were subjected to real-time PCR in triplicate to generate an external standard curve. Melting curves were generated after each assay to check the specificity of amplification. PCR efficiencies were 90%-103% (average 92%) and only the results with correlation coefficient above 0.98 were employed. The melt curve analyses were performed to confirm the specificity of PCR amplifications. All tests were performed in triplicate.
The detection limit of the hzsB gene was determined by amplifying the diluted plasmid DNA (10-fold dilution series) with the hzsB gene as insert. The sample would be assumed to be at the detection limit in three situations: (1) the last sample with C t value showed high standard deviation between replicates; or (2) the sample with C t value was more than that of the negative control (normally the C t of negative control was around 40 in this study) as shown in Figure S2; or (3) the melting curve was not run as a single peak. Using the hzsB gene on a plasmid, detection limit was around 1.34 × 10 2 copies judging from the amplified standard curve ( Figure S1). For the environmental samples like the soils in Qinghai-Tibet Plateau which can never be as clean as the plasmid DNA sample, hzsB gene copies would be no longer reliable when the number was lower than 1.34 × 10 3 copies·g −1 . This number was the detection limit used in this environmental investigation.

| Sequencing analysis
Sequencing reads were assigned to each sample according to the unique 6-bp barcode of each sample. The barcode and primers then were removed. Raw sequences from original DNA fragments were merged using FLASH (Magoč and Salzberg, 2011) and then were filtered (i.e., with a quality score <25 and read length <200 bp were filtered using the split_libraries command) using the QIIME software package (Caporaso et al., 2010). Then, the chimeric sequences were removed using UCHIME (Edgar, Haas, Clemente, Quince, & Knight, 2011). To accurately detect and correct frameshifts caused by indel sequencing errors, the FrameBot  tool was used.
Briefly, only the sequences containing no ambiguous bases (N), without any barcode or primer mismatches, and with the corrected frameshifts and length (about 292 bp) were included into the downstream analysis. The unique sequences were obtained by Mothur (Schloss et al., 2009) from the remaining high-quality nucleic acid sequences and then translated to protein sequences. Preprocessed sequences were clustered into operational taxonomic units (OTUs) based on their sequence similarity using UCLUST (identity = 0.90) (Edgar, 2010). A representative sequence for each OTUs was finally aligned using MUSCLE program (Edgar, 2004). A local alignment search was conducted with hzsB protein sequences using Basic Local Alignment Search Tool (BLAST).
Singletons were excluded and resampling according to the minimum sequence numbers across all samples was performed before calculation. A variety of alpha diversity indices including Chao1, Shannon, and Simpson were calculated by Mothur software. The Chao1 measured the richness of phylotypes, Shannon index estimated for both richness and evenness, whereas Simpson index detected dominant OTUs in the samples. The beta diversity was determined by weighted and unweighted UniFrac (Lozupone & Knight, 2005). Shifts in the bacterial community compositions were visualized using a principal coordinate analysis (PCoA) of the pairwise Bray-Curtis dissimilarity matrices of OTUs similarity across the different samples. The beta diversity was determined by additive partitioning approach (Lande, 1996) based on OTUs results.

| Occurrence and abundance of anammox bacteria
The 16S rRNA genes of Planctomycetes and anammox bacteria at all sampling sites in Qinghai-Tibet Plateau were specifically amplified.
Results of sequencing revealed that anammox-related 16S rRNA genes were only detected in 9 sampling sites of the 50 total sites, and were most closely affiliated to Candidatus Brocadia (Table S1) (Table S2).

| Diversity of anammox bacteria
A total of 75,959 raw reads generated through high-throughput pyrosequencing were obtained from the six samples collected from wetland and dryland soils. After quality filtering, a total of 72,099 high-quality sequences were obtained ranging from 7,541 to 16,502 sequences per sample. Clustering of the sequences at 90% similarity by UCLUST generated operational taxonomic units (OTUs) that ranged from 21 to 41 OTUs depending on the sample (Table 2).
Rarefaction curves of alpha diversity reached plateau indicating that the potential true diversity of anammox community had been well captured at the sequencing depth used in this study ( Figure S2). This was consistent with the high taxonomic similarity of the OTUs with known representative of anammox species, ranging from 0.99 to even 1.00 at 90% identity cutoff.
The Chao1, Shannon, and Simpson indices were calculated to estimate the alpha diversity of anammox bacteria in Qinghai-Tibet Plateau (Table 2). Alpha diversity was a comprehensive norm to represent community richness and diversity. The Chao1 index related positively with community richness, Shannon index was positively correlated with community diversity, while Simpson index correlated negatively with community diversity. The dryland soils showed higher Chao1 and Shannon values than those of the wetland soils, but exhibited a lower Simpson diversity. This indicated that the anammox bacterial community had a higher diversity in the dryland soils than the wetland soils. Furthermore, differences were also found between the three sampling sites regardless of the soil types, showing that the anammox bacterial community had the lowest diversity at the midaltitude site. Similar to the abundances, Spearman's rank correlation did not reveal significant correlations between alpha-diversity indices and the investigated soil properties (Table S3).
The beta diversity of anammox bacteria in Qinghai-Tibet Plateau was studied using additive partitioning approach (Lande, 1996).
Furthermore, the C beta , representing contribution to gamma diversity, was calculated to quantify and compare the beta diversity among various samples from different points (Crist & Veech, 2006) (

| Community composition of anammox bacteria
A local BLAST analysis of the amplified hzsB protein sequences was performed to assign identity of the sequences. Five protein sequences were used as references for BLAST (Table S6)  . The data as shown in Figure 4 indicated that the first RDA axes explained 54.1% of the variance in the anammox assemblages and the species-environmental correlation at the first RDA axes was 91.5%.

| DISCUSSION
The present study is the first report about the occurrence, diversity, Knowledge on the diversity of bacteria in ecosystems is key information needed for understanding the underlying mechanisms of global All of the indices were calculated based on OTUs results; γ = the total species richness in a certain region; α = the average species richness in sampling points in the region; β = γ − α; C beta = β/γ × 100%.
T A B L E 3 Beta diversity of anammox bacteria in Qinghai-Tibet Plateau nitrogen cycle (Philippot et al., 2013;Taroncher-Oldenburg, Griner, Francis, & Ward, 2003). Our study illustrated that the anammox bacteria had lower diversity in wetland than dryland soils, indicating distinctive spatial heterogeneity in alpine ecosystems between wetland and dryland. To date, the diversity of anammox bacteria has been studied in various ecosystems. In aquatic ecosystems, widespread occurrence but low diversity of anammox were shown in marine ecosystems . So far, the available anammox 16S rRNA sequences from marine (e.g., Black Sea and the Benguela OMZ, Namibia) and estuarine environments (e.g., Randers Fjord, Denmark) were all related to Candidatus Scalindua. Anammox diversity in freshwater ecosystem was higher than in marine, for example, the 16S rRNA genes in the sediments of Xinyi River were closely related to Candidatus Brocadia anammoxidans (Zhang et al., 2007) and in some freshwater extreme environments, most hzsB gene sequences were closely affiliated to Candidatus Kuenenia . Anammox diversity in terrestrial system was higher in soil than in freshwater environments, for example, in paddy soils, high alpha diversity (Shannon index = 1.84-2.71) was found and anammox bacteria is related to Candidatus Brocadia, Candidatus Kuenenia, and two novel unidentified clusters . In vegetable soils, the Shannon index was also high (2.04-2.59) and three different genera of anammox bacteria are detected, including Candidatus Kuenenia, Candidatus Brocadia, and Candidatus Jettenia (Shen, Wu, Liu, & Li, 2017). The above discussions indicated in natural ecosystems, the highest diversity of anammox bacteria occurred in terrestrial systems, followed by freshwater systems and marine systems, which was in accordance with our finding that anammox bacteria had higher diversity in dryland than wetland of alpine system.  the dominant species in the anammox community (Suto et al., 2017;Wang, Peng, Ma, Wang, & Zhu, 2015). This observation suggests environmental selection of anammox bacteria in natural and engineered ecosystems. However, in the present study, it was not found that the anammox bacteria showed a variation with altitude at large scale, and anammox bacteria community showed a very strong heterogeneity between different points in each altitude. This indicates small-scale environmental heterogeneities are important in shaping the community composition and abundance of anammox bacteria.
Certain habitat-specific studies have shown that anammox bacterial distribution is affected by different environmental factors, for example, temperature (Hou et al., 2013), and available nitrite and ammonium concentrations (Li & Gu, 2013). A recent study using global data ordination demonstrated that salinity is one of the key factors driving the biogeography of the anammox bacteria (Sonthiphand et al., 2014 (Wr et al., 2008). In this study, RDA and Spearman's correlation analysis indicated that the community composition of anammox bacteria significantly correlated with the concentrations of NO 3 − . Particularly, NO 3 − was correlated well with the relative abundance of Candidatus Kuenenia (p < .05). However, the role of other unidentified ecological parameters cannot be ruled out (Shehzad et al., 2016).
The Qinghai-Tibet Plateau received low load of reactive nitrogen. In addition, the surveyed regions were perennially in low temperature. In the year we conducted the study (2015), the temperature was −12°C to 15°C in 8 months of the year (http://www.tianqi.com/), while anammox bacteria were more abundant at higher temperatures, that is, their optimum temperature in temperate shelf sediments was 15°C . As a result, the growth of anammox bacteria was limited in such adverse conditions causing the observed low abundance of anammox bacteria in Qinghai-Tibet Plateau.
In this study, RDA and Spearman's correlation analysis indi-