Nitrification inhibitors effectively target N2 O-producing Nitrosospira spp. in tropical soil.

The nitrification inhibitors (NIs) 3,4-dimethylpyrazole (DMPP) and dicyandiamide (DCD) can effectively reduce N2 O emissions; however, which species are targeted and the effect of these NIs on the microbial nitrifier community is still unclear. Here, we identified the ammonia oxidizing bacteria (AOB) species linked to N2 O emissions and evaluated the effects of urea and urea with DCD and DMPP on the nitrifying community in a 258 day field experiment under sugarcane. Using an amoA AOB amplicon sequencing approach and mining a previous dataset of 16S rRNA sequences, we characterized the most likely N2 O-producing AOB as a Nitrosospira spp. and identified Nitrosospira (AOB), Nitrososphaera (archaeal ammonia oxidizer) and Nitrospira (nitrite-oxidizer) as the most abundant, present nitrifiers. The fertilizer treatments had no effect on the alpha and beta diversities of the AOB communities. Interestingly, we found three clusters of co-varying variables with nitrifier operational taxonomic units (OTUs): the N2 O-producing AOB Nitrosospira with N2 O, NO3 - , NH4 + , water-filled pore space (WFPS) and pH; AOA Nitrososphaera with NO3 - , NH4 + and pH; and AOA Nitrososphaera and NOB Nitrospira with NH4 + , which suggests different drivers. These results support the co-occurrence of non-N2 O-producing Nitrososphaera and Nitrospira in the unfertilized soils and the promotion of N2 O-producing Nitrosospira under urea fertilization. Further, we suggest that DMPP is a more effective NI than DCD in tropical soil under sugarcane.

Introduction

