Multiple role of diazotrophs in coupling nitrogen and iron biogeochemical processes through vertical soil proles of different paddy soil types

Background: Soil microbiota exert fundamental functions in maintaining ecosystem functioning and services, including pedogenesis, biogeochemical processes and plant productivity, especially for agriculture system. Despite their ubiquitousness from the epipedon to deep soil, the vertical characteristics of microbiomes (especially for functional microorganisms) and their contribution to soil element cycling when considering soil developmental features are poorly understood. Here, nine proles (0~135 cm) of two canonical paddy soil types (Fe-accumuli- and Hapli-stagnic anthrosols; 111 samples in total) at a local scale were collected, which represented relative long- and short-term water ooding history, respectively. The vertical variations in edaphic characteristics and assemblies of soil bacterial and diazotrophic communities, and microbial contribution to element cycling were explored. Results: Across soil proles, Hapli-stagnic anthrosol was characteristic of higher concentrations in free iron oxides and total iron in the epipedon, and contained higher amounts of ammonia along the subsurface layers, as compared with acidic Fe-accumuli-anthrosol. Community assemblies of bacteria and diazotrophs, as well as edaphic properties, were mainly shaped by soil depths, followed by soil types. Furthermore, random forest analysis revealed that, for Fe-accumuli-anthrosol, available Fe could best predict nitrogen cycling index and nitrogen status was signicantly related to iron cycling index; while in Hapli-stagnic anthrosol, available sulfur was the most important variable in predicting nitrogen and iron cycling indices. Among the dominant genera, some distinctive biomarkers that varied remarkably between the two soil types were noticeable for their contributions to both nitrogen and iron transformation, including iron-reducing diazotroph Geobacter and iron-oxidizing bacterium Rhodanobacter that characterized Fe-accumuli type, and sulfur reducing diazotroph Desulfobacca as main discriminant clades for Hapli-stagnic type. Conclusions: A novel perspective was proposed on the vertical characteristics of edaphic properties and bacterial and diazotrophic communities in the two paddy soil types. The ndings indicated the nitrogen-iron cycling processes for Fe-accumuli-anthrosol and nitrogen-iron-sulfur coupling interaction for Hapli-stagnic anthrosol, advancing our understanding of the signicant multiple role played by soil microorganisms, especially for diazotrophs, in element biogeochemical cycles.


