Metagenomic Insight Into Patterns and Mechanism of Nitrogen Cycle During Biocrust Succession

Wang, Qiong; Han, Yingchun; Lan, Shubin; Hu, Chunxiang

doi:10.3389/fmicb.2021.633428

ORIGINAL RESEARCH article

Front. Microbiol., 17 March 2021

Sec. Terrestrial Microbiology

Volume 12 - 2021 | https://doi.org/10.3389/fmicb.2021.633428

Metagenomic Insight Into Patterns and Mechanism of Nitrogen Cycle During Biocrust Succession

$\r\nQiong Wang,$ Qiong Wang^1,2

Yingchun Han^1,2

Shubin Lan¹

Chunxiang Hu^1*

¹Key Laboratory of Algal Biology, Institute of Hydrobiology, Chinese Academy of Sciences, Wuhan, China
²University of Chinese Academy of Sciences, Beijing, China

The successional ecology of nitrogen cycling in biocrusts and the linkages to ecosystem processes remains unclear. To explore this, four successional stages of natural biocrust with five batches of repeated sampling and three developmental stages of simulated biocrust were studied using relative and absolute quantified multi-omics methods. A consistent pattern across all biocrust was found where ammonium assimilation, mineralization, dissimilatory nitrite to ammonium (DNiRA), and assimilatory nitrate to ammonium were abundant, while denitrification medium, N-fixation, and ammonia oxidation were low. Mathematic analysis showed that the nitrogen cycle in natural biocrust was driven by dissolved organic N and NO₃^–. pH and SO₄^2– were the strongest variables affecting denitrification, while C:(N:P) was the strongest variable affecting N-fixation, DNiRA, nitrite oxidation, and dissimilatory nitrate to nitrite. Furthermore, N-fixation and DNiRA were closely related to elemental stoichiometry and redox balance, while assimilatory nitrite to ammonium (ANiRA) and mineralization were related to hydrological cycles. Together with the absolute quantification and network models, our results suggest that responsive ANiRA and mineralization decreased during biocrust succession; whereas central respiratory DNiRA, the final step of denitrification, and the complexity and interaction of the whole nitrogen cycle network increased. Therefore, our study stresses the changing environmental functions in the biocrust N-cycle, which are succession-dependent.

Introduction

The nitrogen cycle is not only a major element cycle related to climate change but also a burgeoning field of research for the global materials cycle (Delgado-Baquerizo et al., 2013; Weber et al., 2015; Keren et al., 2020; Li et al., 2020). Particularly, our understanding of the nitrogen cycle due to breakthrough discoveries of new microorganisms and techniques has been remarkable (Reverey et al., 2016; Stein and Klotz, 2016; Kuypers et al., 2018; Zehr and Capone, 2020). The existence of new metabolic pathways is the result of the balance and distribution of nitrogen in ecosystems (Rütting et al., 2011) and directly related to the stability of ecosystem function (Isobe and Ohte, 2014). However, the distributional patterns of each pathway and their links to ecological processes remain unclear (Isobe and Ohte, 2014).

The known nitrogen cycle involves 9 inorganic redox states and more than 14 transformation pathways (Stein and Klotz, 2016; Kuypers et al., 2018), although they are often habitat dependent (Lüke et al., 2016; Stein and Klotz, 2016; Li et al., 2018; Wang et al., 2020). Our understanding of the nitrogen cycle remains mostly at the chemical and individual organism level: from the classical view including nitrogen fixation, nitrification and denitrification, to the classification of assimilation [ammonium assimilation, nitrogen fixation and assimilatory nitrate to ammonium (ANRA)] and dissimilation [mineralization, ammonium oxidation, nitrite oxidation, denitrification and dissimilatory nitrate reduction to ammonium (DNRA)] (Thamdrup, 2012), to the division of ammonification, nitrification, denitrification, anammox and NO₂^–-NO₃^– interconversion (Stein and Klotz, 2016). At the community level, although inferences have been made related to ecological processes (Rütting et al., 2011; Isobe and Ohte, 2014), there is limited knowledge from the perspective of ecology because of unknown nitrogen abundance and correlations in the nitrogen transformation processes. Environmental factors such as redox potential (Eh), pH, NH₄⁺, NO₃^–, and C:(N:P) can affect different transformation pathways; however, most studies focus on a single pathway (Pett-Ridge et al., 2006; Burgin and Hamilton, 2007; Rütting et al., 2011; Sgouridis et al., 2011; Bouskill et al., 2012; Jiang et al., 2015; Zheng et al., 2020), but how these factors drive overall patterns in N-cycle is poorly understood.

Biological soil crusts (biocrusts) are the biological-soil mosaic layers on dryland surface and consist of Cyanobacteria, lichens, mosses, algae, and heterotrophic organisms in varying proportions. Biocrusts have been recognized as an ideal model system for soil eco-function research, where different types of biocrusts represent different developmental and successional stages, which significantly affect the stability and eco-functions of the local soil ecosystem (Bowker et al., 2014; Ouyang and Hu, 2017). Similar to other soil types, N-cycling in biocrusts is susceptible to precipitation or water (Zaady, 2005; McCalley and Sparks, 2008; Gu and Riley, 2010; Delgado-Baquerizo et al., 2014) and biocrusts are often in a nitrogen-limited state under drought conditions (Hooper and Johnson, 1999; Schulz et al., 2013). In response to episodic water, the N-cycle in biocrusts occurs in pulses (Austin et al., 2004; Collins et al., 2008), with N₂O, NO, and HONO emissions (Abed et al., 2013; Weber et al., 2015). Although measurements of process rates are limited and inconsistent, typically, nitrogen fixation and ammonia oxidation are higher (Zaady, 2005; Johnson et al., 2007; Strauss et al., 2011), but occasionally denitrification is stronger (Abed et al., 2013). It is likely that no matter the soil or biocrust, they all perform differently depending on the intensity and frequency of water events (Reverey et al., 2016), leading to the inconsistencies in different N-cycle studies. Biocrusts under different successional stages may have distinct microbial community structures and redox states, thus, the potential nitrogen transformation rates measured under similar experimental conditions are difficult to compare, and may not reflect biological reality.