Anthropogenic inputs of N fertilizers to agriculture have stimulated agricultural soils to contribute up to 59% of anthropogenic N2O emissions (Fields, 2004; Robertson and Vitousek, 2009; Ciais, 2013; Signor and Cerri, 2013). Because N2O has a global warming potential 298 times that of CO2 (Ravishankara et al., 2009) and diverts N that would otherwise be used by the crop, reducing N2O emissions is a major target for sustainable management practices (Venterea et al., 2012). The N2O emitted from a soil is the cumulative result of abiotic and biotic N2O‐generating pathways (Graham et al., 2014; Hu et al., 2015). The two main biotic processes contributing to N2O in agricultural soils are nitrification (oxidation of NH4 + to NO2 − to NO3 −) and denitrification (reduction of NO3 − to NO2 − to N2O to N2; reviewed in (Guo et al., 2013; Ruser and Schulz, 2015). Nitrification is carried out by a few bacterial and archaeal genera; ammonia oxidation is mediated by the ammonia‐oxidizing archaea (AOA), such as the Thaumarchaeota Nitrososphaera and Nitrosopumilus (Walker et al., 2010; Offre et al., 2013) and the ammonia‐oxidizing bacteria (AOB), such as the Betaproteobacteria Nitrosomonas and Nitrosospira; nitrite oxidation is carried out by nitrite oxidizing bacteria (NOB), including the Nitrospirae Nitrospira and the Alphaproteobacteria Nitrobacter. Denitrification is carried out by microorganisms widely dispersed over the bacterial, archaeal and fungal domains and denitrification genes can also be carried by nitrifiers in what is termed nitrifier denitrification (Kool et al., 2010). Further, the process of complete nitrification by the recently discovered comammox bacteria, which have so far been found in the NOB Nitrospira genus, might also contribute to N2O emissions (Liu et al., 2017).

Nitrification and denitrification processes are regulated by the abiotic factors temperature, oxygen availability, moisture, ammonia and nitrate availability, carbon availability and pH (Wallenstein et al., 2006; Butterbach‐Bahl et al., 2013). These factors also affect the distribution and niche differentiation of nitrifiers; for example, the AOB numerically dominate in neutral soils with high NH4 + concentrations while the AOA numerically dominate in acidic soils with low NH4 + concentrations (Di et al., 2009, 2010a). These general ranges can vary within taxonomic groups; the AOB Nitrosospira are more common in acid soils (Pommerening‐Röser and Koops, 2005) compared to the AOB Nitrosomonas, which are uncommon in acidic environments (pH 4–5) (Song et al., 2016; Li et al., 2018). The AOB Nitrosospira isolated from acid soils in general have urease enzymes catalysing the breakdown of urea to ammonia (De Boer and Kowalchuk, 2001), and these ureolytic AOB characteristics allow them to grow at low pH with urea source (Pommerening‐Röser and Koops, 2005). However, there are also exceptions to the general rule, for example, a Gammaproteobacteria AOB, Candidatus Nitrosoglobus, was recently isolated from acidic soils with survival in conditions down to pH 2 (Hayatsu et al., 2017). Further, the nitrite oxidizer bacteria Nitrobacter and Nitrospira have optimal growth under higher and lower nitrite supplies, respectively, which is linked to their ecological niches (Attard et al., 2010; Nowka et al., 2015).

Nitrification is doubly implicated in N2O production, either directly or indirectly by producing NO3 − as the basis for denitrification, and has been shown to be the main process involved in N2O emissions in some Brazilian sugarcane soils (Liu et al., 2016; Soares et al., 2016; Wu et al., 2017a, 2017b; Lourenço et al., 2018a, 2018b). The addition of nitrification inhibitors with nitrogen fertilizers is currently being explored as a sustainable management practice in Brazilian sugarcane (Signor et al., 2013; Soares et al., 2015, 2016). In agriculture, dicyandiamide (DCD) and 3,4‐dimethylpyrazole phosphate (DMPP) are commercially‐used nitrification inhibitors which are thought to be Cu‐chelating agents acting on ammonia monooxygenase (Morales et al., 2015). These inhibitors have been shown to effectively reduce N2O emissions by 40%–95% in temperate and tropical soils (Misselbrook et al., 2014; Gilsanz et al., 2016; Soares et al., 2016).The effects of DCD and DMPP are generally shown to act on ammonia oxidizing bacterial abundances, perhaps because this group increases with increasing N fertilization and likely contributes to the N2O emissions under N fertilization in these studies (Morales et al., 2015; Soares et al., 2016). Beyond lowering N2O emissions, nitrogen fertilization with nitrification inhibitors might have complex effects on ammonia and nitrite oxidizer distributions (Stempfhuber et al., 2016). It is yet unknown which nitrifiers are specifically affected by nitrogen fertilization with the nitrification inhibitors DCD and DMPP, especially in tropical conditions.

Here, our objectives were to identify the AOB species linked to N2O emissions in our previous experiment (Soares et al., 2016) and to compare the effects of urea fertilization with or without nitrification inhibitors on amoA‐containing nitrifier abundances. We sequenced amoA AOB amplicons and mined the previous 16S rRNA amplicons from a 258 day field experiment encompassing treatments with urea and two nitrification inhibitors, DCD and DMPP, on soils growing ratoon sugarcane. To gain insight into the environmental niches of different nitrifiers, wee included the previously generated 16S rRNA gene, nitrification and denitrification gene copy numbers and soil environmental variable datasets (Soares et al., 2016). To our knowledge, no studies to date have examined the dynamics of the nitrifier community growing in tropical soil nitrification inhibitors in urea fertilized soils.

Results

amoA AOB community sequencing coverage and composition

Processing of the amoA AOB amplicon data resulted in 68 211 sequences, which were clustered into 54 OTUs. The number of sequences ranged between 121 and 3019 across the 127 samples (4 treatments × 8 time points × 4 replicates with one outlier sample removed). The samples had average Good's coverage of at least 94% (Supporting Information Table S1), which was supported in the rarefaction curve results, with more sequences not adding more species in the samples (Supporting Information Fig. S3). At the genus level, the AOB community was composed of unclassified Betaproteobacteria, Nitrosomonas and Nitrosospira (which included the Nitrosovibrio classification; Fig. 1A). The phylogenetic tree of the amoA AOB OTUs with reference sequences indicated that these aligned with Nitrosospira (52/54 amoA OTUs) and Nitrosomonas (2/54 amoA OTUs; Fig. 2). In support of the low diversity of the amoA AOB communities, the 16S rRNA gene dataset revealed only two Nitrosospira OTUs (abundant OTU 30 and rare 16S rRNA OTU 1102) and one Nitrosomonas OTU (rare 16S rRNA OTU 2875). Further, the Nitrosospira 16S rRNA OTUs had similar absolute abundances as the Nitrosospira amoA AOB OTUs across the treatments (Fig. 1B and D), supporting the idea that both amplicon datasets adequately covered the AOB communities in the samples.

Figure 1

Taxonomic distributions of the amoA AOB amplicon samples by (A) relative abundances or (B) absolute abundance estimates within genus and (C) the amoA AOB gene copy numbers. Also included are the (D) taxonomic distributions of the 16S rRNA amplicon samples by absolute abundance estimates within the Nitrosomonadaceae family and (E) the gene copy numbers of 16S rRNA gene sequences. Mean values within treatments and time points are shown. Treatments were the unfertilized control (C), urea (UR), urea with dicyanimide (UR + DCD) and urea with 3,4‐dimethylpyrazole phosphate (UR + DMPP). Day = days after fertilization. [Color figure can be viewed at http://wileyonlinelibrary.com]

Figure 2

Phylogenetic analysis of amoA AOB OTUs and reference amoA sequences from the FUNGENE database based on the maximum likelihood distance method. Bootstrap values (1000 replicates) of less than 90% are depicted by the blue dots on the branches. The Nitrosospira are depicted with purple colour, Nitrosomonas with green colour and the outgroup Nitrososoccus with the yellow colour. [Color figure can be viewed at http://wileyonlinelibrary.com]

Treatment effects on amoA AOB community beta diversity

Beta dispersion analysis on all the samples revealed that treatment, but not time point, had a significant effect on the AOB community dispersions (F = 3.6529, p < 0.05). Subsequent beta diversity analysis revealed that time point, considering all treatments, had no effect on the AOB community structures (Supporting Information Table S2). Ordination plots showed that the amoA AOB communities overlapped between treatments, considering all time points, according to 95% confidence intervals (Fig. 3A). Within time points and treatments, the beta dispersions of the amoA AOB communities were unaffected by treatment and time point, respectively. Treatment had a significant effect on the amoA AOB community structures only within Days 7 and 16 (PERMANOVA; p < 0.1; Supporting Information Table S2). However, pairwise comparisons revealed that no amoA AOB community structures were significantly different between treatments within these time points. Time point had no effect on amoA AOB community structures within any treatment. Ordination plots within time point revealed that the amoA AOB communities did not cluster separately for treatments nor time points at 95% confidence intervals (Fig. 3C and D).

Figure 3

Ordination plots of the ammonia‐oxidizing bacterial communities using PCoA on Bray–Curtis sample distances based on amoA AOB OTU relative abundances (n = 127) (A) across all time points (n = 127) and (B) across all treatments, (C) by treatment within each time point (n = 16) and (D) by time point within each treatment (n = 24). Time points were 7, 16, 18, 27, 35, 42, 82 and 158 days after fertilization. Treatments were unfertilized (Control, black), urea (UR, red), urea with dicyanimide (UR + DCD, purple), urea with 3,4‐dimethylpyrazole phosphate (UR + DMPP, green). Confidence intervals of 0.95 are drawn around the treatments or days as ellipses. [Color figure can be viewed at http://wileyonlinelibrary.com]

Treatment and time point effects on amoA AOB community alpha diversity

The alpha diversities of the amoA AOB communities ranged from 1 to 3 based on Shannon index (Fig. 4). Considering all time points, treatment had an effect on the alpha diversity of the amoA AOB communities (chi‐squared value 33.884, p value = 2.096e−07), but time point had no effect on the alpha diversities when considering all treatments. Post hoc testing over all time points found that the amoA AOB communities in the DMPP treatment had higher alpha diversity compared to the other treatments (Dunn's test, p < 0.05; Fig. 4). Within time point, treatment had an effect on the amoA AOB alpha diversities for Days 7, 18, 27, with chi‐squared values of 7.6103 (p value 0.05479), 4.7792 (p value 0.1887) and 6.7721 (p value 0.07953), respectively. However, post hoc testing revealed no different pairs. Within treatment, time point had an effect on the amoA AOB community alpha diversities only for the Control treatment (chi‐squared 12.534, p value = 0.08431); further, pairwise post hoc tests revealed no difference in alpha diversity between treatments.

Figure 4

Alpha diversity of the amoA AOB communities calculated from rarefied raw abundances as affected by (A) treatment, for all time points, (B) time point, for all treatments, (C–J) treatment, within each time point and (K–N) time point, within each treatment.Treatments were unfertilized (C), urea (UR), urea with dicyanimide (UR + DCD) and urea with 3,4‐dimethylpyrazole phosphate (UR + DMPP); time points were 7, 16, 18, 27, 35, 42, 82, 158 days after fertilization. The y‐axis label includes the result of a Kruskal–Wallis chi‐squared test (‘*’ for p < 0.05, ‘.’ for p < 0.10); the letters above the plots represent the results of Dunn's post hoc tests at alpha <0.05 in which similar letters denote no difference between groups.

Differential abundance of nitrifier 16S rRNA OTUs and treatment group indicators

From the 16S rRNA variance‐stabilized abundances, four genera of nitrifiers were represented: Nitrosomonas (1 OTU), Nitrososphaera (37 OTUs), Nitrosospira (2 OTUs) and Nitrospira (11 OTUs). The variance‐stabilized trajectories of Nitrosospira, Nitrososphaera and Nitrospira 16S rRNA OTUs across the four treatments can be seen in Supporting Information Fig. S4. The two 16S rRNA Nitrosospira OTUs showed a similar trend across the treatments, with higher abundances in the urea and urea with DCD treatments compared to the control and the urea with DMPP treatments. The 16S rRNA Nitrososphaera OTUs showed three trends, with OTUs 11 and 429 having lowest abundances in the control treatment and higher abundances in the treatments with urea, with the highest abundances in the urea with DMPP treatment; OTUs 40 and 45 having highest abundances in the control treatment, lower abundances in the treatments with urea and the lowest abundance in the urea treatment; and OTUs 112 and 39 having highest abundances in the control and urea with nitrification inhibitor treatments and the lowest abundance in the urea treatment. The 16S rRNA Nitrospira OTU followed the last trend with the highest abundances in the control and urea with nitrification inhibitor treatments and the lowest abundance in the urea treatment.

The results of the differential abundance and indicator species analyses, based on the variance‐stabilized abundances and absolute abundance estimates, respectively, generally agreed. Differential abundance analysis revealed the nitrifier 16S rRNA OTUs that were significantly over‐ and under‐represented between pairwise comparisons of treatments based on variance‐stabilized abundances (Supporting Information Table S3). Of the Nitrosospira 16S rRNA OTUs, OTU 30 was an indicator of the control, urea and urea with DCD treatments, while OTU 1102 was an indicator of only the urea and urea with DCD treatments (adjusted p value < 0.1; Supporting Information Table S3). Of the Nitrososphaera 16S rRNA OTUs, OTU 45, OTU 112, OTU 40, OTU 39 and OTU 11 were indicators of the control, urea with DCD and urea with DMPP treatments. Of the Nitrospira 16S rRNA OTUs, OTU 79 was an indicator of the control, urea with DCD and urea with DMPP treatments.

Nitrifier amoA and 16S rRNA OTU and environmental correlations

The correlations of the environmental variables with the gene copy numbers of AOB, AOA, nirK, nirS, 16S rRNA total bacteria and total Archaea (Supporting Information Fig. S5) of all eight sampling days depict the positive links between AOB, N2O, NO3 −, NH4 +, WFPS and pH and AOA, nirS, nirK, total archaea and total bacteria; and the negative links between CO2, CH4 + and WFPS and AOA, total Archaea, NH4 +, N2O and NO3 − (Fig. 5A). As can be seen in Fig. 5B which depicts correlations including the normalized abundances of 16S rRNA OTUs, N2O emissions were correlated with amoA AOB copy numbers, water‐filled pore space (WFPS), NO3 −, NH4 + and pH. Interestingly, the 16S rRNA and amoA AOB OTU correlations clustered with the previous variables with the exception of NH4 + and pH, which nevertheless suggests that Nitrosospira (16S rRNA OTU 30 and OTU 1102) were the N2O‐producing AOB in these soils. Other interesting clusters were the 16S rRNA Nitrososphaera OTUs 429 and 11 with NO3 −, NH4 + and pH; the nirS, nirK, total archaeal and 16S rRNA gene copy numbers; and the amoA AOA, 16S rRNA Nitrospira OTU 79, the 16S rRNA Nitrososphaera OTUs 45, 112, 40 and 39. These clusters were found in all the correlations with absolute abundance estimates and relative abundances of the amoA AOB and the 16S rRNA gene sequence data (Supporting Information Fig. S5).

Figure 5

Cluster plot visualizing Spearman's correlations (A) between environmental variables and gene copy numbers and (B) between environmental variables, gene copy numbers and the variance‐stabilized abundances of the 16S rRNA nitrifier OTUs. Normalization was carried out using DeSeq2. Only significant correlations are shown (p < 0.01). Clusters were determined using complete linkage clustering. NOS = Nitrosospira, NOP = Nitrososphaera, NOB = Nitrospira. [Color figure can be viewed at http://wileyonlinelibrary.com]

Discussion

From our previous work, we found that bacterial amoA (AOB) but not archaeal amoA (AOA) nor denitrification gene copy numbers (nirK, nirS) were correlated with nitrous oxide emissions from tropical soil growing sugarcane (Soares et al., 2016). Here, we found evidence that the AOB responsible for the N2O emissions was most phylogenetically similar to the Nitrosospira spp. (Nitrosovibrio RY3C), based on the decrease in abundance of these OTUs in soils with the nitrification inhibitors in comparison with the urea treatment and the correlation of these OTUs with N2O emissions. The Nitrosovibrio RY3C species was originally isolated from avocado rhizosphere and its nitrifying activity was susceptible to DCD (Matsuba et al., 2003). To our knowledge, just one other study has identified Nitrosospira spp. as the N2O‐generating AOB in tropical soil under sugarcane, and that study applied NH4NO3 as the N source (Lourenço et al., 2018a). Interestingly, Lourenço et al. (2018a) also found a Nitrosovibrio spp. RY3C‐like OTU as a probable responsible for N2O emissions, along with OTUs similar to Nitrosospira multiformis and Nitrosospira spp. PJA1, using the same FUNGENE amoA AOB database used here. The Nitrosospira in general are widespread spiral soil bacteria with generally low specificity for ammonia and, thus, found in soils under high levels of ammonia (Jia and Conrad, 2009; Di et al., 2010a; Sterngren et al., 2015). The other AOB identified here was Nitrosomonas, which was present in low abundance in the soils and was not linked to N2O emissions. The Nitrosomonas are also generally found in soils with high N inputs; moreover, Nitrosomonas europaea has a 3.5‐fold higher Vmax compared to Nitrosospira sp., suggesting that these AOB might compete with Nitrosospira in soils regularly fertilized with N (Taylor and Bottomley, 2006). Given our results, we suggest that the conditions of tropical soils used in the present study (generally low N with occasional high N inputs from fertilization) selects for the Nitrosospira, and perhaps a Nitrosomonas species adapted to these conditions but without contributing to N2O emissions was also present. Further work can focus on culturing the Nitrosospira spp RY3C‐like nitrifier identified here to verify their N2O‐production and apparent inhibition by DCD and DMPP.

The AOB are widely implicated in N2O emissions under conditions favouring nitrification in tropical and temperate soils, in contrast to the AOA (Di et al., 2010b; Liu et al., 2016; Hink et al., 2017; Theodorakopoulos et al., 2017; Meinhardt et al., 2018). This is thought to be linked to the enzymatic capabilities of different AOB and AOA species, with the former generating higher amounts of N2O through both abiotic (nitric oxide oxidation by O2) and biotic (incomplete hydroxylamine oxidation and nitrifier denitrification) mechanisms, while the latter likely emits lower N2O using only an abiotic (nitric oxidation by O2) mechanism (Harper et al., 2015; Kozlowski et al., 2016). While the AOB Nitrosospira was abundant in the soils under urea and urea with DCD treatments, we found that in the unfertilized and in the urea with DMPP treatment, the AOA Nitrososphaera were more abundant. More than five AOA Nitrososphaera 16S rRNA OTUs were identified compared to the two AOB Nitrosospira 16S rRNA OTUs; this supports the idea that the conditions in these unfertilized soils normally support the AOA Nitrososphaera rather than the AOB Nitrosospira or Nitrosomonas as the main ammonia oxidizers. Moreover, these native Nitrososphaera appeared to be non‐N2O‐producing AOA. These results support observations that the AOA Nitrososphaera is associated with low concentrations of ammonia linked to the stronger affinity of the archaeal ammonia monooxygenase for ammonia (Sterngren et al., 2015).

Interestingly, we identified two types of Nitrososphaera (AOA): one cluster of Nitrososphaera OTUs was more abundant in the soils with urea and DMPP, while the other cluster was more abundant in the unfertilized soils and co‐varied with the NOB Nitrospira. The Nitrospira was the only nitrite‐oxidizer found in our soils according to the 16S rRNA gene sequence data; interestingly, this was most abundant in the unfertilized soils and co‐varied with AOA Nitrososphaera OTUs. The Nitrospira are thought to be adapted to low NO2 − availability (Nowka et al., 2015), which might explain their presence in our soils instead of Nitrobacter (Attard et al., 2010; Gruber‐Dorninger et al., 2015). Further, perhaps the Nitrososphaera and Nitrospira naturally interact in these unfertilized soils, as has been suggested for unfertilized grassland soils and at the field level in agricultural soils (Stempfhuber et al., 2016). Future work could focus on this hypothesized interaction between non‐N2O‐generating Nitrososphaera and Nitrobacter, which appears to be selected for by low levels of available substrate and might be enhanced by adding organic residues with high C:N (Levičnik‐Höfferle et al., 2012; Hink et al., 2017).

The inhibitors DCD and DMPP are both thought to inhibit ammonia monooxygenase by chelating the Cu cofactor in the enzyme (Ruser and Schulz, 2015). The limitation of Nitrosospira but not Nitrososphaera by DCD has been shown also in a paddy field soil and in microcosms of Nitrosospira multiformis but not Nitrososphaera viennensis (Shen et al., 2013; Fu et al., 2018). Based on gene copy numbers, the AOB but not the AOA were inhibited by DMPP in a sandy soil (Duncan et al., 2017); and the AOB but not the AOA were inhibited by DCD in a grazed grassland system (Di et al., 2010b). In a Chinese vegetable soil, DMPP rather than DCD was revealed to be the more effective inhibitor of N2O‐producing AOB rather than AOA, although the N source urea was also amended with manure (Kou et al., 2015). In studies of nitrification in agricultural soils, DMPP inhibited AOB expression under neutral pH conditions (Shi et al., 2016, 2017). The different success of the nitrification inhibitors appears to be a function of temperature, Cu‐levels, and variation in abundance, genetic potential and/or expression levels of the targeted nitrifiers (Ruser and Schulz, 2015). The different effects of DCD and DMPP on the abundance of the AOB Nitrosospira and the AOA Nitrososphaera found here suggests that evaluating the nitrification dynamics of these species in culture would be interesting for future work.

In contrast to our hypothesis that the nitrification inhibitor treatments would decrease the amoA AOB community alpha diversity, this diversity remained largely unchanged across treatments. There overall was low alpha diversity of the amoA AOB community, which was supported in both the amoA AOB and 16S rRNA sequence results. Nitrifiers occupy a specific functional niche in the soil environment, and the nitrifying functions are restricted to a handful of genera; new AOB are not likely to appear at least over the relatively short duration of this experiment (in total 258 days, subset presented here was 158 days). Moreover, the sugarcane plant competes with microbes for NH4 + and NO3 − and these substrates are not likely to remain immobile long in this soil (Hajari et al., 2014). The highly weathered soils have high soil drainage capacity and have been under more than 20 years of sugarcane cultivation. Due to the long time of cultivation by sugarcane, likely the nitrifiers found in this soil are those that are adapted to the natural unfertilized conditions, to the brief high inputs of ammonia through urea fertilization, and to the competition with the sugarcane plant for ammonia. We speculate that the overall low nitrifier diversity and the selection of the nitrifiers that are present in these soils are driven by the generally low N levels.

Caveats to our methods should be mentioned, as well as how we minimized biases inherent in amplicon metagenomics (Weiss et al., 2017). To reduce noisy OTUs in both the amoA AOB and 16S rRNA datasets, we used the data to guide cutoffs for OTU inclusion. Because of the low coverage and diversity of the amoA AOB dataset, we turned to the previously generated 16S rRNA amplicon dataset, first ensuring that the same genera identified in the amoA dataset were present past the cutoffs. The previous qPCR data was generated using independent, duplicate reactions, improving our confidence in the accuracy of these copy number estimates (Soares et al., 2016). Further, the precision of the OTU classification was dependent upon the coverage of the databases used; for example, for our 16S rRNA dataset we were only able to confidently classify to the genus level. This prevented us from directly comparing the classification results between the amoA AOB and 16S rRNA datasets at the species level. However, the congruence of the amoA and 16S rRNA sequence data relative to the absolute abundances of the amoA AOB at genus level was convincing to us. Further, the low diversity of the amoA bacterial communities was echoed in the 16S rRNA data with just a few OTUs identified as Nitrosospira and only one as Nitrosomonas. Last, though the 16S rRNA samples had high Good's coverage values between 85% and 99%, there is a possibility that the nitrifying subset of the community did not have such high coverage values. However, the focus of this study was the amoA AOB nitrifiers, although future studies could target in more depth and with more specificity the nitrifying network in these soils.

In summary, the nitrification inhibitors in our experiment were revealed to target the N2O‐producing bacterial ammonia‐oxidizer Nitrosospira spp. in the soils. The low N availability appeared to drive the nitrifier community found in these soils, which should be explored in future studies. Treatment with urea and DMPP appeared to favour one functional type of AOA Nitrososphaera while the unfertilized soils revealed potentially interdependent, non‐N2O‐producing AOA Nitrososphaera and NOB Nitrospira. Our results support the use of DMPP and especially DCD as inhibitors of N2O‐producing Nitrosospira spp. in tropical soils under sugarcane. The DMPP treatment may also increase the amount of NH4 + in the soil, allowing the sugarcane crop to uptake this N source while blocking the N2O from Nitrosospira. Furthermore, we provide evidence that the nitrification process in these soils is controlled by a few bacterial and archaeal species, driven mainly by the overall low N levels and which have contrasting functional potentials for N2O emission rates.

Experimental procedures

Experimental design and sampling summary from the original experiment

A field experiment on Typic Hapludox soil (also known as Red Latosol) was set up at the Agronomic Institute in Campinas, Brazil at 22°52′15″ S, 47°04′57″ W, as described previously (Soares et al., 2015, 2016). Briefly, four treatments containing four replicate plots each were established in the 2013/2014 season on a third cycle of ratoon sugarcane (cultivar SP791011). The treatments were (1) no N fertilizer (control), (2) urea (UR), (3) urea with dicyandiamide (UR + DCD), (4) urea with 3,4‐dimethylpyrazone phosphate (UR + DMPP). Urea was incorporated into the first 5 cm of soil and applied at a rate of 120 kg N ha−1. The DCD (Sigma Aldritch) was added at 5% DCD‐N per unit N from urea (v/v) while powdered DMPP (Sigma Aldritch) was added as 1% DMPP (w/w). Gas emission rates of CO2, CH4 and N2O were measured daily to monthly using static chambers fixed 5 cm deep in between sugarcane rows (Soares et al., 2016). Gas samples were taken at three‐time intervals: 1, 15 and 30 min; the gas samples were stored in pre‐evacuated Extainers® vials (Labco Limited, Ceredigion, UK) and analysed in a Shimadzu gas chromatograph (GC‐2014). In summary of the previous soil sampling scheme, soil samples were taken of the top 10 cm of soil such that three subsamples were combined per plot. The soil samples were collected at eight time points: 7, 16, 18, 27, 35, 42, 82 and 158 days following fertilizer application and stored at −80 °C. Last, total DNA was extracted from the composite soil samples using a Power Soil kit from Mobio (Carlsbad, CA, USA) without modifications (Soares et al., 2016). One plot was dropped at random due to cost constraints, resulting in a total of 96 DNA samples (3 replicate plots × 8 timepoints × 4 treatments). Further, pH, NO3‐N and NH4‐N were measured from the soil samples and water‐filled pore space (WFPS) and temperature was previously calculated (Soares et al., 2016). In addition to the previously generated dataset of 16S rRNA sequences, the previously generated qPCR dataset of 16S rRNA, nirK, nirS, nosZ and amoA (AOB) and amoA (AOA) gene abundances, which were made from the same DNA samples as the 16S rRNA and amoA AOB (described below), were used for the current study (Soares et al., 2016).

To investigate the dynamics of ammonia‐oxidizing bacterial abundances over the experiment, we sequenced and annotated the amoA AOB amplicons resulting in a dataset containing the relative abundances of ammonia‐oxidizing bacterial species. Because of the challenges inherent in interpreting compositional data such as relative abundances – for example, an increase in relative abundance does not necessarily point to an increase in cell density (Props et al., 2016Weiss et al., 2017) we attempted three transformations of this dataset, resulting in three versions: relative abundances, variance‐stabilized abundances and absolute abundance estimates. Absolute abundance estimates were calculated using the previously generated dataset of amoA AOB qPCR gene abundances (Soares et al., 2016). The variance‐stabilized transformation for the amoA AOB dataset was not possible due to low number of sequences. To support the amoA AOB analyses, we mined the previously generated 16S rRNA dataset and repeated the analyses using the 16S rRNA‐based relative abundances of nitrifiers. These abundances were analysed using the three transformations.

amoA AOB amplification and sequencing

Amplification of the partial amoA bacterial gene (491 bp) was performed on the 96 DNA samples previously generated (Soares et al., 2016). Here, the amoA amplification was carried out using a two‐step barcoding approach. The first PCR from the total DNA samples was carried out using forward primer H‐AmoA1F‐mod (5’‐GCTATGCGCGAGCTGCGGGGHTTYTACTGGTGGT‐3′) and reverse primer H‐amoA2R (5’‐GCTATGCGCGAGCTGCCCCCTCKGSAAAGCCTTCTTC‐3′) (Rotthauwe et al., 1997; Herbold et al., 2015). In the second PCR, the amoA amplification products were amplified with primers that consisted of a 16 bp head sequence and included at the 5′ end a library‐specific 8 bp barcode (Hamady et al., 2008). Each PCR reaction (20 μl in first step, 50 μl in second step) consisted of 0.025 units of FastStart Taq DNA Polymerase (Roche), 1x reaction buffer with MgCl2 (Roche), 0.5 mM dNTPmix (Fermentas), 0.125 μM of the forward and reverse primers, 0.1 mg/ml bovine serum albumin and 1 μl of DNA template. Thermocycler (C1000 Touch Thermal cycler, Biorad) conditions were as follows: (1) 5 min at 95 °C; 35 times 30 s at 95 °C, 30 s at 53 °C, 30 s at 72 °C; and 7 min at 72 °C and (2) 5 min at 95 °C; 10 times 30 s at 95 °C, 30 s at 53 °C, 1 min at 72 °C; and 10 min at 72 °C. The first PCR reaction was performed in duplicate, screened by gel electrophoresis and pooled for use as a template in the second step, which used one primer (5’‐BARCODE‐HEAD‐3′). Second step PCR products were checked by agarose gel electrophoresis and the concentration and quality determined using a fragment analyser (Advanced Analytical). The bar‐coded PCR products from all samples were normalized in equimolar amounts before sequencing. The amoA amplicon pool was sequenced using MiSeq V3 (2x300bp) technology (LGC, Germany). To complement the analysis of the amoA amplicon sequences, we mined the previously published dataset of 16S rRNA partial gene amplicons (Soares et al., 2016). The amoA AOB amplicons were obtained from the same total DNA samples as the 16S rRNA amplicons.

amoA AOB amplicon sequence processing

Bioinformatics steps were performed on a multi core server with 64 threads running Linux Ubuntu 16.04. Processing was accomplished through a Snakemake pipeline and bash and perl scripts. The amoA AOB sequences were clipped of primers and barcodes using bbduk (bbmap version 35.82) and the paired‐ends were merged with the ‘join_paired_ends.py’ script from ea_utils version 1.1.2–537. The AOB merged sequences were dereplicated and clustered into 97% AOB OTUs with minimum size of 2 using USEARCH version 9.2.64 (commands: derep_fulllength and cluster_otus; Edgar, 2010). These parameters were chosen based on the recommendation found in the USEARCH manual (see also Supporting Information Fig. S1). To confirm the functional potential of the OTUs as amoA (KEGG pathway K10944), the centroids were compared to the KEGG database (2014‐03‐17 version) using uproc‐dna (UPROC v1.2.0; (Meinicke, 2015)). The table of OTU abundances across samples was created with the usearch global command based on 97% identity of sequences to the OTUs. Taxonomy was assigned to OTU centroids by diamond blastx v0.8.20 against the 2016‐10‐04 NCBI‐nr database (Buchfink et al., 2015). When this step yielded only classifications in the category ‘environmental samples’, taxonomy was assigned instead by best blastn (e‐value cutoff 0.02; blast v2.6.0) comparison against the custom amoA database described below. The raw sequences were submitted to the European Nucleotide Archive (ENA) under accession numbers ERS3128792 to ERS3129046.

To support the taxonomic classification results, a phylogenetic tree was created to depict the relationships between the 54 amoA OTUs and their closest matches in the custom amoA database. The latter was constructed as follows. High‐quality amoA AOB sequences were downloaded from the FUNGENE RDP database (v9.4.1) with score above 350, HMM coverage above 80% and a minimum amino acid size of 270. Duplicates were reduced to one entry. The amoA OTU centroids and reference amoA AOB sequences along with an outgroup amoA sequence from Nitrosococcus oceani C‐27 were aligned using ClustalW and used as input to make a phylogenetic tree in MEGA7 (Saitou and Nei, 1987; Kumar et al., 2016). The Maximum Composite Likelihood method was used to calculate phylogenetic distances, and bootstrap tests with 1000 replicates were performed (Felsenstein, 1985). The iTOL was used to create the final tree with bootstrap values of at least 90% depicted on the branches (Letunic and Bork, 2016).

amoA AOB OTU processing and beta and alpha diversity analyses

Statistical analyses were carried out in R version 3.3.1 using R‐Studio version 1.0.136. The R package phyloseq was used to handle the amplicon datasets. To remove undersequenced samples, the amoA AOB samples with less than 120 sequences were filtered out. To evaluate the sequencing coverage of the AOB communities, Good's coverage was calculated (package jfq3/QsRutils) and rarefaction curves were produced. Three normalization methods were attempted as follows: raw abundances were converted to percentages of total sample abundances (relative abundances). The variance‐stabilizing transformation of DeSeq2 was applied to the raw abundances (variance‐stabilized abundances; Love et al., 2014), but this was not successful due to low number of sequences. Last, relative abundances of the amoA AOB OTUs were converted to absolute abundance estimates by multiplying by sample the relative abundances by the relevant gene copy numbers previously obtained (Soares et al., 2016, as in Datta et al., 2016).

To ascertain the effect of treatment on the AOB community structure, we ordinated the amoA AOB samples using Bray–Curtis distances based on the relative abundances. Multivariate homogeneity of dispersion was checked with function ‘betadisp’ in the vegan R package. If dispersions were homogeneous, the effects of time point, treatment within time point and time point within treatment were assessed through PERMANOVA analyses (‘vegan’ R package). Post hoc tests of different pairwise group means were carried out using the ‘pairwiseAdonis’ R package (Arbizu, 2017).

To determine the effect of treatment and time point, treatment within time point and time point within treatment on the AOB community alpha diversity, the raw abundances were first rarified to 120 sequences across samples using random seed 42. After confirming that all the data were not normal using the Shapiro–Wilk test and visual check of quantile plots, two‐way crossed analyses of treatment and time point and one‐way analyses of treatment within time point and time point within treatment were evaluated using Kruskal–Wallis tests. These were supplemented with Dunn's post hoc tests.

16S rRNA OTU processing and differential abundance and indicator species analyses

We supplemented the analyses of the amoA AOB dataset using the previously published 16S rRNA gene sequence dataset (Soares et al., 2016). Good's coverage was calculated and rarefaction curves were produced as described for the amoA AOB OTU dataset. The 16S rRNA OTU abundance dataset was processed as follows. Samples with less than 3000 sequences and 16S rRNA OTUs with less than 23 sequences across all samples were filtered out. The same three transformations were applied to the 16S rRNA data as described above for the amoA AOB data. To determine significantly different nitrifier 16S rRNA OTUs between treatments, differential abundance analysis was applied between treatment pairs considering all time points. The DeSeq2 package, which applies a negative binomial and total sum scaling transformations of the filtered abundance data to stabilize variances, was used for the differential abundance testing (Love et al., 2014). The Wald test with local model fit was applied; orthogonal contrasts of the control and all other treatments, and of the urea against the treatments with a nitrification inhibitor, were carried out using Bonferroni–Hochberg correction for multiple tests. Significantly different 16S rRNA OTUs with Bonferroni‐adjusted p values of less than 0.05 were identified.

In parallel to the differential abundance testing using the variance‐stabilized abundances, the 16S rRNA relative abundances were converted to absolute abundance estimates using the 16S rRNA copy numbers previously obtained by real‐time PCR (Soares et al., 2016). To examine the 16S rRNA OTUs that were potential indicators of combinations of up to three treatments, we used the multipatt function from the ‘indicspec’ R package to apply Legendre's indicator species analysis on the 16S rRNA absolute abundances. Multiple comparison p values were adjusted using the Benjamini–Hochberg correction.

Spearman correlations of amoA AOB and nitrifier 16S rRNA OTUs with environmental variables

Correlations between nitrifier abundances and environmental variables of published environmental data were revealed using Spearman correlations (Soares et al., 2016). Log transformations of the gene copy numbers obtained by qPCR (nirS, nirK, amoA AOB, amoA AOA, total Archaeal, 16S rRNA) were carried out leaving the other variables (CO2, N2O, CH4, soil NH4–N, soil NO3–N, soil pH and WFPS) untransformed (Supporting Information Fig. S2). The nitrifier 16S rRNA and amoA AOB relative abundances and absolute abundance estimates, and the nitrifier 16S rRNA variance‐stabilized abundances, were independently correlated with the environmental variables. Significant correlations (p < 0.01) were kept; for visualization the correlations were clustered using complete linkage clustering through the ‘corrplot’ package.

Conflict of interest

The authors declare no conflict of interest.

Authors’ contributions

EEK and NAC designed the current study. JRS and HC designed the original experiment. JRS and KSL carried out the field work and DNA extractions. JRS and AP ran the qPCRs. AP prepared the amoA amplicon library for sequencing. NAC performed all bioinformatics processing and NAC and JRS performed the statistical analyses. NAC wrote the paper. JAV, HC, KSL and EEK contributed to the interpretation of the results and discussion of the article. All authors read and approved the final version of the manuscript.

Appendix S1: Supporting Information

Click here for additional data file.