Changes of microbial population and N-cycling function genes with depth in three Chinese paddy soils

Microbial communities play critical roles in soil nitrogen (N) cycle; however, we have limited understanding of the distribution of N-cycling microbial groups in deeper soil horizons. In this study, we used quantitative PCR to characterize the changes of microbial populations (16S rRNA and 18S rRNA) and five key N-cycling gene abundances involved in N fixation (nifH), ammonia oxidation (amoA) by ammonia-oxidizing bacteria (AOB) and ammonia-oxidizing archaea (AOA), and nitrite reduction (nirS and nirK) along profiles (0–100 cm depth) of different paddy soils from three regions (Hailun, Changshu, Yingtan) across China from north to south. We found that most microbial and N-cycling functional genes significantly decreased with soil depth; however, AOA were enriched in deeper soil layers (20–40 cm). The abundances of microbial and N-cycling functional genes generally decreased by one to two orders of magnitude in the deeper horizons relative to topsoils. The AOA gene abundance was higher than that of AOB in the paddy soil profile, and the nirS and nirK abundances were dominant in topsoil and deeper soil, respectively. All N functional genes except AOA were more abundant in Changshu than Hailun and Yingtan. High abundances and low vertical changes of N-cycling genes in Changshu suggest more dynamic N-transformations in this region. Correlation analysis showed that soil properties and climate parameters had a significant relationship with N-cycling gene abundances. Moreover, the abundance of different N-cycling genes was affected by different environmental parameters, which should be studied further to explore their roles in N cycling for sustainable agriculture and environmental management.


Introduction
Paddy ecosystems are essential components of agricultural systems, and support more than half of the world's population [1]. Soil microorganisms play an important role in paddy ecosystem, and can significantly affect soil fertility and rice productivity by participating in decomposition processes and nutrient cycles, particularly N cycling [2][3][4][5]. Nitrogen is generally the most limiting nutrient for rice production, and microbiological processes, e.g. nitrogen a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 fixation, nitrification and denitrification, were important in N cycling [6]. Biological nitrogen fixation can reduce atmospheric N 2 to bioavailable NH 4 + , and provide natural N input to cultivated rice cropping systems. Nitrification and denitrification are key processes that determine the efficiency of fertilizer use by rice crop, N loss from rice paddy soils, environmental pollutions such as nitrate leaching and emission of nitrous oxide [7]. However, we have limited understanding of microbial community distribution and functional characteristics in subsoil (below 0.2 m), as most studies have focused on the top 0.2 m of the soil column. Previous studies using DNA fingerprinting [8-10], phospholipid fatty acid [11][12][13] and 16S rRNA Miseq sequencing [13,14] have demonstrated that microbial community composition and diversity vary with soil depth, and microbial biomass declines exponentially with soil depth. However, few studies have attempted to comprehensively assess the change of microbial population and the abundance of nitrogen (N)-cycling function genes with soil depth.
Some key N-cycling genes encoded for enzymes involved in N-transformations including N-fixation (nifH), nitrification (bacterial and archaeal amoA), and denitrification (nirS and nirK) have been widely used to evaluate the functional abundance and diversity of soil microbial communities [15,16]. Previous studies have indicated that the abundances of soil Ncycling genes are associated with soil properties [17,18], soil type [19], and different N fertilizer regimes [20]. The abundances of N-cycling genes vary with soil depth, but the variability is different for different N-cycling genes. Leininger et al. (2006) found that both AOA and AOB gene abundances decreased with soil depth, but AOB amoA gene abundance decreased more strongly with depth than AOA at 0-40 cm soil depth [21]. However, in some wetland sites, the gene abundances of both AOA and AOB increased with soil depth at 0-30 cm depth, especially at high-ammonium sites [22]. Some studies found that nirS genes were more abundant than nirK genes [23], but other studies have shown the opposite trend [24]. These studies suggested that N-cycling genes require further investigation.
China is the world's largest rice producer [25], and rice is planted in a large latitude range and on diverse soil types in China. In this study, we selected three typical paddy soils with distinct soil properties across a latitude gradient from northern to southern China to examine the shifts of microbial population and N-cycling gene abundance with soil depth (0-100 cm) in three different paddy soils with the q-PCR method, and to explain the environmental factors controlling this variation.

Materials and methods
I state clearly that no specific permissions were required for these locations/activities, which are concerned about scientific study in these regions in China. I confirm that the field studies did not involve endangered or protected species.