The soil N-cycle is sensitive to water availability, whereas the patterns in metagenomic sequencing are relatively stable (Nelson et al., 2016). We speculate that the stable N-cycle genetic potential of different biocrusts can be obtained through repeated sampling, with distinguishable phenotypic and relatively stable eco-physiological characteristics of biocrusts (Grote et al., 2010; Lan et al., 2017; Maier et al., 2018). We can determine the relevance of different N pathways by coupling them with gene abundance relationships. Furthermore, the relationship between environmental variables and N pathways will help us understand the driving mechanisms of biocrust succession and relevant changes in ecological functions. In this study, four types of natural biocrusts representing successional stages with Cyanobacteria being replaced by lichens or moss, were repeatedly sampled, and three types of simulated biocrusts in which Cyanobacteria gradually become the dominant microorganism were also analyzed as controls, in order to explore the successional variability of N-cycle and its linkage mechanism to ecosystem functions. We hypothesized that the N-cycle which included diverse transformation processes in biocrusts is not a simple balance of nitrogen fixation and denitrification, in the environment of biocrusts with similar pH (Lan et al., 2015), those N-transformation processes will be significantly regulated by the variables closely related to redox condition, and the ability of the microbial community to coordinate the nitrogen input, loss, retention, recycle and utilization will improve with succession.

Materials and Methods

Sites Description and Sampling

Four types of natural biocrusts (algae crusts, A; cyanolichen crusts, C; chlorolichen crusts, G; and moss crusts, M) were collected in 2015, 2016, 2017, and 2018 from the Shapotou area located near the southeastern edge of the Tengger Desert (37° 32′ N, 105° 02′ E), China, by carefully sampling pieces of crusts in their natural thickness (4–15 mm; without the soil beneath) with a sharp shovel as described previously (Lan et al., 2015). Among them, algae crusts, lichen crusts and moss crusts represent the early, transitional, and later stages of succession, respectively, and the two types of lichen crusts have distinct survival strategies (Lange et al., 1998; Grote et al., 2010; Lan et al., 2017). Simulated crusts were cultured in situ in September and October 2017 at the east edge of the Hobq Desert (40° 21′ N, 109° 51′ E) (Ouyang et al., 2017). Microcoleus sp. isolated from biocrusts of Shapotou were inoculated onto the shifting sandy surface in September 2017, watered twice daily (08:00 and 18:00) during the first week only, and collected on the 11th, 35th, and 52nd day post-inoculation to represent early (Ae), middle (Am), and later (Al) stages of development, respectively. For 2015 natural biocrusts, the analysis of GeoChip and sequencing of 16S rRNA gene and N-cycle genes were completed in four replicates, while all other analyses (including metagenomics) were completed in triplicate. Natural biocrusts in 2016 and 2017 have only one replicate, but cyanolichen crusts included two types. Two batches of natural biocrusts were sampled in 2018 with four replicates for each type. For simulated biocrusts, there were four replicates for the three stages. A schematic of the experimental and analytical workflow is provided in Supplementary Figure S1.

Soil Analyses and Characterization

Immediately after the samples reached the laboratory, the water content (WC) (Zhang et al., 2018) and texture (Olmstead et al., 1930) were measured. Electrical conductivity, pH, and Eh were measured in a soil suspension with a soil:water ratio of 1:5 (w/v). Total carbon, organic carbon and dissolved organic carbon (DOC) were determined by carbon analyzer (vario TOC, Elementar, Germany). Total nitrogen and total phosphorus content were analyzed using spectrophotometry (Lan et al., 2015). Biomass carbon was analyzed via chloroform fumigation (Vance et al., 1987), and remnant carbon was calculated by subtracting biomass carbon from organic carbon. The contents of CO₃^2– and HCO₃^– were measured by double indicator titrations, that of NO₃^–, NO₂^–, NH₄⁺, PO₄^3–, Ca²⁺, and Mg²⁺ by ion chromatography. Dissolved organic nitrogen (DON) was calculated as the difference between total N and inorganic N (Zhou et al., 2020).

DNA Extraction and qPCR

DNA was extracted using the PowerSoil DNA Isolation Kit (MO BIO Laboratories, Carlsbad, United States), and isolates were stored at −80°C. DNA was quantified and examined for purity with NanoDrop 2000c (Thermo Fisher Scientific, United States). Absolute copies of bacterial, fungal, archaeal, and cyanobacterial rRNA genes, and nifH and nrfA genes were quantified by qPCR (Applied Biosystems StepOnePlus^TM Real-Time PCR system). The qPCR conditions and primers are given in Supplementary Table S1. All the qPCR reactions were run in triplicate.

Illumina MiSeq Sequencing of 16S rRNA Gene and N-Cycle Genes

Sequences associated with nitrogen fixation (nifH), denitrification (nirK), and bacterial 16S rRNA were amplified using the following primer pairs: nifH_F/nifH_R (Rosch et al., 2002), nirK_F1aCu/nirK_R3Cu (Throback et al., 2004), 338F/806R (Castrillo et al., 2017). DNA libraries were sequenced on the Illumina MiSeq platform with 2 × 250 pair-ends. For 16S rRNA analysis (Supplementary Table S2), sequences of chloroplasts were removed. All quality-filtered sequences were clustered into OTUs at 97% identity and randomly resampled to the depth of minimum number of sample sequences. Representative sequences of 16S rRNA gene were annotated taxonomically against the SILVA 138 database with a confidence cutoff at 0.7. Other sequences were classified by Blastn against the non-redundant protein database¹.

Metagenomic Sequencing and Assembly

Metagenomic libraries were sequenced on Illumina platforms (Supplementary Table S2). Raw sequence data quality was assessed using FastQC², trimmed using Seqprep³, and high quality reads were extracted by filtering low-quality reads with “N” base, reads with < Q20, and short reads of < 50 bp using Sikle⁴. De novo assembly via de Brujin graph approach was performed using the multiple mixing assembly strategy with IDBA-UD (Peng et al., 2012) (k-mer = 47–97) and Newbler⁵. Contigs (≥300 bp) were analyzed for open reading frames (ORFs) prediction using MetaGene (Noguchi et al., 2006). Non-redundant gene catalog was constructed with predicted ORFs (≥100 bp) using CD-HIT (Li and Godzik, 2006) at 95% identity and 90% coverage. All high-quality reads were aligned (95% identity) against the gene catalog via SOAPaligner (Li et al., 2008) to obtain the gene abundance in each sample. The abundance of each unigene was normalized with gene length as previously described (Qin et al., 2012; Zheng et al., 2016).

