Long-term N-addition alters the community structure of functionally important N-cycling soil microorganisms across global grasslands

,


Introduction
Grasslands are globally important ecosystems that cover up to 25% of the terrestrial land surface and help to support more than 1.3 billion people through livestock production (Lieth, 1978;Rogiers et al., 2010). The productivity of natural grasslands depends on the level of nitrogen (N) inputs (Lebauer and Treseder, 2008;Fay et al., 2015). This is largely determined by biological N 2 fixation as a natural supply of N to grasslands which can reduce the use of N fertilizers. N 2 fixation stems from soil diazotrophs (N 2 -fixing bacteria and archaea) that convert dinitrogen (N2) gas to ammonia (NH3) using the nitrogenase protein complex (Nif) (Reed et al., 2011;Vitousek et al., 2013;Van Groenigen et al., 2015). Diazotrophs are common in soils but their abundance, diversity and community structure are still poorly explored in grasslands in spite of their key role in global N input (Patra et al., 2006;Regan et al., 2017;Che et al., 2018). Variations in the diversity and composition of soil diazotrophs are strongly associated with vegetation type, climatic and soil parameters, whereas soil pH and moisture were the main influencing parameters (Collavino et al., 2014;Tu et al., 2016).
The N fixed by diazotrophs is directly available to plants (Vitousek et al., 2013). Nitrification is the biological oxidation of ammonia to nitrite followed by the oxidation of the nitrite to nitrate occurring through specialized microorganisms (Van Groenigen et al., 2015). Ammonia oxidation is a particularly important N-cycle process because it transforms the relatively immobile ammonium (NH 4 + ) to mobile nitrate (NO 3 − ) (Di et al., 2009;Di and Cameron, 2016), which is then prone to loss via leaching and denitrification (N 2 O emission) (Liu et al., 2018;Tao et al., 2018). Overall, grazed grasslands contribute roughly 28% of global N 2 O fluxes to the atmosphere (Rafique et al., 2011;Du et al., 2021). Ammonia oxidation therefore plays a key role in N losses and environmental pollution in grassland ecosystems (Risch et al., 2019). Ammonia-oxidizing archaea (AOA) and bacteria (AOB) are the most important ammonia oxidizers central to the fate of N in the environment Alves et al., 2013;Aigle et al., 2019). Besides AOA and AOB, complete ammonia oxidation (comammox) can also contribute to the conversion of ammonia to nitrite/nitrate (Daims et al., 2015). Apart from soil pH, vegetation and soil texture, N availability has been widely reported as a key driver of AOA and AOB community structure Gubry-Rangin et al., 2011;Liu et al., 2018). Ammonia oxidizers (AOA and AOB) are tightly linked to the diazotrophs. Diazotrophs play a critical role in global N inputs into ecosystems, and ammonia oxidizers are essential for N turnover. The latter group of microorganisms are particularly important because they are involved in a bottleneck process in the global N-cycle (Lehtovirta- Morley, 2018). However, these three important N-cycling soil functional communities are rarely studied together. In particular, the possible ecological implications of shifts in the abundance, alpha-diversity and composition of these communities in response to increased N inputs in global grasslands are still largely unknown. For example, soil diazotrophs become less important and N 2 fixation less prevalent in environments with anthropogenic N inputs, since N 2 fixation is known to be a high energy process (Fan et al., 2019;Shi et al., 2021). This may imply that N-addition could replace N 2 fixation. On the contrary, higher N availability by N-addition could increase the abundance of ammonia oxidizers, in particular increasing the AOB activity, and consequently leading to higher N losses (Di and Cameron, 2016;Cardenas et al., 2019;Du et al., 2021). From local site grassland studies responses of soil diazotrophs and ammonia oxidizers to N-addition were inconsistent. N-addition reduced (Zhang et al., 2013;Hu et al., 2021: N added as urea) or increased the abundance of soil diazotrophs (Ning et al., 2015: ammonium nitrate). N-addition as urea increased the abundance of AOA and AOB in grassland Loess Plateau soils , whereas urea addition in an alpine grassland soil of the Qinghai-Tibetan Plateau only increased the AOB abundance but AOA abundance remained unchanged (Xiang et al., 2017). In addition, most of these studies have focused on high-yield grasslands with high mean annual precipitation. Hence, these results limit our ability to extrapolate how N addition affects N-cycling microorganisms and underlines the need for global studies across a large range of climatic and edaphic factors. Recent global surveys on total microbial communities in grassland soils revealed strong regional scale patterns and nutrient addition effects at short-term , but global studies focusing on soil diazotroph and ammonia oxidizer responses to N-addition remain limited, especially longer-term responses.
Here, we sampled soils from 30 natural and semi-natural grasslands that span a globally relevant range of mean annual temperature and mean annual precipitation. We examined the impact of repeated Naddition (added as urea) on the three functionally important soil Ncycling microbial communities (soil diazotrophs, AOA, and AOB). We added granular, non-sulfur coated urea (NH 2 ) 2 CO because it is a slowrelease N fertilizer and causes less soil acidification (Cai et al., 2014). We also explored the main biotic and abiotic factors driving the geographic distribution of the functional gene abundance, observed richness, Shannon diversity index and community structure of soil diazotrophs, AOA, and AOB. This study is part of a larger effort within the Nutrient Network Global Research Cooperative to study soil microbial communities in global grasslands Prober et al., 2015;Risch et al., 2019;Lekberg et al., 2021;Nepel et al., 2022). This cross-continent experiment allows the assessment of guild-level responses across treatments and global distributions. The advantage and power of this global study lies in the standardized treatments and exact replication across sites around the world. We are not aware of any previous studies in which diazotrophic and ammonia-oxidizing communities in grasslands were concomitantly assessed across diverse grassland communities. We hypothesized (H1) that long-term N-addition would reduce the abundance of soil diazotrophs, because N-addition can replace N 2 fixation (Fan et al., 2019). Further, we expected (H2) that AOA abundance would decrease with N-addition, because AOA thrive in nutrient-poor soils with low NH 4 + concentrations and that AOB abundance would increase because of enhanced NH 4 + availability leading to a decrease in the AOA:AOB ratio (Gubry-Rangin et al., 2011). We also hypothesized (H3) that alpha-diversity (observed richness and Shannon diversity index) and community structure of AOB will respond more strongly to N-addition than the diazotrophic and AOA communities because AOB perform better in N-rich environments (Sterngren et al., 2015;Trivedi et al., 2019). These findings will allow us to better understand the responses of microbial-mediated N-cycling processes in global grassland ecosystems to N-fertilization and atmospheric N deposition.