Introduction
Pedosphere plays indispensable role in supporting functions and services of above-and belowground communities in terrestrial ecosystem [1] . Along the pedogenesis, soil gradually differentiate into different phenotypes, depended on soil-forming factors, including climate, parent material, organisms, relief and time [2] . The characteristics and functionality of varied phenotypes are highly variable [3][4][5][6] . Soil types impose strong selection pressures for microbiome with traits re ecting different resource needs and tolerances to environmental conditions [6,7] . In turn, microbiota serve as key drivers and indicators for soil multiple ecosystem functions ('multifunctionality'), including element biogeochemistry cycles, primary production and natural attenuation of pollutants [1,8,9] . However, numerous studies on successional trajectories and functionality for microbial community in different soil types mainly focused on epipedon [6-8, 10, 11] , overlooking the majority of microbiota in the subsurface ecosystems, especially for functional microorganisms [9] .
Soil-forming processes frequently leave markers in the forms of diagnostic horizons and features that are discontinuous with depth [12] . As the representative for pedogenesis, a soil pro le represent the vertical section of the soil that extends from the soil surface to the parent rock material [5,12] . Researchers showed that the subsurface horizons compass 35 ~ 50% of the total microbial biomass, and the microbial activity and interactions are as intense in the subsoil as those in the topsoil [13] . Divergent assembly pattern of microbial taxonomic and functional compositions have been observed between the topsoil and deeper layers, which directly coupled with soil biogeochemical processes [9,14,15] . Previous studies showed that soil horizons exerted overwhelming in uence on bacterial and fungal community structures than sampling locations [13,16,17] , depending on geographical distance, climate and edaphic properties [3,6] . Bai et al. (2017) reported the strong effects of soil parent material characteristics over soil horizons [18] . However, most of these studies mainly focused on the super cial subsurface layers (e.g. 15 30 cm) [17] or arti cially divided the soil horizons into equal layer (e.g. for each 20 cm) [14,18] , neglecting the edaphic developmental characteristics along pedogenesis and weakening the differences among horizons. Particularly, the assemblage pattern of functional microorganisms community and their functioning along soil pro les are rarely de ned. N 2 xing microorganisms provide the predominant source of active N to the biosphere (e.g. ~50-70 Tg N year − 1 to agricultural systems) [19,20] . As the key driver for global nitrogen (N) dynamics, diazotrophs also participate other major biogeochemical cycles directly or indirectly, such as for carbon sequestration, iron (Fe) and sulfur (S) cycling [19,[21][22][23] . The abundance, composition and N 2 -xing capacity of diazotrophs varied greatly across soil types, while only studied in the plough layer (0-15 cm), which processes are determined by soil properties, climate and geographical distance [3,24] . Previous studies presented that diazotrophic population was more abundant in surface soils due to larger organic matter contents and available N limitation [25,26] . However, the diversity and structure of diazotrophic community from deep soil horizons and their contribution to element biogeochemical cycles with depth remains unknown.
With more than 50% of world`s population depending on rice for subsistence, it is vital to study soil nutrient status, functionality and constraining factors that control rice productivity [27] . As a key limiting nutrient for rice growth, bioavailable N input in paddy soils is heavily depending on diazotrophic community [24] , which also participated in other biogeochemical cycling (e.g. Fe and S) [28][29][30][31][32] . Across the pro le of paddy soils, the alternation between waterlogging and drying results in a vertical oxygen gradient and periodically occurring redox reactions, including intense activity in N, Fe and S cycles. These unique characteristics and intense element biogeochemical processes renders paddy soil as one model system between terrestrial ecosystems and aquatic ecosystems. Owing to distinct ooding history, different paddy soil types with varied redox dynamics are formed that are supposed to harbor divergent microbiota assembly and biogeochemical processes along soil pro les, whereas the microbial mechanism remains obscure.
The objective of this study was to determine bacterial and diazotrophic community assemblies, as well as edaphic properties, across the pro le of different paddy soil types, and the function of these microbial taxa. Our sampling zone, one major rice production area in the Yangtze River basin (China) [33] , offers proper sites at a local scale with the attempt to minimize the effects of spatial distance and climate ( Fig. 1A). Nine pro les from two typical paddy soil types were used: 1) Fe-accumuli-anthrosol with relative longer waterlogging history is characteristic of redoximorphic status, whereas 2) Hapli-stagnic anthrosol with short ooding management is an oxidative type [12,34,35] . According to the redoximorphic features and degree of paddy soil development, the soil pro les were divided into 3-5 horizons. The abundance and structure of bacterial and diazotrophic community were assessed using quantitative polymerase chain reaction (qPCR) and Illumina MiSeq sequencing. We posited that 1) both paddy soil type and depth matters for edaphic characteristics and microbial community assembly, 2) the variation in microbial community assembly across soil pro les would generate differences in soil functioning between distinct paddy soil types.

Results
Vertical pro les of edaphic properties in different soil types The vertical variations in morphological characters and edaphic properties with depth were observed across the pro les of paddy soils ( Figure 1B, S1). It was noted that pH elevated and approached 7 with the increase in soil depths in all samples (P < 0.001). When compared with topsoil, the deeper layers were depleted in partial elements, including OM, N (total N, available N and NH 4 + -N), P (total and available P), Fe (total and available Fe), Mn (total and available Mn) and available Cu/Zn/Mo (P < 0.01), while K, Ca and total Mg/Cu/Zn were insensitive to soil depths.
With regard to the two paddy soil types, pH, Fe and N were the most differentiated edaphic properties ( Figure 1C, S2). In Fe-accumuli stagnic anthrosols, the upper layers (w1~w2) appeared to be more acidic with lower pH (4.90~6.61) when compared with those in Hapli-type (pH 5.47-7.27). Fe-accumuli type that possessed typical Fe-Mn cutans in the hydragric horizon (Fe-accumulate diagnostic horizon), was characteristic of more than 1.5 times higher levels of free Fe oxides in the deep layers than that in epipedon ( Figure S2). The total Fe and Fe oxides in Fe-accumuli type were also signi cantly lower than that in surface layers of Hapli-type (P < 0.0001), while possessing higher available Fe content ( Figure S2).
The lower concentrations of total N and NH 4 + -N were also exhibited in the deep layers (w3 and w2-w4, respectively) of Fe-accumuli soils than those of Hapli-stagnic soils (P < 0.05 and P < 0.01, respectively). Besides, Fe-accumuli soils tended to be rich in total K (w2-w4), total Zn (w3), available Zn (w2) and total Mo (w4), whereas Hapli-stagnic soils possessed higher nutrient levels with available Ca/Mg and total Mg/Cu/Zn in the super cial layers (w1) and available Mo in the deepest layers (w4) (P < 0.05). The differentiations were also revealed by Random Forest (RF) classi cation analysis, which showed that soil type could be best distinguished by free Fe oxides for the topsoil, while NH 4 + -N and Mo (total and available Mo) was of the most importance in the deep layers ( Figure S3).

Microbial community assembly shifted between soil types
The assembly patterns of total bacterial and diazotrophic communities between the two soil types were analyzed. Along soil pro les, bacterial and diazotrophic α-diversity indices, including observed operational taxonomic unit (Sobs; richness) and Shannon (taxonomic diversity) displayed a decrease trend as soil depth increased, which tendency was also observed in the gene abundances of the total bacteria and nifH genes quanti ed by qPCR ( Figure S4). However, there was no signi cant difference observed between the microbial diversities and abundances of the two soil types, except for signi cant higher values for bacterial and diazotrophic abundances in w2 layer of Fe-accumuli type.
The turnovers of microbial community compositions were visualized with PLS-DA analyses. Samples formed distinct clusters according to soil depths and types ( Figure 2A-B), with statistically signi cant differences being noted at taxonomic levels (P ≤ 0.001, ANOSIM/ADONIS) (Table S2). To determine whether deterministic or stochastic processes best explain the assembly of the bacterial and diazotrophic communities, mean nearest taxon index (NTI; indicates the 'terminal' phylogenetic dispersion of the community) and βNTI were calculated. When considering soil depth, mean NTI scores were all positive, with significantly higher values in the topsoil (Table S3). βNTI values tended to approach the range of -2 to +2 with the increase in soil depths (Table S4). The data suggested that these microbial communities from the epipedon were more phylogenetically clustered than that of the deeper layers. By contrast, the variance of mean NTI values between the two soil types was slight (Table S3). The majority of βNTI values were above +2 (59.76~64.58%) for bacterial assembly, while the scores for diazotrophic community were mainly in the range of -2 to +2 (54.20~58.25%), with Hapli-stagnic soil always exhibiting higher proportions ( Figure S5). These data revealed that the deterministic processes dominated the phylogenetic community dynamics of bacterial assembly, while stochastic processes for diazotrophic community turnover in these paddy soil samples.

Environmental factors shape microbial assembly
As revealed by partial Mantel test between soil types in the whole layers, diazotrophic community was more sensitive to geographic distances (R = 0.28~0.32) with higher regression slopes (R = 0.13~0.19) than environmental heterogeneity did, while spatial and environmental distances contributed similarly to bacterial community dissimilarities (Table S5). The data were further con rmed by distance-decay pattern ( Figure S6). For Fe-accumuli soils, both environmental and geographic variables exerted more impacts on diazotrophic community dissimilarities (R = 0.19 and 0.32, respectively) than those in Haplistagnic group (R = 0.13 and 0.28, respectively). The similar trend was also noted for bacterial community.
When considering edaphic factors only, α-diversities and abundances of bacteria and diazotrophs were positively correlated with OM, total N and available nutrients (N, Fe, Cu, Zn and Mo), while negatively linked with pH ( Figure S7). Mantel test further showed that the turnover in bacterial and diazotrophic community compositions exhibited signi cant and strong correlations with most edaphic properties (e.g. soil pH and available Mg/Ca/N/Fe), except for soil texture (clay, sand, silt) and total Mg/Mo (Table S6).
When considering the two soil types separately, Fe-type soil possessed more signi cant correlations between microbial community assembly and edaphic factors as compared with Hapli-type soil did ( Figure S8). Among the three most differentiated indices (NH 4 + -N, available Fe and free Fe oxides) between the two soil types, their signi cant correlation with bacterial and diazotrophic community structures was observed in Fe-accumuli type, while the correlation was weak in Hapli-type soil, except for the effect of available Fe on diazotrophic community.

Impacts of soil depth and type on edaphic properties and microbial community assembly
To evaluate the effects of soil depth and type on edaphic properties, two-way ANOVAs analysis was applied ( Figure 2E). pH, available N/Fe, NH 4 + -N and total P/Mn were signi cantly in uenced by both soil depth and type (P < 0.05). OM, total N/Fe and available S/P/Mn/Zn/Cu were solely signi cantly in uenced by soil depth (P < 0.05), while soil types exerted impacts speci cally on K contents and partial micronutrients (total Zn/Mo and available Ca/Mg) (P < 0.05).
Microbial populations were also signi cantly sensitive to both soil depth and type. Combined with PERMANOVA, soil depth was the major determinant for rendering bacterial and diazotrophic community structure separation, which explains 13.04% and 9.42% of the observed structure variation (P = 0.001), respectively, followed by soil types (7.58% and 6.72%, P = 0.001) ( Figure 2F).

Accounting for potential drivers of soil N and Fe cycling processes
To comprehensive evaluate soil biogeochemical functions, the differentiated N and Fe cycling indices were calculated by normalizing and standardizing each of the N/Fe-related nutrient properties ( Figure 3A) [36] . The main variables for predicting the variation in soil N and Fe cycles of the two soil types were identi ed by RF analysis. Throughout the whole pro les, the models were highly predictive by edaphic properties according to the determination values (R 2 = 0.94~0.99, P = 0.01), and the values were 69%~80% for microbial variables ( Figure 3B-3C). In Fe-accumuli anthrosol, available Fe was the most important variable for explaining the variation in N cycling index, and Cu was the most important explanatory factor for Fe cycling index. It was notable that varied status of N contents (NO 3 -N, NH 4 + -N and available N) also contributed signi cantly to Fe cycling index in Fe-accumuli type (P < 0.05). By contrast, available S acted as the most important variable for predicting both N and Fe cycling indices in Hapli-type soil ( Figure 3B). When considering microbial predictors alone, diazotrophic and bacterial abundances best predicted the dynamics of soil N and Fe cycling indices, followed by microbial βdiversity ( Figure 3C).
When considering soil horizons, the contribution of dominant microbial genera to soil N and Fe cycling indices varied remarkably among each layer ( Figure 4). Some bacterial predictors were important for both N and Fe cycling, such as for Rhodanobacter that was characteristic of Fe-accumuli type soil. The genus belonging to family HSB_OF53_F07, Thiobacillus and Oryzihumus were pivotal in predicting N cycling in topsoil (P < 0.01), Paenibacillus, unclassi ed_order_Myxococcales and norank_phylum_SBR1093 in the secondary layers, and norank_order_JG30_KF_AS9, norank_class_S085 and norank_phylum_SBR1093 in the deep layers ( Figure 4A). For Fe cycling index, bacterial unclassi ed_family_Intrasporangiaceae and Gaiellanorank_order_JG30_KF_AS9 best explained the variations in epipedon, subsurface and bottom layers, respectively ( Figure 4B).
With respect to diazotrophic genera, Dechloromonas and Paludibacter contributed signi cantly to soil N cycling index, while Derxia/Methylocystis and Nitrospirillum/Burkholderia were among the best predictable variables for secondary and deep layers, respectively ( Figure 4C). Notably, some diazotrophs were also vital in predicting soil Fe cycling process, including Leptothrix, Derxia, Anaeromyxobacter and SRB (Desulfobacca, Desulfotomaculum and Desulfatibacillum, which were characteristic in Hapli-type soil) ( Figure 4D). Geobacter that was enriched in Fe-accumuli type was important for predicting for N cycling index (P < 0.05) in the bottom horizons ( Figure 4B-C), and the signi cant correlation between Geobacter abundance and available Fe concentration in Hapli-type soils was further noted by linear regression analysis ( Figure S13).

Discussion
Given the importance of maintaining a productive soil environment for sustainable crop production, improving the understanding of the driving mechanism for soil quality and functioning is highly bene cial [8] . With clues provided by soil pro le, we can draw a detailed picture of what a soil is and what a soil does from the point of soil development [5] . Here, we found that soil nutrient status (especially for higher N and Fe cycling indices in Hapli-type soil) and microbial community assemblies varied signi cantly and vertically in paddy soils at a local scale, depending on soil horizons, followed by paddy soil types. For Fe-accumuli-anthrosol with long-term waterlogging history, available Fe could best predict N cycling index, and meanwhile varied forms of N were signi cantly related to Fe cycling index, where iron-reducing diazotroph Geobacter and iron-oxidizing bacterium Rhodanobacter might be involved in these processes. By contrast, available S was the most important variable in determining N and Fe cycling indices in the short-term waterlogging Hapli-stagnic anthrosol, where the main discriminant clade sulfur reducing diazotroph Desulfobacca might act as the driver for the turnover of N and Fe cycling. Our results proposed two distinctive pictures depicting the divergent biogeochemical processes in varied paddy soil types that could be predicted by microbial community compositions (Fig. 5).
First of all, our data revealed that soil pro les exhibited more signi cant impacts on edaphic characteristics and microbial community structure as comparing with soil types. Generally, the epipedon possessed higher nutrient levels than deep layers did, especially for limiting nutrients for plant growth (e.g. OM, total N, NH 4 + -N, total P, available P and available microelements) ( Figure S1). The enrichment might be attributed to high microbial activity and agricultural practices (e.g. fertilization) in the plough layers [9] , which could be evidenced by the higher microbial abundances and diversities in the topsoil ( Figure S4). Microbial communities from the epipedon were more phylogenetically clustered than that of the deeper layers, which might be due to the strong heterogeneous selection and progressive complexity in epipedons [7] . Along soil pro les, the heterogeneous distribution and changes in valence states of elements could be mainly in uenced by vertical oxygen gradient, moisture content, light, pH, parent materials and plant characteristics (e.g. biomass cycling rates and root distributions) [37] , thereby determining the succession of microbial community at the bottom soil [38] . A few researches have reported the decrease in soil diazotrophic abundances with depth [26] . Our data went further, representing the detailed changes in diazotrophic diversity and structure along soil pro les, probably due to the decreased oxygen which is unfavorable for the most microbiomes [9] .
Paddy soil type acted as the secondary explaining variable for the variation in edaphic factors and microbial community assembly. For a certain layer, the abundance and diversity of microbial community were generally consistent between the two paddy soil types, indicating the relative semblable and stable successional trajectories of microbiome under similar soil-forming processes (paddy management). However, the variations in structure and composition of bacterial and diazotrophic communities between the two paddy soil types were signi cant while slight, which could be mainly attributed to the different rice cultivation years and underground water table. Furthermore, our data showed that geographic distance exerted greater impacts on microbial community structure than environmental dissimilarity did ( Figure S6, Table S5). By contrasting to the overwhelming contribution of soil depth over paddy soil groups to microbial community studied here, the great impact of soil types on microbiomes was reported previously at larger scales [3,6,18,24] . Han et al. reported that soil types possessed signi cant divergence in diazotrophic community in the topsoil of farmland, which was signi cantly affected by spatial distance, climate and soil properties [3] . Bai et al. (2017) showed that soil type exerted more signi cant impacts on bacterial and fungal diversity than soil pro le did with distant geographic distances (e.g. more than 1000 km), and the manifest variance in functional structure was only observed among soil types, but not across pro le depths [18] . Combined with previous results, our data highlighted the need for sampling multi-layered soil pro les in order to comprehensively understand the structure and function of microbiomes at a local scale.
The two representative paddy soil types identi ed here, displayed divergent morphology and pedological features (especially for soil pH, Fe and N contents) (Fig. 1B-1C). In Fe-accumuli-anthrosols, the epipedon tended to be bleaching and was characteristic of 1.5 times fewer free Fe oxides than the subsoils, which was also lower than that in Hapli-stagnic topsoil ( Figure S2). The loss of Fe oxides in Fe-accumuli soil might be due to intense ferrolysis (reduction of Fe 3+ to Fe 2+ ) in the topsoil and thereby reductive eluviation of more soluble Fe 2+ under long-term anoxic waterlogging condition [34] , which was indicated by the enriched available Fe in the epipedon. This was also evidenced by the notable deposition of extremely insoluble Fe/Mn oxides in the hydragric horizon of Fe-accumuli type, which imparted a reddish color to the soil [39] . Besides, the acidic condition (pH 5.34 ± 0.49) of Fe-accumuli type might also facilitate Fe activity and availability by Fe reduction [40] . By contrast, weak vertical migration of Fe oxides is one typical feature for Hapli-stagnic anthrosol which undergoes relative short-term rice planting period [12,35] . N cycling index was another differentiated edaphic factor between the two soil types, with a more sharply decreasing trend in NH 4 + -N content of Fe-accumuli soil with depth as compared with Hapli-stagnic soil. In paddy soils, NH 4 + -N xed by diazotrophs serves as the major source of reactive N supply for rice and heterotrophic organisms [20] . The depletion of NH 4 + -N in the deep layers of Fe-accumuli anthrosol might be attributed to the detrimental impact on diazotrophic growth and nitrogenase activity under acidic condition in ferrosols [3] . Previous studies showed that N 2 xing ability in acidic soils was weaker than that under neutral conditions [24] . It is notable that acidic Fe-accumuli soils tended to harbor high abundance of nifH gene ( Figure S4) that was consistent with previous studies, which might represent a compensation strategy for N de ciency utilized by indigenous diazotrophs [24] . Along the chronosequence of rice cultivation, the relative long-term waterlogging management in Fe-accumuli soils would be conducive to anoxic ammonia oxidation, nitrate leaching and denitri cation, thereby resulting in N loss, especially for the anaerobic bottom layers [41,42] . Besides, as the core element for N cycle (e.g. N 2 -xation, nitri cation denitri cation and dissimilatory nitrate reduction to ammonium) [43] , Mo was noted to be positive related with N cycling index and diazotrophic community indices (Fig. 3B, S7). Mo de ciencymediated suppression in N 2 -xation capacity is a prevalent phenomenon in farmland soils, with a threshold of 0.15 mg kg − 1 for available Mo concentration [44] . In acidic Fe-accumuli soils, it is prone to form recalcitrant Mo induced by Fe-and aluminium precipitation [24,40] . As compared with Hapli-type soil, the lower level in available Mo content of Fe-accumuli soils (0.02 ± 0.01 mg kg − 1 ; Table S1) might also limit nitrogenase activity and thereby impede Mo-dependent N cycle.
The importance of microbial predictors to Fe and N cycling indices was indicated here across soil pro les, including diazotrophic and bacterial abundances, followed by microbial β-diversity. The predicting taxa for element cycling varied depending on soil depth, probably due to the vertical heterogeneity in environment factors (e.g. oxygen, moisture, light and nutrient poor) [27,38] . Some bacterial and diazotrophic predictors were responsible for both Fe and N cycling processes. In the bacterial community, Fe reducers (FeRB, e.g. Paenibacillus and Acidobacteria [45,46] ), Fe oxidizers (FeOB, Rhodanobacter) and sulfur-oxidizers (SOB, Thiobacillus and Dyella [47,48] ) were involved in N cycling process (Fig. 4A-4B), suggesting the interplays between N and Fe biogeochemical processes in paddy soils. For N 2 -xing bacteria that were important for predicting Fe cycling index, anaerobic FeRB (e.g. Anaeromyxobacter and Desulfotomaculum) [28,29] and autotrophic microaerophilic and neutrophilic FeOB (e.g. Leptothrix and Sideroxydans) [49] were identi ed. Although Fe bioavailability has been thought to govern diazotrophic distributions and activities due to its link to the activity of MoFe-dependent nitrogenase [50] , the evolvement of diazotrophs in Fe biogeochemistry in paddy soils has been largely overlooked. Besides, some Fe cycling related bacterial predictors, such as Gaiella (indicator for Hapli-type topsoil) and De uviicoccus, were reported to be signi cantly correlated with the chemodiversity of dissolved organic matter in paddy soil, although their mechanism underlying Fe transformation was obscure [51] .
Notably, the edaphic and microbial predictors (especially for diazotrophic community) differed markedly between the two soil types, indicating the distinctive biogeochemical processes (Fig. 4). For Fe-accumuli type, available Fe acted as the most important variable for controlling N cycling process, and N pool status (NO 3 − -N, NH 4 + -N and available N) was strongly correlated with Fe cycling index (Fig. 4B), indicting the coupling between N and Fe cycling in this soil type. Diazotrophic Geobacter was identi ed as the biomarker with higher relative abundances for Fe-accumuli type (Figure S10-S12), which was involved in both Fe and N cycling indices (Fig. 4B). The anaerobically respiring Geobacter are a known clade of active dissimilatory FeRB, which have been detected in rice paddy soil [52,53] . Geobacter-mediated Feammox (iron reduction coupled to anaerobic ammonium oxidation) has been shown to cause N loss from paddy soil [22] , which might explain for the negative link between Geobacter population and NH 4 + -N content in the anoxic topsoil of Fe-accumuli soils ( Figure S13A). Denitrifying bacterium Rhodanobacter was another differentiated indicator for Fe-accumuli soil along soil pro le ( Figure S9, S11). The dual-role played by this anaerobic species in N and Fe cycling (Fig. 4B) might be due to its nitrate-dependent Feoxidizing capacity, which dominated in acidic contaminated subsurface environment [54] . Besides, a higher proportion of aerobes were noted in diazotrophic community of Fe-accumuli soils as comparing with that in Hapli-type soils ( Figure S12), which might be linked with the NH 4 + -N storage in anaerobic paddy soils of Fe-type soil with weak N 2 xing ability. The phylogenetic a liation and formation mechanism of the presence of these aerobic diazotrophs will be further explored combined with omics and N 2 xation activity analysis (acetylene reduction assay and 15 N-isotope techniques) [23,24] .
On the other hand, available S best predicted N and Fe cycling indices in Hapli-type soil (Fig. 3B). Consistently, Hapli-stagnic anthrosol was characterized of large amounts of nifH-containing SRB across the pro les, including Desulfobacca and Desulfovibrio ( Figure S10, S12), which contributed signi cantly to both Fe and N cycling indices (Fig. 4B). Sulfur reducing and oxidizing microorganisms have been reported to act as active N 2 xers in ecosystems, such as in sul dic sediment [23,55,56] . The obligate anaerobe, Desulfobacca belonging to Syntrophaceae, has been identi ed as major acetate-degrading SRB and common N 2 xers in paddy soil [24,57] . In ooding soil, electron acceptors are reduced sequentially according to the thermodynamic theory, with the oxidative capacity of oxygen > NO 3 − > sulfate and Fe 3+ oxides [53] . Some SRB (e.g. Desulfovibrio and Desulfotomacufum) are also facultative FeRB or trigger Fe reduction indirectly via sul de production, serve as a dominant force in Fe biogeochemical cycling [30][31][32] .
Combined with previous reports that showed the S-driven Fe reduction coupled to anaerobic ammonium oxidation [21] , these data highlighted that the biogeochemical cycling of S, Fe and N could be coupled in Hapli-type anthrosol. To gain comprehensive and deeper understanding of mechanism coupling multielement biogeochemical processes in paddy soils, metagenomic approaches and hydrogeochemical evidences for element transformation should be further considered [58] .

Conclusions
The deep subsurface remains an environment where there is a limited understanding of microbial assembly and biogeochemical processes, especially for different paddy soils types with distinctive redox reactions. Here our results provided an integrated perspective on vertical assembly of soil bacterial and diazotrophic communities in two typical paddy soils types at a local scale, which were signi cantly in uenced by soil depth, followed by soil type. As compared with Hapli-stagnic anthrosol that experienced relative short-term ooding period, Fe-accumuli-anthrosol was characteristic of lower Fe cycling index in epipedons and N cycling index in the deep layers under long-term ood-drying management. By quantifying the contribution of microbial taxa to the differentiated soil N and Fe cycling indices, our data suggested the distinctive biogeochemical processes in the two paddy soil types, which were characteristic of iron reducing diazotroph Geobacter-mediated Fe-N cycling in Fe-accumuli type and sulfur reducing diazotroph Desulfobacca-mediated S-Fe-N cycling in Hapli-type soil. These ndings further highlighted the signi cant multiple role played by diazotrophs (e.g. Geobacter and Desulfobacca) in coupling multiple element cycling processes at the community-level.

Soil sampling
Soil samples were collected in October 2015 from paddy fields after rice harvest (Anhui Province, China; Figure 1A). For each site, one soil pro le was excavated to parent material (2 m Í 4 m, 80 to 120-cmdepth). Each soil pro le was divided into different horizons according to their pedological characteristics.
Detailed description and soil properties about the sampling sites was listed in Table S1. For each horizon, nine soil cores were randomly collected and every three soil subsamples were pooled into a composite sample. Samples were stored in 4°C in the eld and transported in a cooler to the laboratory immediately. In total, nine pro les containing 111 soil samples (37 subsamples × 3 replications) were collected. According to the diagnostic horizons, the nine paddy soil pro les were classi ed as two typical types (Feaccumuli-stagnic anthrosols and Hapli-stagnic anthrosols) based on the redoximorphic features and degree of paddy soil development [34,35] .

Soil properties
For each sample, soils were sieved through 2 mm meshes to remove visible plant tissues, stones and debris, homogenized and subdivided into three parts. One subsample was stored at -80°C for DNA 3') was used for the amplification of diazotrophic community [59] . PCR ampli cation and product puri cation were conducted based on previous studies [60] , and the amplicons were sequenced using the

Data analysis
Raw data from Illumina MiSeq high-throughput sequencing were analyzed using Quantitative Insight into Microbial Ecology (QIIME) 1.9.1 pipeline [9] . Low-quality data were removed (sequences with an average quality score < 25 for 16S rDNA and < 20 for nifH). Chimeras were removed using the Usearch8 in de novo mode. For 16S rDNA analysis, sequences were assigned to OTUs (operational taxonomy units) at 97% similarity by using UCLUST. The representative sequences for 16S rDNA were aligned referring to the SILVA database release 128 for bacteria at the 80% threshold [36] . For diazotrophs, nifH sequences were translated into amino acid sequences using the FunGene Pipeline of the Ribosomal Database Project in order to remove the failed translated sequences using FrameBot [60] . Sequences were clustered into operational protein units (OPUs) at 95% similarity by using UCLUST. All singleton OTUs/OPUs were deleted. To control heterogeneity of the number of sequences per sample, each sample of the data set was rare ed to minimal sequencing depth before to analyzing alpha and beta diversity indices. Statistical analysis were conducted using R software (Version 3.4.1; R Software for Statistical Computing, Vienna, Austria), with packages (e.g. 'vegan', 'picante' , 'ggcor' and 'labdsv'; details listed in Supplementary texts).

Assessing N and Fe cycling indices
Soil element cycling indices were quanti ed using the standardized average of edaphic properties with "multifunc" package [61] . N nutrient cycling index was calculated using total N, available N, NO 3 --N and NH 4 + -N, and Fe cycling index was evaluated by total Fe, available Fe and free Fe oxides [36,61] . For each index, related edaphic properties were normalized (log 10 -transformed as required) and standardized to a common scale using Z-score transformation [1] .

Random forest (RF)
Explanatory variables (environmental factors and microbial predictors) for soil type and nutrient cycling indices were predicted by RF analysis with "Random Forest" package [9] . The importance of each predictors was estimated by the percentage increases in the mean square error (MSE) of explanatory variable between observations and predictions. Signi cance of the model and cross-validated R 2 values were assessed with 5000 permutations of the response variables, using the "rfUtilities" package. The signi cance of predictor importance on the response variables was assessed with the "rfPermute" package.

Declarations
Note: The designations employed and the presentation of the material on this map do not imply the expression of any opinion whatsoever on the part of Research Square concerning the legal status of any country, territory, city or area or of its authorities, or concerning the delimitation of its frontiers or boundaries. This map has been provided by the authors. for the indicated factors. In each analysis, F, the P-value and the percentage of variation (%) explained by each factor refers to the total variance reported. Signi cance levels: ***, P < 0.001; **, P < 0.01; *, P < 0.05; ns, not signi cant. shown. Gene abundances were quanti ed by qPCR with speci c primers. Chao1 is used to represent the microbial richness. β-diversity indicates the rst principal coordinate of PLS-DA. Signi cance levels: *P < 0.05 and **P < 0.01.

Figure 4
The potential bacterial (A-B) and diazotrophic drivers (C-D) in predicting soil nitrogen (N) and iron (Fe) cycling indices in the whole pro le, super cial, the second (2nd) and deep layers (3rd to 5th), respectively, as revealed by random forest analysis. Signi cance levels: *P < 0.05 and **P < 0.01. Genera that are highlighted in red bold are enriched in Fe-accumuli type soil, while taxa that are highlighted in blue bold are abundant in Hapli-type soil.