Metagenome Annotation

Taxonomic analysis was conducted by searching the putative amino acid sequences against NR databases (June 2018) using BLASTp (BLAST version 2.3.0; E-value ≤ 1e-5). To analyze the microbial nitrogen cycle, and minimize the influence of bryophytes whose genome size is much larger than microbes, the sequences that were “Bryopsida” at the class level were removed, and the remaining sequences were used for analysis of subsequent functional annotations using BLASTp against KEGG (version 59). The corresponding relationship between all the N-cycle associated genes and KEGG Orthology number is shown in Supplementary Table S3. Each gene combination or pathway was summed (multiple pathways performing the same conversion step) or averaged (multiple enzymes/subunits in the same conversion step) and normalized to 1,000,000 (Llorens-Mares et al., 2015).

GeoChip Analysis

GeoChip 5.0 (180K) analysis was performed as described previously (He et al., 2007; Yang et al., 2013) with little modification. Probes were considered positive if their signals were detected in at least 2 of 4 replicate sets. The signal intensity of each gene or functional category was divided by the number of probes included in this gene or functional category to reduce the impact of the number of probes.

Potential Activity of N-Cycle Pathway

Potential N fixation activity was measured by acetylene reduction assays (Belnap, 2002; Strauss et al., 2011). Biocrusts were evenly moistened with glucose-containing solution and incubated in a 100 mL serum bottle at 26°C for 8 h. Potential ammonia oxidation activity was determined according to methods described previously (Johnson et al., 2005; Strauss et al., 2011). Biocrusts were mixed with 1 mM (NH₄)₂SO₄ and 75 mM NaClO₃ and incubated on a shaker at 24°C for 10 h. The change in NO₂^– concentration between 2 and 10 h was measured. The acetylene inhibition assay was used to determine the potential denitrification activity (Strauss et al., 2011). Biocrusts and sterile denitrification media (containing 100 mg/L KNO₃ and 100 mg/L glucose) were mixed in a bottle, purged with nitrogen, and incubated at 24°C in the dark. The increase in N₂O concentration after 4 and 8 h was measured.

Statistical Analysis

The significant differences between various indicators in the samples were determined with SPSS 20.0 software (IBM, United States) using one-way analysis of variance (ANOVA) at p < 0.05. Non-metric multidimensional scaling analysis (NMDS) and redundancy analysis (RDA) were conducted using Canoco 5 (Šmilauer and Lepš, 2014). Analysis of similarity (ANOSIM) were conducted to test the statistical significance of differences between different types of biocrusts with vegan package in R, performing 10,000 permutations. Spearman’s correlation was calculated using psych package in R 3.4.4. Network diagrams were visualized with Cytoscape (Shannon et al., 2003). The co-occurrence network of the N-cycle key genes and high abundance genus (average abundance for each type in each batch above 0.01%, represent the top 833 of all 5,998 detected genus) was clustered with an edge-weighted spring-embedded algorithm with correlation coefficient r score. Each connection stands for a significant (p < 0.001) correlation. ModuLand as a Cytoscape plug-in (Kovacs et al., 2010) was used for network modularization, with a threshold of merging modules at 0.9. The network dissimilarity (β_w) between different successional stages was determined using previously described methods (Koleff et al., 2003; Liu et al., 2019). Details of methods can be found in Supplementary Material.

Results

Microbial Community Structure

The NMDS analysis of 16S rRNA gene sequences (Figure 1A) at the OTU-level showed that bacterial community structure had significant differences among different successional stages of natural biocrusts and also different developmental stages of simulated biocrusts (ANOSIM, R = 0.539, P = 0.001). The diversity index (S_obs or Shannon) increased with the succession of natural biocrusts but decreased throughout the development of simulated biocrusts (Figure 1B). For community composition, in both 16S rRNA gene sequencing (Figure 1C) and metagenomic sequencing (Figure 1D), the top three most abundant taxa were Cyanobacteria, Actinobacteria, and Alphaproteobacteria in both early (A) and transitional stages (C, G) of natural biocrusts, as well as middle (Am) and late stages (Al) of simulated biocrusts. With the succession of natural biocrusts, the relative abundance of Cyanobacteria decreased, while Acidobacteria and Gammaproteobacteria increased, but the trend was reversed through the development of simulated biocrusts. Thus, the abundance of Acidobacteria was higher than that of Cyanobacteria in late stage (M) natural biocrusts, and the abundance of Bacteroides was higher than that of Cyanobacteria in early stage (Ae) of simulated biocrusts.

FIGURE 1

Figure 1. Community structure of biocrusts. (A) Non-metric multidimensional scaling (NMDS) combined with analysis of similarity (ANOSIM) for bacterial community dissimilarity based on 16S rRNA gene sequencing (OTU level). Colors represent different stages of biocrusts; shapes represent sampling batches. (B) Alpha diversity (S_obs and Shannon) of bacterial communities based on 16S rRNA amplicon data set compared across different stages of biocrusts. (C) Mean relative abundance of primary bacterial phyla (Proteobacteria in class) with 16S rRNA gene sequencing. (D) Taxonomic analyses based on the metagenomic data on phyla level (Proteobacteria in class). A, C, G, M represent four types of natural biocrusts at different succession stages. A, algae crusts, n = 13; C, cyanolichen crusts, n = 15; G, chlorolichen crusts, n = 13; M, moss crusts, n = 13. Ae, Am, and Al, early, middle, and later stage of simulated cyanobacterial crusts, respectively.

The Abundance of N-Cycle Pathways

