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

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 NO3–. pH and SO42– 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 2 − -NO 3 − 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 4 + , NO 3 − , 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;. 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 2 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.

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) . 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 postinoculation 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)  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 3 2− and HCO 3 − were measured by double indicator titrations, that of NO 3 − , NO 2 − , NH 4 + , PO 4 3− , Ca 2+ , and Mg 2+ 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 1 .

Metagenomic Sequencing and Assembly
Metagenomic libraries were sequenced on Illumina platforms (Supplementary Table S2). Raw sequence data quality was assessed using FastQC 2 , trimmed using Seqprep 3 , 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 4 . 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 5 . Contigs (≥300 bp) were analyzed for open reading frames (ORFs) prediction using MetaGene (Noguchi et al., 2006). Nonredundant 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 4 ) 2 SO 4 and 75 mM NaClO 3 and incubated on a shaker at 24 • C for 10 h. The change in NO 2 − 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 3 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 2 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.

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.
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.
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).
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.

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.
For the environmental factors that influenced each N-cycle pathway, the correlation network ( Figure 7B) showed that, NH 4 + 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 3 − , and HCO 3 − were negatively correlated with both ANiRA and mineralization, while SO 4 2− was most strongly negatively correlated with denitrification. In short, both pH and SO 4 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 3 − 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 3 − 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., 2015Nelson et al., , 2016Souza 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 4 + -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 4 + -N substrate, hypoxia (Johnson et al., 2005; and other factors also inhibit nitrification (Castillo-Monroy et al., 2010). Denitrification with medium abundance is characterized by multiples steps, substrates, and products (N 2 /N 2 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 2 − 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 4 + -N at 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).
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 3 − . The DON is a highly reducing electron donor, and NO 3 − 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 2 − -NO 3 − 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 3 − transportation, assimilation, and ammonium assimilation ANiRA and mineralization (Figure 7) are mainly affected by the characteristic variables (NH 4 + , Eh, NO 3 − ) 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 3 − (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 2 − (Li et al., 2018). DNiRA and nitrogen fixation are strongly negatively correlated, and they are mainly affected by the opposite effects of NH 4 + , C:(N:P), DON and NO 3 − . NH 4 + and C:(N:P) are subjected to stoichiometric constraints, especially for NH 4 + -N. And, DON and NO 3 − 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 4 + , NO 3 − , 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 4 + , DON, SO 4 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 4 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 4 2− and denitrification have been found previously (Wang et al., 2018b), for the reduction of SO 4 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 4 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 2 in water, adsorption of acetylene on soil colloids and toxicity to microbes, and 15 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 3 − (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).
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 . The other three types of biocrusts displayed common adaptations to water pulse . 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 2 O in the cyanolichen crusts was related to the low DNRN and low concentration of NO 3 − and NO 2 − (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 3 − , 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.

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.