DNA extraction
Soil DNA was extracted from 0.5 g of freeze-dried soil samples using the FastDNA SPIN Kit for soil (MP Biomedicals, Solon, OH, USA), according to the manufacturer's instructions. The DNA quality was assessed in agarose gel (1%) and DNA quantity was measured with a Nanodrop spectrophotometer (NanoDrop, Wimmington, DE, USA). The extracted DNA was stored at −20˚C for subsequent quantification of microbial population and functional genes.

Quantitative PCR
The abundances of bacterial 16S rRNA gene, fungal 18S rRNA gene and N-cycling functional genes (nifH, AOA-amoA, AOB-amoA, nirS and nirK) were measured by qPCR using SYBR1 Premix Ex Taq™ II (Takara, Japan) on a LightCycler1 96 Real-time PCR System (Roche, Switzerland). The qPCR reaction (20 μl) contained 10 μl of SYBR Premix Ex Taq™ II, 0.4 μl each of the forward and reverse primers ( To obtain the standard curves for qPCR assays, the bacterial 16S rRNA gene, fungal 18S rRNA genes and functional genes (AOA-amoA, AOB-amoA, nifH, nirS and nirK) were amplified with the primers listed in Table 2. The PCR products from each gene were purified with a Gel Extraction Kit (CW Biotech co, Beijing, China) and cloned into pMD19-T vector (Takara,

Statistical analysis
Statistical analyses of the effects of region and soil depth on microbial population and N-functional genes were conducted by two-way analysis of variance (ANOVA) using Tukey's honest significant difference (HSD) multiple comparison with a P = 0.05 grouping baseline. All of the statistical analyses were performed with SPSS 16.0 (SPSS Inc., Chicago, IL, USA). Spearman correlation analysis was employed to determine the significant correlations between the abundances of microbial population and N-cycling genes and environmental parameters. The correlation coefficients were calculated and plotted using the "corrplot" package in R software (version 3.3.1). The probability level P < 0.05 was considered to be statistically significant.

Abundance of bacteria and fungi
The bacterial 16S rRNA gene copy numbers in all three regions ranged from 2.70 × 10 10 to 4.31 × 10 8 copies g −1 dry soil, which were more abundant than fungal 18S rRNA gene copy numbers, ranging between 3.56 × 10 8 and 5.71 × 10 6 copies g −1 dry soil (Fig 1A and 1B). The ANOVA revealed that soil depth, region, and their interaction all had significant effects on bacteria and fungi population (P < 0.001, Table 3). Bacterial 16S rRNA gene numbers slightly decreased at 0-20 cm soil depths, then significantly decreased at 20-100 cm soil depths in all nirK-R3Cu nirK GCCTCGATCAGRTTGTGGTT 57 three regions but especially Yingtan. The abundance of 16S rRNA gene throughout soil profile generally decreased in the order Hailun > Changshu > Yingtan. The bacterial 16S rRNA gene copy numbers in Changshu ranged from 1.88 × 10 10 to 3.55 × 10 9 copies g −1 dry soil, and reduced by 81.11% with depth. The ranges for Hailun and Yingtan were 2.70 × 10 10 to 5.44 × 10 9 and 1.94 × 10 10 to 4.31 × 10 8 copies g −1 dry soil, and decreased by 79.85% and 97.77% with depth, respectively. Fungal 18S rRNA gene copy numbers significantly decreased at 0-60 cm depths in the three regions, especially in the Hailun soil, and there were no significant differences among the 40-60, 60-80, and 80-100 cm depths. The abundance of 18S rRNA gene throughout soil profile generally decreased in the order Changshu > Yingtan > Hailun. The fungal 18S rRNA gene copy numbers in Changshu ranged from 3.56 × 10 8 to 2.66 × 10 7 copies g −1 dry soil, and reduced by 92.52% with depth. The ranges for Hailun and Yingtan were 1.78 × 10 8 to 5.71 × 10 6 and 2.30 × 10 8 to 1.37 × 10 7 copies g −1 dry soil, and decreased by 96.79% and 94.04% with depth, respectively. The abundances of both 16S rRNA and 18S rRNA genes were generally one to two orders of magnitude lower in the deeper horizons than the topsoils. The ratio of 18S rRNA to 16S rRNA gene abundance generally decreased at 0-60 cm depth and increased at 60-100 cm depth in the three regions (Fig 2A). At 0-60 cm soil depth, the ratios of 18S rRNA/ 16S rRNA decreased from 0.6-1.9% to 0.05-0.7%, and then the ratios  Table 3. Two-way ANOVA analysis of the abundances of bacteria, fungi and N-cycling functional genes at six soil depths in three sites each with five replicates (n = 90). The data in italics indicated that abundance of N-cycling gene was affected by soil depths, sites or their interaction. 1.18 × 10 6 copies g −1 dry soil, and reduced by 96.30% with depth. The ranges for Hailun and Yingtan were 6.20 × 10 7 to 2.21 × 10 6 and 2.10 × 10 7 to 1.31 × 10 5 copies g −1 dry soil, and decreased by 96.43% and 99.37% with depth, respectively. The AOB amoA gene copy numbers significantly decreased at depths between 0 and 40 cm depths (especially 20-40 cm depth) in all three regions, but there was no significant difference among the 40-60, 60-80, and 80-100 cm depths. The AOB amoA gene copy numbers in Changshu ranged from 2.44 × 10 6 to 2.20 × 10 4 copies g −1 dry soil, and reduced by 99.09% with depth. The ranges for Hailun and Yingtan were 6.72 × 10 5 to 6.72 × 10 3 and 4.33 × 10 5 to 2.18 × 10 4 copies g −1 dry soil, and decreased by 99.00% and 94.96% with depth, respectively. The abundance of both AOA and AOB amoA genes generally decreased by one to two orders of magnitude along soil depth gradients.

Bacteria
The ratio change of AOA to AOB amoA gene abundances was similar with the change of AOA amoA gene abundance with soil depths (Fig 4A). The ratio increased at 0-40 cm soil depth, and then decreased at 40-100 cm soil depths in Hailun and Changshu soils, but the ratio in Yingtan generally decreased with increasing soil depth. The ratio peaked at 20-40 cm depth, and AOA amoA gene abundances were two to five orders of magnitude higher than AOB amoA gene in soil profiles.

Denitrification
The nirS gene copy numbers ranged from 9.72 × 10 8 to 1.11 × 10 6 g −1 dry soil among the three regions and nirK gene abundances were from 4.85 × 10 8 to 2.85 × 10 6 g −1 dry soil (Fig 5A and  5B). Two-way ANOVA showed that soil depth, region, and their interaction all had significant effects on the nirS and nirK gene copy numbers ( Table 3). The changes of nirS and nirK gene abundances with soil depth in all three regions showed a similar trend in which the copy numbers of both genes generally decreased with increasing soil depths. The abundances of both genes throughout soil profile generally decreased in the order Changshu > Hailun > Yingtan. The nirS gene abundances in Changshu ranged from 9.72 × 10 8 to 5.79 × 10 7 copies g −1 dry soil, and reduced by 94.04% with depth. The ranges for Hailun and Yingtan were 7.43 × 10 8 to 6.57 × 10 6 and 6.75 × 10 8 to 1.11 × 10 6 copies g −1 dry soil, and decreased by 99.11% and 99.83% with depth, respectively. The nirK gene abundances in Changshu ranged from 4.85 × 10 8 to 1.84 × 10 8 copies g −1 dry soil, and reduced by 62.06% with depth. The ranges for Hailun and Yingtan were 1.73 × 10 8 to 3.29 × 10 7 and 8.39 × 10 7 to 2.85 × 10 6 copies g −1 dry soil, and decreased by 80.98% and 96.60% with depth, respectively. The abundance of the nirS gene generally decreased by one to two orders of magnitude along soil depth gradients. In contrast, the abundance of the nirK gene generally decreased by one order of magnitude along soil depth gradients.
The ratios of nirS to nirK gene abundances generally decreased with increasing soil depths, except for the 10-20 cm depths in Changshu soil, where the ratio slightly increased (Fig 4B). The ratio ranged from 9.04 to 0.24. The nirS gene abundances were higher than those of the nirK gene at depths between 0 and 20 cm in Changshu and Hailun soils and between 0 and 60

Relationships between microbial population, functional gene abundances and environmental parameters
Spearman's correlation coefficients were used to compare environmental parameters and the abundances of bacterial 16S rRNA, fungal 18S rRNA and N-cycling functional genes. In Fig 6, all environmental parameters were shown in the correlation analysis across the three regions, but for the correlation analysis of single regions, only the environmental factors that significantly correlated with the abundances of bacterial 16S rRNA, fungal 18S rRNA or N-cycling functional genes are presented. The correlation analysis across three regions, integrating the influences of depth and region, revealed that soil physical structure (sand and clay) correlated with all microbial and N-cycling functional gene abundances, except nirK and clay. In addition to soil physical structure, both nirS and nirK gene abundances were significantly correlated with elevation and pH. Furthermore, nirK, AOA and 16S rRNA gene abundances had more significant correlations with environmental parameters than nirS, AOB and 18S rRNA, and were positively correlated with total N, total carbon (C), C/N ratio and negatively correlated with annual precipitation and frost-free period. The 18S rRNA gene abundances were negatively correlated with elevation. In addition, nirK gene abundance was also correlated with soil water content and dissolved organic carbon (DOC) concentration, and AOA gene abundance positively correlated with nitrate concentration. For the correlation analysis of single regions, results similar to the total correlation analysis across three regions were shown (soil physical structure (sand and clay) correlated with most of the microbial and N-cycling functional gene abundances). In addition to soil physical structure, the changes of nirS gene abundances along the soil profile negatively correlated with DOC concentration, and the abundance changes of 18S rRNA, nifH and AOA genes were significantly positively correlated with the C/N ratio.

Discussion
Shifts in the abundances of bacteria and fungi along the soil profile across the three regions Both bacteria and fungi were one to two orders of magnitude higher in the upper 10 cm of soil than at 1 m depth. Decreasing microbial biomass with soil depth has been previously observed in similar studies [34][35][36]. Throughout the soil profiles, we observed significant differences in nutrient concentrations (C and N), soil moisture, pH, and soil texture with soil depth (Table 1). Both bacteria and fungi in surface soils are likely to have greater nutrient resources and environmental conditions that favor their growth. Although both bacterial and fungal abundance decreased with soil depth, fungal abundance decreased more dramatically than bacteria at depths between 0 and 60 cm. This is similar to the findings of previous study about community structure that indicated fungal communities were more strongly influenced by soil depth than bacterial communities [37]. We observed that the change of fungi abundance along soil depth was significantly correlated with soil C/N ratio, and the results are broadly consistent with work by Lauber [38], who found that fungal community was most closely associated with soil nutrient status and the soil C/N ratio. At the large scale across the three regions, the abundances of both bacteria and fungi were significantly different among the three regions. This supports the findings of other studies that soil type or site conditions are also the determining factors of changes in bacterial and fungal communities [38][39][40][41].
Shifts in the abundances of N-cycling functional groups along the soil profile across the three regions In terrestrial ecosystems, the inputs of N to soil mainly depend on biological N fixation. Nitrogen-fixing bacteria use nitrogenase enzymes to catalyze this process. The nifH gene is widely used as a marker for molecular analysis of N-fixing bacteria [42]. In our study, the abundance of the nifH gene substantially declined with depth in the three soil profiles, which is consistent with the findings of other studies [18,43]. Moreover, we observed that the abundance changes of the nifH gene with soil depth were significantly correlated with soil C/N ratio. Some studies have previously indicated that C/N is an important factor influencing the abundance of N-fixing bacteria [44,45]. The increase in C/N ratio that C-rich and N-depleted conditions might provide a competitive advantage for N-fixing bacteria.
Nitrification, the microbial oxidation of ammonia to nitrite and nitrate, is a key process in the global cycling of N. The first and rate-limiting step is catalyzed by the enzyme ammonia monooxygenase (AMO), which exists in both AOB and AOA [46]. In our study, the abundance of AOA was greater than that of AOB in all the soil samples, and only AOA abundances were significantly positively correlated with the ammonium content (substrate of ammonia oxidation), which suggested that AOA are more important nitrifiers in paddy soils than AOB in our study. Our results were in accordance with previous studies, which showed that AOA were the dominant group among ammonia-oxidizing prokaryotes in paddy soil [47,48]. Furthermore, some researchers have previously indicated that AOA are the main drivers of ammonia oxidation in soil by measuring transcriptional activity of AOA amoA genes or based on analysis with stable isotope probe [49,50]. Correlation analysis illustrated that soil total C, total N, C/N ratio, pH, and ammonium content were correlated with abundance of AOA rather than AOB. This suggests that AOA were more sensitive to soil characteristics than AOB, which is consistent with the findings of previous studies [51,52]. Our study indicated that environmental factors might be important in affecting the relative abundances of AOA and AOB.
Denitrification is of major importance to the N cycle and consists of four reaction steps in which nitrate is reduced into dinitrogen gas. The reduction of nitrite to nitric oxide is catalyzed by two structurally different nitrite reductases (Nir): one contains copper (Cu-Nir) encoded by the nirK gene while the other contains heme c and heme d1 (cd1-Nir) encoded by the nirS gene [53]. In this study, denitrification genes (nirS and nirK) decreased significantly with soil depth at all three sites. Both genes accounted for the largest proportion of measured N-cycling genes, which is in accordance with previous findings that their proportions in bacterial community ranged from 0.5 to 5% [54,55]. In our study, the abundance of nirS was greater than that of nirK in 0-20 cm surface paddy soil, but in subsoil nirK abundance was higher than nirS, which suggests that nirK-harboring denitrifiers thrived in the deeper soil layers rather than nirS-harboring denitrifiers. Correlation analysis illustrated that total C, total N, soil water content, C/N ratio, pH, DOC concentration, and sand and clay content were correlated with nirK or nirS abundances among soil properties, and soil pH and sand content had the strongest correlations with nirK gene abundance (P < 0.01). Previous studies have shown that nirS denitrifiers are located mostly in the rhizosphere, whereas nirK denitrifiers are more abundant in the bulk soil [56,57]. This suggests that soil nutrient resource availability is an important factor influencing the denitrifier community. Clark et al. (2012) found that soil properties (soil pH, Fig 6. Correlogram representing spearman's correlation coefficient rank between environmental parameters and abundances of microbial and of N cycling genes. The correlation coefficients ranging from negative to positive are indicated by color intensity changing from red to blue. *P < 0.05, **P < 0.01. clay, and C and N contents) had a major influence on denitrifier communities [58]. In addition, Bárta et al. (2010) found nirK gene abundance was higher in more alkaline soils and nirS abundance was higher in more acidic soils [59]. Dini-Andreote et al. (2016) indicated that nirK gene abundance was positively and significantly correlated with soil pH and sand content (p = 0.70 and 0.69, respectively; P < 0.01) [60]. In this study, resource availability, pH and soil physical properties showed significant changes along the soil depth gradient or among the three regions, which might represent an important effect on the distribution of nirS and nirKcontaining denitrifiers.

Linking the environmental factors and N-cycling abundance
In our study, not only soil properties, but also climate parameters were significantly correlated with N-cycling gene abundances. Previous studies have indicated that soil physicochemical changes are important factors influencing the distribution of N-cycling functional groups including pH [61,62], dissolved oxygen levels [59,63], nutrient resource [60,61], salinity [64], C/N ratio [44,61], and soil texture [58,65]. However, there has been little research focusing on the distribution of N-cycling microbes at a large scale. Hu et al. (2015) found that the biogeographic distributions of both AOA and AOB were significantly regulated by spatial factors (latitude), climatic factors (mean annual temperature and precipitation), and geochemical factors (soil pH, total N, sulfate, clay%, and C/N ratio) [66]. Bru et al. (2011) reported that climate explained around 10-13% of the variability in the relative abundances of nirS and nirK denitrifiers within the total bacterial community when the distribution of N-cycling functional microbial communities was investigated over 107 sites in Burgundy, a 31,500 km 2 region of France [55]. In this study, we found that climate parameters and elevation were significantly correlated with AOA, nirK, and nirS gene abundances at the large scale, which, together with soil properties, affected the changes of N-cycling gene abundances.
In this study, we used q-PCR method to characterize the changes of microbial populations and N-cycling gene abundances in the paddy soil profile. However, this method only detected the abundances of specific genes that we have designed to target. In the future, more novel metagenomics and metatranscriptomics approaches will enable us to connect specific bacterial and fungal populations to specific nitrogen cycling processes.

Conclusions
In this study, we examined changes of microbial populations (16S rRNA and 18S rRNA) and five key N-cycling gene abundances along three soil profiles in paddy fields. Our results showed that the abundances of microbial and N-cycling functional genes (except AOA) generally decreased significantly along the soil depth gradient; however, the AOA were enriched in the deeper soil layers (20-40 cm). The abundance of AOA genes was higher than those of AOB throughout the paddy soil profile, and the abundances of nirS and nirK genes were dominant in topsoil and deeper soil, respectively. Soil properties and climate parameters had significant relationships with N-cycling gene abundances, which should be explored further to assess their roles in N cycling for applications in sustainable agriculture and environmental management.