From metagenomic sequencing data, pathways involved in the N-cycle can be indicated by relevant genes. Based on metagenomic sequences, NMDS and ANOSIM analyses showed that, the N-cycle was significantly different among different successional stages of natural biocrusts, early (Ae), and middle-late (Am, Al) simulated biocrusts (R = 0.775, P = 0.001) (Figure 2A). For nitrogen transformation processes in natural biocrusts (Figure 2B), the relative abundance of ammonium assimilation (gdh, glnA, glt) was the highest, followed by mineralization (ureABC, gls), DNiRA(nirBD, nrfAH) and ANRA (nasAB, narB, NR, nirA, NIT-6), while that of DNRN (narGHI, napAB), nitrite oxidation (nxrAB) and denitrification (nirK, nirS, norBC, nosZ) were lower; and nitrogen fixation (nifDHK) and ammonia oxidation (amoABC) were the lowest. The abundance of genes for denitrification (nirK, nirS, norBC, nosZ) decreased along with the reaction sequences (Figure 2B). Except for cyanolichen crusts which have a special survival strategy, with the succession of natural biocrusts, the abundances of nrfAH and nosZ increased, while other genes such as glt, ureABC, gls, amoABC, hao, narB, nirA, nirBD showed a decreasing trend. The pattern of N-cycle in simulated biocrusts was similar to that of natural biocrusts, with a slightly higher relative abundance of anammox (hzsABC, hdh), nirBD in DNiRA, and nasAB, narB and nirA in ANRA. Along with the development of simulated biocrusts, the abundances of some genes such as narB, nirA, NR, ureABC, gls and glt increased, while the abundances of some other genes such as amoABC, hao, nxrAB, nirBD, nrfAH, nirK, nosZ, narGHI, glnA, and gdh decreased.

FIGURE 2

Figure 2. Nitrogen cycle based on metagenomic sequencing analysis. (A) NMDS and ANOSIM showing the significant difference of N-cycle for different stages of natural biocrusts (nBSCs) or simulated biocrusts (sBSCs) based on KEGG Orthology level. (B) Relative abundance of N-cycle genes or pathways. Different colors represent different nitrogen transformation pathways, and the values after pathways or genes are the mean relative abundance of nBSCs (5 batches) and sBSCs. Arrow width is proportional to the abundance of the nitrogen pathways. Heatmap from left to right indicate the change of pathways or genes with the succession of nBSCs (blue) or the development of sBSCs (yellow), and different letters represent significant differences within nBSCs or sBSCs (Tukey HSD; p < 0.05). Denitrification here indicates the transformation of NO₂^– to N₂. DNRN, dissimilatory nitrate to nitrite; DNiRA, dissimilatory nitrite to ammonium; ANRN, assimilatory nitrate to nitrite; ANiRA, assimilatory nitrite to ammonium; ND, not detected. DON, dissolved organic N.

From the absolute quantitative analysis using GeoChip (Figure 3), the normalized average signal intensity of amoA (amoA, amoA2), nirA and ureC gradually decreased along the succession of natural biocrusts, while nrfA and nosZ showed the opposite trend. With the development of simulated biocrusts, an increased signal intensity was observed for genes nirA, narB, and ureC, while a decreasing trend was observed for nrfA and p450nor. In addition, genes involved in anammox (hzsA and hzo) were detected both in natural biocrusts and simulated biocrusts.

FIGURE 3

Figure 3. Normalized average signal intensity of key genes in nitrogen cycle for natural biocrusts in 2015 (A) and simulated biocrusts (B). Error bars represent standard errors, and different letters represent significant difference (Tukey HSD; p < 0.05) within natural biocrusts (n = 4) or simulated biocrusts (n = 4). amoA2, glnA’ represent the gene names of amoa_quasi and glnA_fungi in GeoChip, respectively.

From quantitative PCR analysis (Figures 4A,B), we confirmed that gene copies of nifH and nrfA were highest in cyanolichen crusts. Potential nitrogen fixation and ammonia oxidation rates were both highest in cyanolichen crusts, and potential denitrification rates were highest in moss crusts (Figures 4C–E).

FIGURE 4

Figure 4. The gene copies of nifH (A) and nrfA (B), and the potential rate of classical nitrogen cycle including N₂ fixation (C), ammonia oxidation (AO, D), and denitrification (E) for natural biocrusts in 2015. Error bars represent standard errors, and different letters represent significant difference (Tukey HSD; p < 0.05) within natural biocrusts (n = 4) or simulated biocrusts (n = 4).

The Potential Microbial Contribution to N-Cycle

In natural biocrusts (Figure 5A), marker gene related to nitrogen fixation (nifH) mostly belonged to Cyanobacteria in early (A) and transitional stages (C, G), while Alphaproteobacteria dominated later stages (M). Main taxa involved in ammonium assimilation (gdh, glnA) and mineralization (ureC) were Alphaproteobacteria and Actinobacteria. Ammonia oxidation (amoA) was mostly carried out by Archaea. While Alpha-, Beta-, Delta-proteobacteria, and Actinobacteria were the main taxa contributing to interconversion of NO₂^–-NO₃^– (nxrA, narG, napA); Alpha-, Delta-proteobacteria, Verrucomicrobia, Actinobacteria, Cyanobacteria, and Planctomycetes were the main potential taxa involved in DNiRA (nrfA, nirB). For both assimilatory nitrate to nitrite (ANRN) (nasA, narB, NR) and ANiRA (nirA, NIT-6), microbes such as Actinobacteria, Cyanobacteria, Ascomycota, Planctomycetes, and Bacteroidetes were found. Potential denitrifiers were different depending on the reaction steps, but mainly belonging to Proteobacteria (Alpha-, Gamma-, Delta-), Actinobacteria, Cyanobacteria, Chloroflexi, Gemmatimonadetes, or Bacteroidetes.

FIGURE 5

Figure 5. Relative abundance of the microbial taxa related to each key gene of nitrogen cycle for natural biocrusts in 2015 (A) and simulated biocrusts (B).

In simulated biocrusts (Figure 5B), sequences matching to N-cycle were identified as similar organisms to those in natural biocrusts with a few differences: (1) potential nitrogen fixers were mainly Alpha- and Beta-proteobacteria; (2) as well as Archaea, Betaproteobacteria was also important for ammonia oxidation; (3) there were more Nitrospirae containing gene nxrA or narG; (4) Alphaproteobacteria were the most important bacteria involved in each step of denitrification; (5) Anammox always showed up in Planctomycetes.

Co-occurrence Network of Microbial Community and the N-Cycle