Study sites and experimental design
Thirty sites from the Nutrient Network Global Research Cooperative (https://nutnet.umn.edu: https://nutnet.umn.edu) were used in this study (Fig. 1, Tables S1 and S2). Thirteen sites overlapped the ones used in previous studies conducted earlier and represent a different soil sampling campagin Prober et al., 2015). Most of the 30 sites were located in North America, Europe and Australia. A total of five sites were located in South America, Asia and Africa, where belowground diversity and processes are under-studied (Guerra et al., 2020). Sites were subjected to 1-9 years of urea fertilization (Table S2), building on the previous studies of 2-4 years of fertilization Prober et al., 2015). Our study therefore sheds light on impacts of N-addition across a larger global distribution and range of climatic and edaphic conditions, including mean annual temperature from − 4.1 • C to 22.3 • C, mean annual precipitation from 252 mm to 1592 mm, and elevation from 6 m to 4261 m above sea level (Fig. 1, Table S1). Soil organic carbon (C) varies from 0.3% to 22%, soil total N from <0.1% to 1.3% and soil C:N ratio from 9.1 to 23.6 across our 30 sites (Table S2). Soil clay content (3.0%-53.3%) and soil pH (3.4-7.7) also spanned a large range across the 30 sites ( Fig. 1, Table S2). Climate data for each site were obtained from WorldClim (http://www.worldclim.org/curr ent, Table S1; Hijmans et al., 2005).
The full experimental design of the NutNet global experiment has been described by Borer et al. (2014). Fertilization treatments were applied to 5 m × 5 m plots arranged in a randomized complete design with at least 3 replicates. Individual plots were at least 1 m apart. This study focused on untreated control and N addition, where N was added at 10 g m − 2 per year (equal to 100 kg N ha − 1 per year) as time-release urea [(NH 2 ) 2 CO]. This rate corresponds to a typical fertilization regime of most grasslands but is signficantly higher than natural N deposition in grasslands. Time-release urea was chosen because it causes less soil acidification than that associated with ammonium nitrate (Cai et al., 2014). The slow-release form releases urea during rain events, better approximating natural wet N deposition (Tian and Niu, 2015).
Soils were sampled from a 1 m × 1 m sampling plot specifically designated for soil sampling (Borer et al., 2014). Each site received an identical soil sampling kit and soils were sampled as described in detail by Risch et al. (2019). Briefly, three 5 cm diameter soil samples (0-12 cm depth) per plot were collected. Two samples were composited, and one subsample of approximately 50 g was taken to analyse total C and N content, organic C content, NH 4 + and NO 3 − concentration and pH (Risch et al., 2020). Another subsample was sieved (4 mm) and used to determine the N-cycling functional microbial communities (six samples per site). The third 5 × 12 cm subsample was collected and sealed in a steel cylinder, remaining undisturbed. From this subsample water holding capacity, bulk density and soil texture were assessed. All soil samples were collected 6-12 weeks prior to peak aboveground biomass at each site and were immediately shipped to the Swiss Federal Institute for Forest, Snow and Landscape Research WSL in Birmensdorf, Switzerland for laboratory analyses.

Laboratory procedures
The soil samples for chemical and particle size analyses were dried at 60 • C. Soil particles were fractionated into sand, silt and clay using the pipette method according to Gee & Bauder (1986). Mineral soil pH was measured potentiometrically in 0.01 M CaCl 2 (Risch et al., 2019). A subsample of the dried samples was ground, 40 mg of soil from each sample was weighed into a tin capsule, and total carbon (TC) and total nitrogen (TN) contents were measured with a CHN analyser (NC-2500; CE Instruments, Wigan, United Kingdom). The organic C (Corg) content of samples with a pH (CaCl 2 ) < 6.0 was assumed to be equal to the total carbon (Ctot) content, whereas samples with a pH > 6.0 were assumed to potentially contain carbonates and, thus, were fumigated with HCl vapor prior to the CN analysis (Risch et al., 2019). In addition, an equivalent of 20 g dry soil was extracted with 1 M KCl (1:5 m/w) and NO 3 − and NH 4 + concentrations were measured. Ammonium (NH 4 + ) was measured photometrically with a FIAS 300 flow injection system (Per-kinElmer, Waltham, MA, United States). The NO 3 − concentration was determined using UV-absorption but with CuSO 4 -coated zinc granules as a reducing agent (Risch et al., 2019). Plant species richness and cover (in %) of each plant species were assessed each year in a 1 × 1 m subplot at each site (Borer et al., 2014).

Soil DNA extraction and quantitative PCR
Soil DNA was extracted from 0.30 g of fresh soil using a PowerSoil™ DNA Isolation Kit (MO BIO Laboratories, Carlsbad, CA, USA). To remove foreign DNA and prevent microbial contaminations, working bench surfaces and non-autoclavable materials were cleaned with 5% sodiumhypochlorite and 70% ethanol solutions prior to the DNA  (Table S1). Our sites also cover a wide range of soil edaphic conditions, as described in the main text and shown in Table S2. Numbers refer to # in Tables S1 and S2. extractions. The extracted DNA was quantified with Quant-it Pico Green dsDNA Assay (Thermo Fisher Scientific, Reinach, Switzerland). Soil diazotrophs and ammonia oxidizers have been extensively investigated with molecular marker methods (Levy-Booth et al., 2014). The nifH gene family, encoding the dinitrogenase reductase subunit, is widely used to gain insight into diazotrophic communities (Hsu and Buckley, 2009;Gaby and Buckley, 2012;Collavino et al., 2014). The amoA gene, encoding subunit A of ammonia monooxygenase, is commonly used to investigate AOA and AOB (Tourna et al., 2008;Nicol et al., 2008). The abundance of functional marker genes encoding archaeal and bacterial ammonia monooxygenase (archaeal amoA, bacterial amoA) and the dinitrogenase gene (nifH) were quantified using quantitative PCR (qPCR). The primers and thermocycling conditions for bacterial amoA (AmoA-1F and AmoA-2R; Rotthauwe et al., 1997), archaeal amoA (Arch-amoAF and Arch-amoAR; Francis et al., 2005) and nifH (PolF and PolR; Poly et al., 2001) are reported in Table S3. Negative controls for the DNA extractions (extraction buffer without soil) and PCR amplifications (high-purity water without DNA template) were included. PCR for each sample was performed in triplicate and negative controls were included in each batch of PCR.
All qPCR assays were carried out with equal amounts of template DNA (20 ng DNA) in 20 μl reactions using QuantiTect SYBR Green PCR master mix (Qiagen, Hirlen, Germany) on an ABI7500 Fast Real-Time PCR system (Applied Biosystems, Foster City, CA, USA) as previously outlined (Frey et al., 2011(Frey et al., , 2021. Standards for qPCR were generated by serial dilution of stocks containing a known number of plasmids carrying the respective functional gene as an insert (Frey et al., 2020). The specificity of the amplification products was confirmed by melting-curve analysis, and the expected sizes of the amplified fragments were checked in a 1.5% agarose gel stained with ethidium bromide. The amplification efficiencies of qPCR assays were 95% (±2) for archaeal amoA, 93% (±3) for bacterial amoA and 98% (±2) for nifH. R 2 values were >0.99 for all runs. The gene copy number was calculated using a regression equation that related the cycle threshold (Ct) value to the known number of copies in the standards. For each soil sample, the qPCR reactions were repeated three times. Abundances of bacterial and archaeal amoA and nifH genes refer to copy numbers per gram of dry soil.

MiSeq sequencing and bioinformatic analyses
PCR amplifications of the bacterial amoA and nifH genes were performed as described previously with primers similar to those used in qPCRs but containing adapters for MiSeq Illumina sequencing. Due to the limits of sequence lengths using the MiSeq platform (2 × 300 bp), primer amplification fragments of archaeal amoA were too long (629 bp) to be assembled; therefore, another primer set was used (Arch-amoA-104F and Arch-amoA-616R; Tourna et al., 2008;Alves et al., 2013) which generated a shorter amplicon (Table S3). Each sample consisting of 40 ng DNA was amplified in triplicate and pooled before purification with Agencourt AMPure XP beads (Beckman Colter, Berea, CA, USA) and quantified with the Qubit 2.0 fluorometric system (Life Technologies, Paisley, UK). Amplicons were sent to the Genome Quebec Innovation Centre (Montreal, Canada) for barcoding using the Fluidigm Access Array technology and paired-end V3 MiSeq sequencing (Illumina Inc., San Diego, CA, USA), enabling production of 2 × 300-bp reads.
Raw nifH and bacterial and archaeal amoA reads were each merged using PEAR v.1.9.4  with a minimum overlap of 10 bp between forward and reverse paired end reads. Quality filtering was performed on the paired reads using USEARCH v.9.2.64 (Edgar, 2010) removing sequences with a Phred score below Q30, exceeding a per read maximum expected error (maxEE) of 0.5 or containing ambiguous bases. Primers were subsequently trimmed off using CUTADAPT v.1.9 (Martin, 2011) set with an error tolerance of 0.15 and 0.20 for the forward and reverse primers, respectively. Next, possible frameshifts were detected and corrected to frame 1 using an inhouse Python script. Any reads containing stop codons (TAG, TAA, TGA), indicating erroneous indels, were removed. Reads were then sorted according to abundance after which putative chimeras were removed from the dataset using the UCHIME algorithm (Edgar et al., 2011). High quality sequences were assigned to operational taxonomic units (OTUs) using the UCLUST algorithm (Edgar, 2010) at a 5% nucleotide dissimilarity for amoA and nifH gene sequences (Gaby and Buckley, 2011;Pester et al., 2012;Penton et al., 2016;Carrell et al., 2019), with -minsize set to 2 for the nifH gene. Prior to diversity analyses, OTUs were rarefied to 11025 (nifH), 7165 (AOA amoA) and 9427 (AOB amoA) reads per sample. OTU-based alpha-diversity was reported as observed richness and Shannon diversity index, which includes community richness and evenness in the metric. Differences between community structure (beta-diversity) were estimated based on the Bray-Curtis distance. Raw sequences were uploaded in the NCBI Sequence Read Archive (SRA) with BioProject accession number PRJNA739076.

Statistical analyses
Some of the explanatory variables (Table S4) were highly skewed and were thus log or square-root transformed prior to analyses. All continuous explanatory variables were centred and scaled to have a mean of zero and variance of one. Variables were additionally filtered to avoid collinearity using correlation analysis. If variables were strongly correlated (Pearson's |r| > 0.70) (Dormann et al., 2013), those that allowed us to minimize the number of variables were selected (Fig. S1, Table S4). Consequently, soil total C concentration, total N concentration, sand content and bulk density were removed from the dataset. For climatic variables mean annual temperature, mean annual precipitation, temperature of the wettest quarter and temperature seasonality (standard deviation; hereafter temperature variability; Hijmans et al., 2005) were used, as these variables were not correlated under the definition given above (Fig. S1, Table S4).
To assess N-addition we fit linear mixed-effects models (LMMs). observed richness, Shannon diversity index, gene abundance (gene copy numbers) of all three groups (nifH, AOA, AOB), and the ratio of AOA to AOB as response variables using the 'lme' function in the nlme package (version 3.131.1; Pinheiro et al., 2021) in R version 3.4.4 (R Core Team, 2020). Treatment (N addition, control) was included as a fixed factor and block nested within site was included as a random factor. The number of years since the N-addition treatment began was also included as a fixed effect in all the initial models. However, with the exception of AOA Shannon diversity index and AOB richness, this factor was not significant and therefore was not retained in the models. To examine global biotic, edaphic and climatic drivers (Table S4) of nifH, AOA, and AOB alpha-diversity (observed richness and Shannon diversity index), gene abundance and the AOA:AOB ratio in grassland soils, multi-model inference and LMMs were used. Separate models were run for each response variable. The explanatory variables together with treatment were included as fixed effects. In addition, significant interaction terms between the treatment and explanatory variables were included. Block nested within site was again included as a random effect in these models. First, full models including all predictor variables and relevant interaction terms were fitted and then the R package MuMin (version 1.42.1; cran.r-project.org) was used to select the best models that explained the most variation based on Akaike's information criterion (AIC; 'model. avg' function). The corrected AIC (AICc) was used to account for the small sample size. Model averages are presented in the main text and all best models (models that fell within 2 AICc units of the best model; delta AICc <2) are given in Tables S6 and S7. To assess if the significant explanatory variables contained in the top models had opposite trends of N additions on alpha-diversity (observed richness and Shannon diversity index) and gene abundance of nifH, AOA and AOB as well as the AOA to AOB ratio, e.g., if the response to N additions differed between dry and wet sites, or at sites with high versus low pH we calculated the log response ratio (LRR) for gene abundance and alpha-diversity of nifH, AOA and AOB as well as for the AOA:AOB ratio. The LRR was then plotted against the explanatory variables of the control plots. To assess whether these relationships were significant, we again used LMMs.
PERMANOVA analyses (10 5 permutations) were conducted to test for significant differences in diazotrophic (nifH), AOA and AOB community structure between control and N-added plots, sites and their interaction. Site was included as a random factor using the statistical package of PrimerE v7 (Clarke and Gorley, 2015). Beta-diversity was visualized using canonical analysis of principal coordinates (CAP) based on Bray-Curtis dissimilarity distances. To explore environmental correlations with the beta-diversity of diazotrophic and ammonia-oxidizing communities (AOA and AOB), pairwise Bray-Curtis dissimilarities between all plots were calculated using vegan. Dissimilarities were averaged per treatment level within site-pairs to avoid pseudo-replication . Final correlations between environmental factors and beta-diversity were thus comparisons between sites at the treatment level, averaged from plot-level dissimilarities. Multi-model inference and linear models were then used to assess environmental correlates of between-site dissimilarities in diazotrophic and ammonia-oxidizing communities (AOA and AOB). The model selection procedure was identical to that described above for alpha-diversity, gene abundance and AOA:AOB ratio.

N-addition effects, and drivers of gene abundance and alphadiversity
AOB alpha-diversity (observed richness and Shannon diversity index) and abundance were significantly higher and the AOA:AOB ratio was significantly lower in N-added plots than in control plots (Fig. 2, Table S5, Figs. S2-S6 in the Supplementary Results). N-addition did influence nifH richness but not nifH Shannon diversity index and abundance. AOA abundance and alpha-diversity were not affected by Naddition (Fig. 2, Table S5). We detected higher nifH alpha-diversity in Nadded plots at sites with higher MAP, while larger nifH gene copy numbers were correlated with higher soil organic C content, higher temperature of the wettest quarter, lower temperature variability and higher MAP across the sites (model average results are given in Table 1, whereas climatic, soil physico-chemical and vegetation parameters included in all best models [within an AICc of 2] are given in Table S6).
In contrast, we were not able to explain AOA Shannon diversity index across our sites with the climatic, soil or plant variables available, although higher AOA richness was related to higher soil pH and soil organic C content. Similarly, higher abundances of AOA were correlated to higher soil pH (Table 1, Table S6). AOB Shannon diversity index and AOB abundance were both higher in N-added plots than in controls (Fig. 2). Besides of the N-addition treatment, AOB abundance was further explained by multiple factors (i.e., interaction between N-addition treatment and MAT; interaction between N-addition treatment and temperature of the wettest quarter; and NO 3 − content), whereas AOB Shannon diversity index only by the treatment (Table 1, Table S6). AOB richness was explained by soil pH and NO 3 − content. The N-addition treatment was the only significant predictor for differences in the AOA: AOB ratio across the global grasslands (Fig. 2, Table 1, Table S6). There were no clear trends in the LRR ratios of gene abundance and diversity of nifH, AOA and AOB as well as the AOA to AOB ratio in relationship to the significant variables included in the models (Figs. S7-S10).

N-addition effects and drivers of community structure
PERMANOVA analysis revealed a strong N-addition effect on diazotrophic (nifH) and AOB community structure (Table 2, Fig. 3). AOB communities additionally differed significantly among the 30 sites and also showed a significant treatment × site interaction (Table 2). AOA communities were, in contrast, not affected by N-addition. Betadiversity of the N-cycling communities (dissimilarity among sites) was correlated with both climatic and soil variables, while vegetation properties did not play a large role (Table 3, Table 4, Table S7). Differences between sites in terms of MAP, geographic distance, soil pH, clay, MAT, C:N ratio and temperature of the wettest quarter explained the dissimilarities in diazotrophic communities across the global grasslands. Strongest predictors in the dissimilarities in diazotrophic communities were MAP and geographic distance whereas C:N ratio and temperature were weakest (Table 3). The beta-diversity of AOA communities was related to differences among sites in MAT and in soil pH, water holding capacity, clay content and C:N ratio (Table 3, Table S7). The strongest predictors of AOB beta-diversity were differences in soil pH, MAT and geographic distance, whereas N-addition, the interaction between N-addition and the difference in C:N ratio, MAP and NO 3 − concentration had less, but still significant, explanatory power (Table 3). Fig. 2. N-addition effects on diazotrophic (nifH), AOA (archaeal amoA) and AOB (bacterial amoA) alpha-diversity (richness and Shannon diversity index) and abundance, as well as, AOA:AOB ratio across global grasslands. Results are based on linear mixed-effects models, with block nested in site included as a random effect. * Indicates a significant treatment effect (alpha-level = 0.05). Data points are effect sizes (N-addition) and error bars represent standard errors (SE). For statistics see Table S5.

Impact of N-addition on diazotrophs, AOA and AOB abundance
The abundance of nifH copy numbers was not affected by N-addition, which does not support our first hypothesis (H1). We expected that soil diazotrophs would become less abundant in environments with increased N inputs, because N 2 fixation is an energetically expensive process and the diazotrophs are outcompeted under high N availability by other non-fixing microorganisms (Fan et al., 2019;Shi et al., 2021). In contrast to our study, higher available N content in soils as a consequence of N-addition has been shown to negatively affect nifH gene copy numbers in other grassland soils (Zhang et al., 2013;Hu et al., 2021), as the dinitrogenase enzyme is sensitive to NH 4 + (Reed et al., 2011;Han et al., 2019). In a long-term study (40 years), urea addition was found to be detrimental for obligate N 2 fixers (Fan et al., 2019). The addition of N was shown to decrease symbiotic N 2 -fixation in grasslands, including many of the sites in this study (Vázquez et al., 2022). The negative effect of N addition on symbiotic N 2 fixation by legumes per area was caused exclusively by the negative effect of N addition on legume biomass (Tognetti et al., 2021), and not by an effect on the N 2 fixation rate per unit biomass. In addition to symbiotic N 2 fixing microorganisms, Rui et al. (2022) reported that free-living diazotrophs may also play important roles in biological N 2 fixation of grasslands. The capacity for biological N 2 fixation relies on the nitrogenase enzyme system (Reed et al., 2011). The nifH gene encodes the nitrogenase reductase subunit, which is the most widely sequenced biomarker for identifying diazotrophs (Gaby and Buckley, 2012). However, nifH gene abundance cannot differentiate between the free-living and symbiotic diazotrophs (Rui et al., 2022). Non-fixing plants (grasses and forbs) outcompeted legumes under high N availability and a decrease of N 2 fixation could be expected. Vegetation (diversity and biomass) is a crucial factor affecting soil diazotrophs, which are highly dependent on easily available carbon sources from root exudates (Reed et al., 2011;Collavino et al., 2014). It is possible that a period of one to nine years of N-addition in our current study was not long-term enough to cause such changes. The quantity, form and duration of N-addition play a crucial role in the N availability and acidification of the soils after N-addition (Tian and Niu, 2015). In some of the long-term grassland sites considered here (>7 years of treatment), urea fertilization acidified the soils (pH decrease >0.4 in burrawan.au, cbgb.us, cdcr.us, cdpt.us, valm.ch), but this decrease in pH (p = 0.004) apparently had no influence on the abundance of nifH gene copies. This is in accordance with Leff et al. (2015) who also included soils from some sites analyzed in our study. Overall, Table 1 Effects of significant environmental drivers on diazotrophic (nifH), AOA (archaeal amoA) and AOB (bacterial amoA) alpha-diversity (observed richness and Shannon diversity index) and abundance, and the AOA:AOB ratio across global grasslands. Model averaged coefficients for the full average are shown. For all best models for each variable see Table S6.    Effects of nitrogen (N) addition and differences in site characteristics (climatic, soil and vegetation variables) on dissimilarities of diazotrophic (nifH), AOA (archaeal amoA) and AOB (bacterial amoA) community structure (beta-diversity). Differences in site characteristics and dissimilarities in microbial communities (Bray-Curtis distance) were averaged for each treatment level within site-pairs. Model averaged coefficients (full average) are shown for those predictors that occurred in at least 50% of all models. For all best models see Table S7. Significant effects are highlighted in bold (P < 0.05). No. models = number of models where a variable was included. Rel. imp. = relative importance, Δ = difference between two sites. Definition of predictors see Table S4. changes in the soil environment were surprisingly small which might explain the lack of responses in terms of the nifH gene abundance in our study. Although we did not observe N-addition effects on nifH gene abundance, soil organic C, together with climatic properties, had the strongest positive effects on nifH gene abundance. Edaphic variables, such as organic C, NH 4 + , NO 3 − , and available P content have been widely recognized as crucial factors affecting the abundance of soil diazotrophs (Reed et al., 2011;Reardon et al., 2014;Han et al., 2019). AOA and AOB responded differently to N-addition. The gene abundance of AOB increased whereas AOA were not affected by N-addition. We expected that AOA abundance would decrease with N-addition (H2) because AOA thrive in nutrient-poor soils with low NH 4 + concentrations and that AOB abundance would increase because of enhanced NH 4 + availability. Hence, we only found partial support for our hypothesis. Yet, our findings are in agreement with other research where N-addition has been shown to positively affect (Verhamme et al., 2011) or have no influence on AOA abundance , whereas AOB abundance flourished with N enrichment (Xiang et al., 2017;Zhang et al., 2019). However, these earlier findings are mostly from local-scale studies. Dominance of AOA relative to AOB in the amoA gene pool, as we found here, has also been reported in grasslands and agricultural environments (He et al., 2007;Gubry-Rangin et al., 2011;Zhou et al., 2014;Sterngren et al., 2015). This could be attributed to a difference in the affinity for ammonia caused by the unique physiological characteristics and preferred habitat conditions of AOA and AOB (Prosser and Nicol, 2012). AOB are typically copiotrophic (Hatzenpichler, 2012) and are most abundant in soils that have a relatively high N content (Sterngren et al., 2015;Trivedi et al., 2019). AOA perform better in habitats with more acidic, organic matter-depleted soil conditions (Banning et al., 2015;Trivedi et al., 2019) which is strongly supported by our findings in some of the grassland soils (e.g. spin.us; sgs.us; kibber.in). AOA also appear to be less constrained by the availability of N substrates and prosper in oligotrophic environments (Martens-Habbena et al., 2009;Trivedi et al., 2019), and our results confirm this pattern. In addition, AOA have lower membrane permeability, smaller cell and genome size, and higher energy efficiency than AOB, which enables AOA to better adapt to and perform their physiological activities in nutrient-limited environments (Prosser and Nicol, 2012;Banning et al., 2015). Burton and Prosser (2001) suggested that AOB in contrast to AOA can produce urease and thus respond directly to urea addition which might be responsible for the stronger response of AOB compared to AOA to N-addition in our study. These and potentially other physiological or metabolic differences between AOA and AOB likely lead to niche differentiation between these two groups (Prosser and Nicol, 2012;Carey et al., 2016;Aigle et al., 2019), and therefore different responses to N-addition. Overall, our results suggest that AOB perform better in nutrient-rich environments than AOA when faced with an increased N input (H2).
As observed in our study, N-addition may change the dynamics of microbial-mediated N-cycling processes. Diazotrophs and both AOA and AOB are tightly linked (Van Groenigen et al., 2015). Therefore, shifts in the abundance of these functionally important microbial groups with the high N dose applied in this study are likely to have ecological implications in global grasslands. N-addition may disrupt the soil N-cycling by increasing the abundance of AOB. Larger population sizes of AOB may promote higher rates of ammonia oxidation (Banning et al., 2015;Carey et al., 2016), and subsequently change the availability of oxidized forms of N and thus their mobility in soils (Yang et al., 2018). In addition, higher AOB abundance may also increase N 2 O emissions, which is a potent greenhouse gas (Hink et al., 2017;Cardenas et al., 2019;Du et al., 2021). Taken together, our findings suggest that N-fertilization has important ecological implications for the abundance of functionally important microbial groups. Specially, an increased abundance of ammonia oxidizers with N-addition may increase nitrification activity, which likely leads to N losses and environmental pollution.

Impact of N-addition on diazotrophs, AOA and AOB alpha-diversity and community structure
We also hypothesized (H3) that the alpha-diversity of AOB respond more strongly to N-addition than the diazotrophic and AOA communities because AOB have a more copiotrophic lifestyle than the other two groups. Our findings partially confirmed our expectations. AOB alphadiversity (observed richness and Shannon diversity index) increased with N-addition, but N input had no significant effect on the AOA alphadiversity and 'Treatment' was not included in any of the best models we developed for AOA (Table 4). The AOB response that we observed contradicted previous work that showed that N input had no significant effect on the AOB (or AOA) alpha-diversity in a semi-arid grassland on China's Loess Plateau . Furthermore, we did find a treatment effect on nifH richness, which contrasts Zhu et al. (2018) who Table 4 Summarizing the variables for predicting diazotrophic (nifH), AOA (archaeal amoA) and AOB (bacterial amoA) abundance, alpha-diversity (richness and Shannon diversity index) and community structure, and the AOA:AOB ratio. For details about individual variables see Table S4. Tvar = temperature seasonality, T = temperature, P = precipitation, log = transformed using log, sqrt = transformed using square root. reported no effect on diazotroph richness with urea applied at the same rate as our study. However, at higher rates (150 and 300 kg urea ha − 1 per year), diazotroph richness did decrease (Zhou et al., 2021) suggesting a certain threshold of N addition levels for the nifH richness. The best predictor of nifH alpha-diversity in our study was MAP, with higher nifH diversity at sites with higher MAP, whereas soil pH and NO 3 − concentration were the best variables in explaining nifH diversity in agricultural land (Zhou et al., 2021), but we are not aware of a study that investigated these interactions in global grassland ecosystems. The AOB community structure responded strongly to N-addition whereas the AOA community structure did not (Table 4). This finding is in agreement with our third hypothesis (H3) and provides new insight into the responses of AOA and AOB community structure to N-addition in global grassland ecosystems. Earlier findings are mostly from single site studies. Xiang et al. (2017) showed a significant shift in AOB but not AOA community structure within short-term N urea amendment in an alpine grassland soil of the Qinghai-Tibetan Plateau. Similarly, Zhang et al. (2019) reported that N-fertilization (as NH 4 NO 3 at 10 g N m − 2 year − 1 for 11 years) significantly affected AOB community structure over the AOA counterparts in a semi-arid grassland system in Northern China. The most likely explanation for these findings is that most of the added N enters the soil in the form of NH 4 + , which AOB and AOA transform to NO 3 − thereby producing H + (nitrification) (Patra et al., 2006;Van Groenigen et al., 2015). Fertilization-induced changes in pH may partly explain the N-fertilization effects on AOB community structure. However, changes are more likely due to direct effects of increased N availability and interactions between AOA and AOB . Overall, these findings demonstrate that AOB were more sensitive than AOA in response to urea suggesting that N fertilizer may be an important controller of ammonia oxidizers. Besides that, AOA communities were not affected by N-addition. We could also not detect an obvious site effect on AOA community structure in our global grassland study as found for AOB. Overall, the AOA richness was low (193 OTUs; Supplementary Results) compared to the higher AOB richness (1418 OTUs) and some of the dominant AOA taxa are known to have a cosmopolitan distribution (Alves et al., 2018). Furthermore, we found that N-addition changed the diazotrophic community structure, even though the Shannon diversity index and abundance of nifH gene copy numbers were not affected (Table 4). This agrees with previous findings that diazotrophic communities are sensitive to N-addition in arable soils (Wang et al., 2017a;Zhou et al., 2021). Shifts in the nifH community structure under longer-term N-addition may lower the capacity for biological N 2 fixation, particularly in nitrogen-poor soils (Feng et al., 2018;Fan et al., 2019;Tognetti et al., 2021). This may have long-term ecological implications for soil diazotrophs and microbial-mediated N-cycling processes in grassland ecosystems.
Soil pH was one of the main predictors of the community structure of our three key microbial groups. All models included soil pH as the main influencing soil factor (Table 4). The importance of pH is consistent with previous work demonstrating that soil pH is the primary determinant for shaping total microbial community structure at the continental scale (Lauber et al., 2009;Fierer, 2017). In our study, both AOA and AOB communities were strongly affected by soil pH, which is in line with findings from single site studies (Ying et al., 2017;Zhang et al., 2019;Vázquez et al., 2020;Wang et al., 2021). Recent work by Wang et al. (2017b) and Zhou et al. (2021) also highlighted the importance of soil pH as a fundamental driver shaping the diazotrophic community structure. We reached the same conclusion here with the global grassland soils.

Conclusions
Diazotrophic and ammonia-oxidizing microorganisms play crucial roles in global N cycle. Our study provides new insight into the responses of three functionally important N-cycling microbial groups to N-addition in grassland ecosystems. Annual addition of 10 g N m − 2 year − 1 for up to 9 years did not affect nifH gene abundance or Shannon diversity but changed the nifH community structure across the global grassland sites. AOB responded more dramatically to N-addition than diazotrophs or AOA, with increases in gene copy numbers and alpha-diversity and a shift in community structure. N-addition favoured copiotrophic AOB over AOA, which were unaffected by the treatment. The generality of these patterns in a globally coordinated nutrient addition experiment allows, as far as we know for the first time, to draw conclusions that consistent, guild-level responses occur worldwide despite substantial differences in plant communities, climatic conditions, and soil properties across grasslands. Long-term N-fertilization at the dose applied in this study may have considerable ecological implications for the abundance of functionally important microbial groups and effects on nutrient cycling in general. More abundant AOB under N-addition may accelerate nitrification activity in grasslands, which has the potential to increase NO 3 − in soils, thereby leading to higher environmental pollution in the form of NO 3 − leaching or N 2 O emission. Long-term N-addition is also likely to change above-and belowground grassland biodiversity , due to soil acidification and preferential uptake of NO 3 − , and could therefore threaten the sustainability of grassland functioning.

Author contributions
AF and BF designed the study, all authors contributed to the soil sampling and collection of data. AR, BF, BM and BT analyzed the data, AR and BF drafted the manuscript to which all authors contributed.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability
Data will be made available on request.