Based on the correlation (p < 0.001) of highly abundant genera and N-cycle marker genes, different networks were drawn for the four successional stages of natural biocrusts (Figure 6). Pairwise dissimilarity increased throughout biocrust succession. In the early (A) and transitional stages (C, G), the number of positive and negative edges, average degree, and network heterogeneity were lower; while they were higher in the late stages (M). For algae crusts (Figure 6; A), genes related to denitrification (nirK, norB, nosZ) and genera from several phylum formed a functionally independent module, as did the ammonium assimilation gene (glnA) and genera mainly from Proteobacteria; communities dominated by Cyanobacteria associated with genes such as nifH and nrfA. For cyanolichen crusts (Figure 6; C), there were more functional modules, such as the group dominated by Cyanobacteria and Alphaproteobacteria and associated with genes such as nifH, nrfA; the group coupled Alphaproteobacteria with NR and NIT-6. For chlorolichen crusts (Figure 6; G), there were also more modules, such as genera mainly from Alphaproteobacteria clustered with nrfA and nosZ; nifH, nirB and many other genes, clustered together. Communities dominated by Ascomycota and Alphaproteobacteria grouped together with NR, NIT-6. For moss crusts (Figure 6; M), nrfA and genera from several phyla formed the largest independent module. Another group was dominated by Actinobacteria and associated with genes such as nifH, nirB. With the succession of natural biocrusts, network complexity increased, and the associated functions such as denitrification, DNiRA, glnA, and nifH showed obvious shifts.

FIGURE 6

Figure 6. Co-occurrence network of the N-cycle key genes and high abundance of species for different successional stages of natural biocrusts (n = 13 for each stage). The colors of nodes represent different module, which are in line with the colors of abscissa label of histogram on the right (M stands for “module”). The community compositions of modules that contain few nodes are not shown in histogram. The size of the nodes corresponds to the degree of nodes. The connection stands for a significant spearman correlation (p < 0.001).

Relevance of Pathways and Influence of Environmental Variables

To reveal the relevance of N-cycle pathways, correlation networks (Figure 7A) were used. We found a strong positive correlation between DNRA and nitrite oxidation (R = 0.896, P < 0.001). A similar relationship was also observed between ANiRA and mineralization (R = 0.864, P < 0.001). As the hub of the network, DNRN and nitrite oxidation with the highest degree were significantly positively correlated with the most pathways. ANiRA and mineralization were positively correlated with nitrogen fixation and ANRN; while nitrogen fixation was negatively correlated with DNiRA, nitrite oxidation, and DNRN.

FIGURE 7

Figure 7. The spearman correlations between N-transformation pathways (A), as well as between environmental factors and N-transformation pathways (B), were depicted using Cytoscape. Environmental factors that shaped nitrogen cycle significantly (p < 0.05) are shown on pathway level (C) or gene level (D) for natural biocrusts with forward selection of RDA. For network (A,B), the connection stands for a significant spearman correlation (p < 0.001), and the edge width is proportional to the magnitude of correlation coefficient (|r|). EC, electrical conductivity; WC, water content; Eh, redox potential; DON, dissolved organic nitrogen; DOC, dissolved organic carbon; BC, biomass carbon.

For the environmental factors that influenced each N-cycle pathway, the correlation network (Figure 7B) showed that, NH₄⁺ was significantly and positively correlated with nitrogen fixation, ANiRA, and mineralization. In contrast, the remaining variables [C:(N:P), DON, C:N, WC, and DOC] were positively correlated with nitrogen fixation, but negatively correlated to the other six intermediate pathways (denitrification, DNiRA, DNRN, nitrite oxidation, ANRN and ammonium assimilation). Compared with other variables, C:(N:P) had the strongest correlation with DNiRA, followed by DNRN, nitrite oxidation, and ANRN, while pH was negatively correlated with nitrogen fixation and positively correlated with denitrification, followed by DNiRA, DNRN, nitrite oxidation, and ANRN. In addition, Eh, NO₃^–, and HCO₃^– were negatively correlated with both ANiRA and mineralization, while SO₄^2– was most strongly negatively correlated with denitrification. In short, both pH and SO₄^2– were the most important variables affecting denitrification. However, pH and C:(N:P) were not only the most important variables affecting nitrogen fixation but also the most important variables affecting DNiRA, nitrite oxidation, and DNRN.

Furthermore, to explore what environmental factors shaped the whole N-cycle function in natural biocrusts, RDA was used and showed that DON and NO₃^– were significant variables explaining the variance of the N-cycle in natural biocrusts (p < 0.05), either at the pathway or gene level (Figures 7C,D). DON was positively correlated with the relative abundance of nrfA, but negatively correlated to some other genes, such as nirB, nxrA, and narG; while NO₃^– was positively correlated with the relative abundance of nasA and negatively correlated with the relative abundances of genes, such as nirA and norB (Figure 7D).

Discussion

In our study, we investigated microbial community structure using 16S rRNA (Figure 1C), nifH or nirK (Supplementary Figure S2) gene amplicon sequencing and metagenomic sequencing (Figures 1D, 5A). Community structures obtained from both methods were similar, although we cannot rule out issues caused by the limitation of primers and unknown microbes. We analyzed 5 batches of data from 4 different years (Figure 2) and found that the relative abundances of genes involved in denitrification decreased as reaction steps proceeded (nirK > norBC > nosZ), further supporting decreased growth yields per molar substrate (Koike and Hattori, 1975), as well as increases in oxygen sensitivity of denitrification related enzyme activities (Kraft et al., 2011). Although changing sequencing platforms and multiple batches may have biased our results, the differences between different successional stages are quite large (Figure 2A), confirming the stability and reliability of our outcomes. Thus, the genetic characteristics of the four types of natural biocrusts in the present study represent the different successional stages very well and fully meet our experimental expectations. Furthermore, in both natural and simulated biocrusts, we confirmed from the metagenomic data that the basic pattern of the N-cycle was that ammonium assimilation, mineralization, DNRA, and ANRA were higher, while denitrification and nitrite oxidation, nitrogen fixation, and ammonia oxidation were lower, with anammox being the lowest. This is consistent with the general pattern that ammonium assimilation is common, and nitrogen fixation and ammonia oxidation are less common in most soils (Varin et al., 2010; Nelson et al., 2015, 2016; Souza et al., 2015; Tu et al., 2017), but it is notably different from that in sandy saline soils (Ren et al., 2018), marine oxygen minimum zones (Ganesh et al., 2014), and sediment (McTigue et al., 2016; Rasigraf et al., 2017) with abundant denitrification. The organic-rich biocrusts dwell on the soil surface, in which the material cycle is mainly driven by light energy. In organic-poor deep seas (Rasigraf et al., 2017) or stratified lake bottoms (Llorens-Mares et al., 2015), as well as water environments rich in nitrate (Lüke et al., 2016; Li et al., 2018; Keren et al., 2020), the material cycle is mainly driven by chemical energy. We speculate that the difference in material cycling pattern is determined by the most abundant electron donors and receptors available in the different habitats.

From the perspective of ecological function, the nitrogen cycle in biocrusts is also as expected. In the DON-rich biocrusts, nitrogen fixation would not be abundant because of the high energy investment, polymetallic cofactor, and strict anaerobic conditions. In line with nitrogen fixation, nitrification is not expected to be high, especially NH₄⁺-N, as the inorganic nitrogen form is preferred by dominant organisms Cyanobacteria, lichens, and mosses (Lin and Stewart, 1997; Dahlman et al., 2004; Muro-Pastor et al., 2005; Liu et al., 2013) and is volatile under high temperature and strong radiation. Except for the lack of NH₄⁺-N substrate, hypoxia (Johnson et al., 2005; Ouyang and Hu, 2017) and other factors also inhibit nitrification (Castillo-Monroy et al., 2010). Denitrification with medium abundance is characterized by multiples steps, substrates, and products (N₂/N₂O/NO). However, biocrusts are rich in organic matter and have changeable redox states, which would increase the probability of denitrification (Tiedje et al., 1982). Nitrite oxidation is limited by low NO₂^– in biocrusts, but it is a necessary process to convert inorganic nitrogen to the highest oxidation state and balance the redox potential, which helps explain why nitrite oxidation is in moderate abundance in biocrusts. Finally, ammonium assimilation, mineralization, and nitrate ammonification (DNRA and ANRA) are abundant in biocrusts because all of them play important roles in biocrust development. Among them, ammonium assimilation is an indispensable potential for the growth of most organisms. Mineralization is not only an important source of nutrients for organismal growth but also a necessary process during species replacement (Reverey et al., 2016). ANRA has the function of supplying NH₄⁺-N at low cost, and DNRA has the capacity of energy yield (Cole and Brown, 1980). Both ANRA and DNRA use a two-step reaction to transform nitrogen from the highest oxidation state to the lowest reduction state and increase the nitrogen retention in a given ecosystem (Tiedje, 1988). Roughly, ANRA and DNRA match the abundances of ammonium assimilation and mineralization and eventually reduce the nitrogen loss in biocrusts, particularly when they frequently experience dry-wet cycles in the field. Therefore, the patterns of N-cycle in biocrusts can be regarded as an ideal ecological strategy with low energy cost and high nitrogen retention. RDA analysis (Figures 7C,D) suggest that the nitrogen cycling pattern of biocrusts is mainly driven by DON and NO₃^–. The DON is a highly reducing electron donor, and NO₃^– is a highly oxidizing electron acceptor. Their influence is a multidimensional control of redox potential from high to low or vice versa.

Among the correlation of pathways (Figure 7A) in biocrusts, nitrite oxidation and DNRN occupy a hub position in the network with highest degree, and key genes have close genetic relationships (Kirstein and Bock, 1993; Simon and Klotz, 2013). The two processes are related to the substrates competing via denitrification, DNRN and ANRN, and are mainly affected by C:(N:P) and DOC (Figure 7), which supports the proposal of Stein and Klotz (2016) to define NO₂^–-NO₃^– interconversion as an independent nitrogen-transformation flow. Compared with DNRN, ANRN had a much higher relative abundance (6x). Coupled with similar microbial communities between NO₃^– transportation, assimilation, and ammonium assimilation (Figure 5 and Supplementary Figure S3), our results suggest DNRN is the true hub of nitrogen-transformation network. ANiRA and mineralization (Figure 7) are mainly affected by the characteristic variables (NH₄⁺, Eh, NO₃^–) during hydrological changes (Reverey et al., 2016), while the nutrition strategies are significantly different among the potential bacteria (Cyanobacteria, Proteobacteria, Chloroflexi, Planctomycetes, and Ascomycota). Furthermore, the replacement of ecosystems is often accompanied by both decomposition and synthesis of organic matters caused by mineralization (Reverey et al., 2016), all of which indicate that these two processes closely respond to frequent dry-wet changes. Finally, DNiRA is a more advantageous process than denitrification under high C:NO₃^– (Tiedje et al., 1982; Rütting et al., 2011). The nmo gene (encoding nitronate monooxygenase) in biocrust is comparable to nirB which may be employed for supplying NO₂^– (Li et al., 2018). DNiRA and nitrogen fixation are strongly negatively correlated, and they are mainly affected by the opposite effects of NH₄⁺, C:(N:P), DON and NO₃^–. NH₄⁺ and C:(N:P) are subjected to stoichiometric constraints, especially for NH₄⁺-N. And, DON and NO₃^– are regulated by redox state. Therefore, nitrite oxidation and DNRN are the hubs connecting N-cycle pathways, DNiRA and nitrogen fixation are the most important ways to maintain the balance of elemental stoichiometry and redox states, while ANiRA and mineralization are processes that closely respond to dry-wet cycles.

Environmental variables such as Eh, pH, NH₄⁺, NO₃^–, C:(N:P) influence N-cycle pathways (Pett-Ridge et al., 2006; Burgin and Hamilton, 2007; Rütting et al., 2011; Sgouridis et al., 2011; Bouskill et al., 2012; Jiang et al., 2015; Zheng et al., 2020). In this study, we found that pH played the opposite role in the N-cycle compared to many other variables (C:(N:P), NH₄⁺, DON, SO₄^2–, C:N, WC, DOC). The effects of pH and C:(N:P) were the most significant, with pH negatively correlated to nitrogen fixation and positively correlated to denitrification. Different diazotrophic taxa respond differentially to soil pH (Wang et al., 2017), and free-living N fixers prefer neutral pH (Limmer and Drake, 1996), thus leading pH negatively correlate with N fixation in the slightly alkaline biocrust (pH 7.43–8.09), which is consistent with previous studies (Keshri et al., 2015; Wang et al., 2018a). The significant effects of C:(N:P) and C:N on nitrogen fixation were similar to those in forest succession, although single nutritional factors such as DON and SO₄^2– did not affect nitrogen fixation significantly in forest ecosystems (Zheng et al., 2020). The increase in soil pH would solubilize organic matter and increase denitrification potential (Anderson et al., 2018). The negative correlation of SO₄^2– and denitrification have been found previously (Wang et al., 2018b), for the reduction of SO₄^2– can enhance DNRA and inhibit denitrification (Tiedje, 1988). More interestingly, the negative correlation of C:(N:P) with DNiRA, nitrite oxidation, DNRN, and ANRN decreased sequentially, and the positive correlation of pH with the above four pathways also decreased sequentially. This order is consistent with the order of thermodynamics-driven redox reactions (Falkowski et al., 2008; Ettwig et al., 2010), and, thus, we speculate that these four are the central pathways of N-cycle in biocrusts. Many variables significantly affect multiple pathways, but the key variables affecting the main pathways are different, in which pH and SO₄^2– are the strongest variables affecting denitrification, and C:(N:P) is the strongest variable affecting nitrogen fixation, DNiRA, nitrite oxidation, and DNRN.

In order to further improve the reliability of metagenomic sequencing results (Figure 2 and Supplementary Figure S4) and compare gene abundances between different successional stages, we adopted the same result as GeoChip (Figure 3), qPCR, or potential rate (Figure 4). It has been recognized that acetylene reduction assay has limitation to examine potential N fixation rate such as differences in solubility of acetylene and N₂ in water, adsorption of acetylene on soil colloids and toxicity to microbes, and ¹⁵N stable isotopes should be used instead (Barger et al., 2016). In this way, with the succession of natural biocrusts, the increase of nrfA and nosZ was associated with the increase of Deltaproteobacteria, Verrucomicrobia, Gemmatimonadetes (Figure 5); the decrease of nirA and ureC was related to the decrease of Cyanobacteria and the increase of Bacteroidetes, Alphaproteobacteria and Actinobacteria (Figure 5). However, with the development of simulated biocrusts, the increase of nirA and ureC was consistent with an increase in Cyanobacteria; the decrease of nrfA was associated with a decrease in weakly resistant bacteria such as Deltaproteobacteria and Acidobacteria. The changes of nrfA and nirA may be related to the abundance of electron donors and receptors under the influence of DON and NO₃^– (Figure 7). Therefore, along with the succession of biocrusts, the changes of N-cycle follow the direction that ANiRA and mineralization gradually decreasing, respiratory DNiRA increasing, and more thorough denitrification, which indicated the improvement of nitrogen balance ability.

Our co-occurence network analysis (Figure 6) showed that the complexity and interaction degree between microbes and N-cycle pathways were significantly increased with succession of biocrusts, which supports the proposal of enhancement of microbial interactions that promote ecological succession (Datta et al., 2016) and fully displays ecosystem stability and functional enhancement of the N-cycle (Gravel et al., 2016). From the perspective of coupled relationships between species and genes, there were several obvious functionally independent groups, including two large groups in the early stage (A) coupled with denitrification, or ammonium assimilation (glnA) with higher affinity, and a group coupled with nrfA in the later stage (M), which is consistent with the harsher conditions of early stage (Geisseler et al., 2010; McTigue et al., 2016) and more nitrogen retention with increased organic matter in later stage (Moreno-Vivian et al., 1999; Rütting et al., 2011). The assimilation nitrate reduction group in the transitional stage (C, G) suggests that competition for nitrogen resources between fungi and bacteria exists in lichen crusts, which is likely accomplished by the process of nitrate reduction. However, nifH was clustered with nrfA in algae crusts and cyanolichen crusts, and nirB in chlorolichen crusts and moss crusts (Figure 6). Our data indicate that the functional groups and transformation processes of regulating nitrogen fixation are mainly the coordinated effects of Cyanobacteria and respiratory DNiRA function in the early (A) and early-middle stages (C), but the result of Actinobacteria and fermentive DNiRA in middle-late (G) and later stages (M).

As biocrusts in our study have adapted to dry-wet cycles for a long time, their functional microbial communities should be the result of long-term adaptation. In early natural biocrusts, nitrogen-fixing organisms were mainly dominated by Cyanobacteria (Yeager et al., 2012; Wang et al., 2016), while Proteobacteria played an important role in other biocrusts where nitrogen-fixing Cyanobacteria were not abundant (Pepe-Ranney et al., 2016). For aerobic ammonia oxidation, except for the reported Thaumarchaeota and Betaproteobacteria (Marusenko et al., 2013; Stein and Klotz, 2016), the proportion of unknown archaea in our study was also high. Microbes contributing to NO₂^– and NO₃^– interconversion were mainly Proteobacteria (Alpha-, Beta-, Delta-) and Actinobacteria (Alcantara-Hernandez et al., 2009; Daims et al., 2016), and Nitrospirae contributed more to simulated biocrusts. Compared with the microbes involved in DNiRA (Seitz and Cypionka, 1986; Guo et al., 2016; Cerqueira et al., 2018; Li et al., 2018; Keren et al., 2020) or ANRA (Nelson et al., 2015; Li et al., 2018; Ren et al., 2018) in other environments, there were more Actinobacteria, Cyanobacteria, and Verrucomicrobia in biocrusts associated with DNiRA, and more Cyanobacteria, Planctomycetes, and Ascomycota associated with ANRA. In denitrification, in addition to the reported Proteobacteria, Bacteroidetes, and Chloroflexi (Saarenheimo et al., 2015; Guo et al., 2016; Ren et al., 2018), there were abundant Actinobacteria (nirK, norB), Cyanobacteria (norB), and Gemmatimonadetes (nosZ) involved in the biocrusts. Overall, Alphaproteobacteria, Actinobacteria, and Cyanobacteria have the highest nitrogen transformation potential, while microbes that contained genes such as amoA, nosZ, and nirK were a higher proportion of unknown taxa.

From the perspective of physiology and ecology, the cyanolichen crust (C) being dominated by Collema in this study is a classic adaptive balance strategy for rapid biosynthesis under high temperature after sufficient rainfall (Lange et al., 1998; Grote et al., 2010). However, a deficit in material accumulation will occur when only a small amount of water is supplied (Ouyang et al., 2017). The other three types of biocrusts displayed common adaptations to water pulse (Lan et al., 2017; Ouyang et al., 2017). From our metagenomic analyses, the abundances of nifDHK, nrfAH, amoABC, and norBC were extremely high in the cyanolichen crusts, and gene copies of nifH, nrfA, and potential rate of nitrogen fixation and ammonia oxidation were also the highest. The only lower production of N₂O in the cyanolichen crusts was related to the low DNRN and low concentration of NO₃^– and NO₂^– (Supplementary Table S4). This suggests that the difference in nitrogen transformation rate is only comparable when the samples have similar redox background, and the biocrusts with special survival strategies are not suitable for the comparison, which may be the reason for the divergences in the results from previous studies on potential nitrogen-transformation rates.

Conclusion

This study conducted repeated sampling and metagenomic sequencing throughout 4 years and points toward the overall patterns and successional variation of N-cycle, the microbial contribution, the driving mechanism of environmental variables, the relevance between N-transformation pathways and between N-cycle and ecosystem functions. Our results highlight the consistent N-cycle pattern of low energy cost and high nitrogen retention ability in biocrusts, and an obvious successional variability of decreased responsive ANiRA and mineralization, increased respiratory DNiRA, and more thorough denitrification. Our results demonstrate that N-cycle in biocrusts is driven by DON and NO₃^–, which indicates a multidimensional control of redox potential. The linkage of N-cycle processes and the relevance with environmental variables show that nitrogen fixation and DNiRA are the most important pathways to maintain the balance of elemental stoichiometry and redox states, while ANiRA and mineralization are processes responding closely to the wet-dry cycles. In addition, our study suggests that the coordination of DNiRA and nitrogen fixation is crucial for N-cycle, and the role of two types of DNiRA changes with succession. This study indicated that biocrusts were an ideal natural model to research the function and succession of microbial ecosystems with metagenomic methods.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The name of the repository and accession numbers can be found below: National Center for Biotechnology (NCBI), https://www.ncbi.nlm.nih.gov/, PRJNA562239, PRJNA548650, PRJNA646196, PRJNA647185, PRJNA544139, PRJNA321565, PRJNA647637, PRJNA647192, PRJNA647699, PRJNA543295, PRJNA636669, and PRJNA636680.

Author Contributions

CH designed the research and aided in data analysis. QW performed major experiments, carried out the analysis, prepared most figures. YH carried out parts of experiments of 2017 simulated biocrusts. CH and QW wrote this manuscript. SL provided advice and critically reviewed the manuscript. All authors read and approved the final manuscript.

Funding

This work was financially supported by the Strategic Priority Research Program of Chinese Academy of Sciences (XDA17010502), the National Natural Science Foundation of China (41877339), and Forestry Research Special Funds for Public Welfare (201404204).

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.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmicb.2021.633428/full#supplementary-material

Supplementary Figure 1 | Schematic of the experimental and analytical workflow. A, C, G, M represent four types of natural biocrusts at different succession stages. A: algae crusts, C: cyanolichen crusts, G: chlorolichen crusts, M: moss crusts. Ae, Am, and Al represent early, middle, and later stage of simulated cyanobacterial crusts, respectively.

Supplementary Figure 2 | Phylogeny of nifH (A) and nirK (B) genes from natural biocrusts in 2015. The heatmap and percentage from left to right are responding to the four types of natural biocrusts: algae crusts, cyanolichen crusts, chlorolichen crusts, and moss crusts, respectively. The evolutionary history was inferred using the Neighbor-Joining method in MEGA7.

Supplementary Figure 3 | The potential microbial community involved in the transport of nitrate and nitrite in natural biocrusts in 2015. NRT, corresponding to K02575, encoding nitrate/nitrite transporter; nrtA, corresponding to K15576, encoding nitrate/nitrite transport system substrate-binding protein. A: algae crusts, C: cyanolichen crusts, G: chlorolichen crusts, M: moss crusts.

Supplementary Figure 4 | The relative abundance of N-cycle pathways (A) or genes (B) based on metagenomic analysis in the natural biocrusts (5 batches) and simulated biocrusts.

Supplementary Table 1 | The qPCR primers and conditions used in this paper.

Supplementary Table 2 | Assembly and annotation statistics of metagenome and 16S rRNA gene sequences (ORF, open reading frame; QC, quality control).

Supplementary Table 3 | The corresponding relationship of N-cycle pathways, genes and KEGG Orthology (KO) number.

Supplementary Table 4 | The properties of biocrusts. Data are reported as average ± SD. Values with different letters are significantly different at 0.05 level within natural biocrusts (n = 3) or simulated biocrusts (n = 4). WC: water content; EC: electrical conductivity; Eh: redox potential; TC: total carbon; TOC: total organic carbon; BC: biomass carbon; RC: remnant carbon; TN: total nitrogen; TP: total phosphorus; DON: dissolved organic nitrogen. ND: not detected; NC: not collected.

Footnotes

References

Abed, R. M., Lam, P., de Beer, D., and Stief, P. (2013). High rates of denitrification and nitrous oxide emission in arid biological soil crusts from the Sultanate of Oman. ISME J. 7, 1862–1875. doi: 10.1038/ismej.2013.55

PubMed Abstract | CrossRef Full Text | Google Scholar

Alcantara-Hernandez, R. J., Valenzuela-Encinas, C., Marsch, R., and Dendooven, L. (2009). Respiratory and dissimilatory nitrate-reducing communities from an extreme saline alkaline soil of the former lake Texcoco (Mexico). Extremophiles 13, 169–178. doi: 10.1007/s00792-008-0207-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Anderson, C. R., Peterson, M. E., Frampton, R. A., Bulman, S. R., Keenan, S., and Curtin, D. (2018). Rapid increases in soil pH solubilise organic matter, dramatically increase denitrification potential and strongly stimulate microorganisms from the Firmicutes phylum. PeerJ 6:e6090. doi: 10.7717/peerj.6090

PubMed Abstract | CrossRef Full Text | Google Scholar

Austin, A. T., Yahdjian, L., Stark, J. M., Belnap, J., Porporato, A., Norton, U., et al. (2004). Water pulses and biogeochemical cycles in arid and semiarid ecosystems. Oecologia 141, 221–235. doi: 10.1007/s00442-004-1519-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Barger, N. N., Weber, B., Garcia-Pichel, F., Zaady, E., and Belnap, J. (2016). “Patterns and Controls on Nitrogen Cycling of Biological Soil Crusts,” in Biological Soil Crusts: An Organizing Principle in Drylands, eds B. Weber, B. Büdel, and J. Belnap (Cham: Springer International Publishing), 257–285. doi: 10.1007/978-3-319-30214-0_14