The experiments were carried out at perennial ryegrass and white clover swards at Lincoln University’s Ashley Dene Research and Development Station, Canterbury, New Zealand (−43.6468, 172.3467). Urea was applied to separate sampling sites during the farm season (1 June 2017 to 31 May 2018) to achieve designed nitrogen treatment [no-nitrogen (0 N kg/ha/yr), medium-nitrogen (150 N kg/ha/yr), and high-nitrogen (300 N kg/ha/yr)]. The three levels of nitrogen treatment were selected based on the range of nitrogen application rates in pastoral dairy farms in Canterbury as previously described (Beukes et al., 2020).

The experimental swards were grazed by 30 lactating dairy cows (10 cows per nitrogen treatment) to allow the normal grazing practise for pasture and animals. Same amount of herbage intake was controlled across cows (∼30 kg DM/cow/day) by break feeding. To stabilise the environmental microbiomes, the experimental sites were grazed by cows 2 weeks before sampling. The allocated swards were examined daily before and after grazing using a rising-plate meter (RPM) to estimate pasture intake during sampling by using the equation RPM reading ×140 + 500 kg DM/ha.

Seven types of samples (bulk soil, cow rumen content, cow faeces, white clover leaf, white clover root, ryegrass leaf, and ryegrass root) were taken from three ecological niches (soil, plant, and animals). A total of 169 samples were used in the current study where metadata details can be found in Supplementary Table 1.

Soil cores were taken using a 110 mm wide (internal measurement) corer, with a 100 mm wide (internal measurement) plastic sleeve inserted. Plant materials were gently twisted into the centre of the core and sleeve to minimise damage during sampling. The corer was driven into the soil to a depth of 150 mm to ensure that at least the top 100 mm soil is recovered intact. Once the depth was reached, the corer was gently twisted and lifted out of the ground. Two smaller plastic sleeves were used to transfer the soil from the corer. The soil sample was then wrapped in clingfilm and stored at 4°C. Bulk soil samples were taken after removing stones and plant materials.

Rumen samples were collected by stomach tubing of animals. Faecal samples were taken with digital collection from cow anuses at the same time of rumen content collection. Both types of animal samples were flash frozen with liquid nitrogen and stored at −80°C until further progress.

For both white clover and ryegrass, individual plants were gently removed from the soil core and split into root and leaf parts with a sterile scalpel blade. After removing loose soil particles, the plant samples were then snap froze with liquid nitrogen and stored at −80°C until DNA extractions.

DNA was extracted from individual sample (300 mg/sample) with a Macherey-Nagel NucleoSpin 96 Soil kit. The manufacturer’s protocol was followed but with the following modifications. Buffer SL1 (700 μL) and enhancer SX (300 μL) was pre-combined and added to a 1.7ml tube. Grounded sample was then transferred to the tube and vortexed for 1 min. The mixture was then centrifuged for 4 min at 11,000 x g. After removing the supernatant, 150 μL buffer SL3 was added to the tube. The tube was then incubated at 4°C for 20 min and centrifuged for 2 min. The clear supernatant was then processed following manufacture protocols and eluted twice with 70 μL of buffer SE each time.

We examined microbiomes using 3 approaches focusing on specific microbiome subpopulations, namely prokaryotes (16S rRNA), eukaryotes (18S rRNA), and fungi (ITS rRNA). Next-generation sequencing data were captured from 7 farm niches, specifically 2 animal-associated niches [Rumen (R) and Faecal (F)], 3 below-ground niches [WhiteCloverRoot (WCR), RyeGrassRoot (RGR), and Bulk Soil (BS)], and 2 above-ground niches [WhiteCloverLeaf (WCL) and RyeGrassLeaf (RGL)]. For each of these microbiomes, community responses were compared across three different soil nitrogen treatments [control or no-nitrogen (0 N kg/ha/yr), medium-nitrogen (150 N kg/ha/yr), and high-nitrogen (300 N kg/ha/yr)].

Next generation sequencing was conducted with an Illumina MiSeq V2 platform at Massey University. The total community genomic DNA for each sample was sequenced targeting 16S ribosomal RNA (rRNA), 18S rRNA, and internal transcribed spacer (ITS) regions with primers following the Earth Microbiome Project.

The demultiplexed paired end sequence reads were then quality filtered and merged using dada2 package V1.12 (Callahan et al., 2016) in R V3.6.1 (R Core Team, 2019) following the Dada2 tutorial pipelines (16S and 18S, ITS) with default settings. The merged sequence reads were then denoised, chimera removed, and dereplicated into Amplicon Sequence Variant (ASV). Silva v132 (Quast et al., 2013), and UNITE databases (Nilsson et al., 2019) were used to assign taxonomy to each ASV. For each amplicon, ASV table, metadata (Supplementary Table 1), and taxonomy table were converted into a phyloseq object with phyloseq version 1.28.0 (McMurdie and Holmes, 2013). The raw amplicon (16S, 18s, and ITS) sequencing files were uploaded to the NCBI SRA database under the accession number PRJNA702151. For all sample of each dataset (16S, 18S, and ITS) sequences were rarefied to an even depth, that 10000, 25000, and 10000 reads per sample were randomly selected for 16S, 18S, and ITS datasets, respectively. For 16S phyloseq object, chloroplasts sequences were removed. Singletons were discarded prior to from downstream diversity analyses. R scripts of all data processing were accessible online from GitHub.

Amplicon sequence variant counts for each sample across farm niches were generated based on rarefied and trimmed (i.e., singletons removed) datasets (16S, 18S, and ITS). The alpha-diversity measurements across farm niches and in response to nitrogen treatments was calculated by function estimate_richness using the “Observed” method in R with packages dplyr version 0.8.3 (Wickham et al., 2019), plyr version 1.8.4 (Wickham, 2011), reshape2 version 1.4.3 (Wickham, 2007), and vegan version 2.5.6 (Jari Oksanen et al., 2019). The violin plots were created with function plot_richness facet by niche type and ecosystem type using packages ggplot2 version 1.28.0 (Wickham, 2016), Biostrings version 2.52.0 (Pagès et al., 2019), forcats version 0.4.0 (Wickham, 2019), ggpubr version 0.2.3 (Kassambara, 2019), lemon version 0.4.3 (Edwards, 2019), phyloseq version 1.28.0 (McMurdie and Holmes, 2013), RColorBrewer version 1.1.2 (Neuwirth, 2014), and scales version 1.0.0 (Wickham, 2018). Kruskal-Wallis tests were used to compare species richness across farm niches and nitrogen treatments that statistical significances were masked on the violine plots with function stat_compare_means of package ggpubr version 0.2.3 (Kassambara, 2019). Microbiome dissimilarities across all samples were determined using a phyloseq function ordinate with the Bray-Curtis distance method (i.e., distance = “bray”). Beta-diveristy of all samples across treatments were visualised with NMDS plots using function plot_ordination with stat_ellipse option. Changes in community structure and variance were tested via functions anosim and adonis, respectively, with vegan version 2.5.6 (Jari Oksanen et al., 2019). Both functions were conducted using Bray-Curtis distance matrix as the input data matrix and N treatment (i.e., no nitrogen, medium nitrogen, and high nitrogen) as grouping.

Exact tests were used to determine differentially abundant ASVs for each individual niche between nitrogen treated and untreated microbiomes illustrated with a dot plot. For each phyloseq dataset (16S, 18S, and ITS), Packages edgeR version 3.26.8 (Robinson et al., 2010), ggplot2 version 1.28.0 (Wickham, 2016), ggpubr version 0.2.3 (Kassambara, 2019), ggrepel version 0.8.1 (Slowikowski, 2019), lemon version 0.4.3 (Edwards, 2019), and phyloseq version 1.28.0 (McMurdie and Holmes, 2013) were used. Heatmap and volcano plot were used to supplement Exact tests. Heatmaps provided the presence and absence information on responsive ASVs across individual niches between samples. Volcano plots illustrated the changing direction (positive or negative) for ASVs responses across niches and between nitrogen treatments. R packages ggplot2 version 1.28.0 (Wickham, 2016), ggpubr version 0.2.3 (Kassambara, 2019), lemon version 0.4.3 (Edwards, 2019), and phyloseq version 1.28.0 (McMurdie and Holmes, 2013) were used.

Network correlations for each niche under each nitrogen treatment were constructed with packages igraph version 1.2.4.1 (Csardi and Nepusz, 2006), networkD3 version 0.4 (Allaire et al., 2017), and tidyr version 1.0.0 (Wickham and Henry, 2019) to provide insight into the microbial interactions responses to different nitrogen managements. In each niche, ASVs were chosen based on two criteria: (1) were present in at least 60% of samples; and (2) were significantly (p ≤ 0.05) differentially abundant in response to nitrogen treatments in any niche, or were the top 200 most abundant across niches. Spearman correlations were conducted between all chosen ASV pairs in each niches with package psych version 1.8.12 (Revelle, 2018). ASV pairs having strong significant correlations (r ≥ 0.8 or r ≤ −0.8, p ≤ 0.05) were included in the networks. The force-directed layout algorithm was then used to calculated the network layout. Nodes were coloured by phylum represent ASVs and edges were coloured by correlation positiveness represent significant microbial coefficients.

Niche microbiome metapopulation was investigated by analysing occupancy-frequency distribution of ASVs in each niche was made with packages ggplot2 version 1.28.0 (Wickham, 2016), vegan version 2.5.6 (Jari Oksanen et al., 2019), and scales version 1.0.0 (Wickham, 2018) following previous study (Lindh et al., 2017). In the current study, the occupancy of an ASV is defined as the fraction of samples it occupied in each niche. The frequency denotes the existence number of unique ASVs in a certain proportion of samples in each niche.

To compare the microbiome richness across niches and nitrogen treatment, Alpha diversities were calculated and analysed across niche microbiomes with ASVs (amplicon sequence variants). For all three amplicons, richness of microbiomes were distinct across farm niches (p < 0.001, Kruskal-Wallis rank sum test) (Table 1 and Supplementary Table 2), but were not significantly differentiated between nitrogen treatments (Figure 1A, Supplementary Figures 1A–2A, and Supplementary Table 2). In general, leaf microbiomes were the least diverse where BS microbiomes richness under no-nitrogen were the highest across amplicons, with an average of 788, 357, and 385 observed ASVs from 16S, 18S, and ITS, respectively. In contrast, phyllosphere microbiomes richness were relatively low. For instance WCL microbiome richness under no-nitrogen were low across amplicons, with an average of 36, 20, and 149 observed ASVs from 16S, 18S, and ITS, respectively. Nitrogen treatment did not alter microbiome richness except for ITS of below-ground niche WCR (Supplementary Figure 2A), where microbiome richness was decreased with soil nitrogen treatment (p = 0.045, Kruskal-Wallis rank sum test).

Taxonomic compositions were distinct across niches, but did not change in response to soil nitrogen treatment (Figure 1B, Supplementary Figures 1B–2B, and Supplementary Tables 3A–C). For 16S communities, regardless of nitrogen treatment, the phyllosphere microbiomes in both RGL and WCL were dominated by Proteobacteria (82 and 79% on average, respectively). Similarly in soil and rhizosphere niches (BS, RGR, and WCR), Proteobacteria dominated at 28, 36, and 47% on average, respectively. Animal-associated microbiomes R and F were both dominated by Bacteroidetes and Firmicutes (Supplementary Table 3A). Trends (i.e., changes across niches, with no treatment effects) were similar for 18S and ITS microbiomes (Supplementary Figures 1B–2B and Supplementary Tables 3B,C) with a couple of exceptions. For 18S, Phragmoplastophyta dominated in all phyllosphere and rhizosphere niches, while Cilliophora dominated the 2 animal-associated niches (Supplementary Figure 1 and Supplementary Table 3B). For ITS, Ascomycota dominated across all niches except for R, where Neocallimastigomycota was the most abundant (Supplementary Figure 2 and Supplementary Table 3C).

Bray-Curtis distances were calculated and plotted with NMDS ordinations to examine microbiome beta diversities across niches. Significant differences in community composition were observed across niche microbiomes (16S: ANOSIM: R = 0.758 and ADONIS: R² = 0.505. 18S: ANOSIM: R = 0.246 and ADONIS: R² = 0.248. ITS: ANOSIM: R = 0.647, ADONIS: R² = 0.325, p < 0.01 for all cases), but not in response to nitrogen treatment (Figure 1C, Supplementary Figures 1C–2C, and Table 1) except for WCR (16S community) and R (18S community). For 16S, animal-associated microbiomes were more distinct compared with other niches (Figure 1C). Samples from the rhizosphere were clustered primarily by niches, but partial overlaps were observed between BS, RGR, RGL, WCR, and WCL microbiomes. Microbiome compositions for 18S and ITS were also significantly different across niches, but broad overlaps across niches were found. Nitrogen treatment had minimal impact on microbiome compositions within each individual niche except for 16S microbiomes in WCR (ANOSIM p = 0.004, R = 0.317 and ADONIS p = 0.01, R² = 0.189), and for 18S microbiomes in R (ANOSIM p = 0.027, R = 0.043 and ADONIS p = 0.048, R² = 0.056).

To compliment the microbiome compositional and structural change in response to nitrogen treatment, microbiome structural dynamics were also measured for each amplicon across niches with frequency and occupancy plots (Figure 2 and Supplementary Figures 3, 4). Microbiomes occupancy-frequency distributions across amplicons in both of the phyllosphere niches showed similar skewed patterns. The number of shared ASV declined with increasing number of samples within niche, followed by an increase in species occupying all or most sites, illustrating a core-satellite pattern (Hanski, 1982; Lindh et al., 2017). Compared to the phyllosphere microbiomes, rhizosphere microbiomes showed variations in distribution dynamics across amplicons. Core-satellite patterns were found in ITS and 18S communities, but not in 16S, suggesting variations in community assembly mechanisms and selective pressure across sub-populations and niches.

To identify changes in ASV abundance in response to nitrogen treatment, an Exact test was performed for each niche between two nitrogen treatment pairs: no-nitrogen verses medium-nitrogen, or no-nitrogen vs. high-nitrogen. For each niche, responsive ASVs (Exact test with logFC ≥ 2 or logFC ≤ −2, and BH adjusted p < 0.05) (defined as N-responsive ASVs) were identified. Scattered plots were used to illustrate abundance changes of N-responsive ASV across niches (Figure 3 and Supplementary Figures 5, 6). Volcano plots (Supplementary Figures 7–9) and heatmaps (Supplementary Figures 10–12) were used to provide complementary details on fold changes of N-responsive ASVs across niches and nitrogen treatments.

All N-responsive ASVs were linked with animal-associated niches (Figure 3 and Supplementary Figures 5, 6), but their presences were detected in other niches (Supplementary Figures 10–12). In 16S communities for example, 126 ASV were originally identified as N-responsive in animal-associated niches (Supplementary Table 4A and Figure 3), but none in other niches. Interestingly, over one third (47 out of 124) of N-responsive 16S ASVs were classified under the genus Prevotella. Moreover, all of the N-responsive Prevotella ASVs had positive foldchanges under medium- or high-nitrogen (Supplementary Table 4A). Similarly in 18S and ITS communities, the majority of N-responsive ASVs of (18S: 1692 out of 1694, ITS: 139 out of 140) were only responsive in animal-associated niches (Supplementary Figures 8, 9 and Supplementary Tables 4B,C). Surprisingly, no ASV from BS niche was N-responsive across amplicons. In addition, close to a third of 18S N-responsive ASVs (495 out of 1694) were unclassified, suggesting a lack of reference sequences for N-responsive 18S sequences.

To investigate microbiome network changes in response to nitrogen treatment across niches, network analyses were performed for each niche individually with a pool of unique ASVs. The ASV pool was formed by a conjunction of N-responsive ASVs and the 200 most abundant ASVs (Supplementary Table 5) from each niche. Pairwise Spearman correlations were conducted between selected ASVs (i.e., ASVs found in ASV pool) to calculate their correlations. Only strong correlations (r ≥ 0.8 or r≤ −0.8, p ≤ 0.05, defined as connections) were included in the networks (Supplementary Table 6). Networks in each niche were formed by all connections found with each nitrogen treatment, each network connection was formed by two ASVs (i.e., nodes) and a line (i.e., edge) in between. Overall, microbiome networks were responsive to nitrogen treatment across niches, but responses were inconsistent between niches and across amplicons.

Nitrogen treatment drastically reduced or eliminated network connections in BS across amplicons (Figure 4, Supplementary Figures 13–14, and Supplementary Table 7). On average, BS networks had 4.94, 9.43, and 16.21 edges per node in 16S, 18S, and ITS, respectively, under no-nitrogen. BS network connections were drastically reduced for ITS under medium- and high-nitrogen (2.86 and 3.03 edges per node, respectively), and were eradicated for 16S and 18S under medium- and high-nitrogen.

Nitrogen treatment had knock-on effects on other microbiome networks, but resulted in random changes across niches and amplicons (Figure 4 and Supplementary Figures 13, 14). For 16S, faecal networks under medium-nitrogen (108 nodes and 149 edges) were the most dense compared to networks under no- or high-nitrogen (no-nitrogen: 22 nodes and 20 edges; high-nitrogen: 76 nodes and 59 edges, Figure 4 and Supplementary Table 7). Similar trends were noted in the other 16S animal-associated networks, but not for other amplicons (Supplementary Figures 13, 14 and Supplementary Table 7). Rhizosphere networks responded to nitrogen treatment differently between niches and amplicons. Network connections of WCR in both 16S (18 nodes and 24 edges) and ITS (94 nodes and 700 edges) under high-nitrogen were the most dense compared to that under no- or medium-nitrogen (Figure 4, Supplementary Figure 14, and Supplementary Table 7). In contrast, RGR networks only showed notable changes in 18S, where networks (39 nodes and 55 edges) under medium-nitrogen were the most dense (Supplementary Figure 13). Phyllosphere microbiome networks across amplicons showed no response to nitrogen treatment.

Microbiome diversity and community composition had no response to nitrogen treatment in the soil niche itself or in any other niche within the same farm. This is consistent with previous studies (Eo and Park, 2016; Banerjee et al., 2019; Zhao et al., 2019), but not with many others (Ramirez et al., 2012; Zeng et al., 2016; Wang et al., 2018). Zhou et al. (2016) suggested that long term nitrogen fertilisation significantly decrease fungal diversity and alter fungal compositions in soil, but our results suggested much weaker effects. The disagreement could be caused by the heterogeneity of soil microbiomes between different research sites. Unlike rumen or gut microbiomes, which were relatively more defined, soil microbiomes varied both spatially and temporally (Ettema and Wardle, 2002; Martiny et al., 2011; Pasternak et al., 2013). Besides, certain organisms could be more sensitive to nitrogen treatments than others. Zhao et al. (2019) suggested that the overall fungal community composition and diversity were not affected by nitrogen treatments, but around 10% of taxa were sensitive to nitrogen treatment. Furthermore, only a small number of soil samples were taken in the current study. This potentially reduced the statistical power to detect significant distinctions between microbiomes under different nitrogen treatments.

Plant associated (phyllosphere and rhizosphere) microbial community compositions and diversities were not responsive to nitrogen treatments except 16S microbiomes in WCR. Ledgard et al. (2001) suggested that nitrogen additions (400 kg/ha/yr) reduced the productivity of nitrogen fixation in white clovers. The symbiosis between white clover root and Rhizobium are fundamental for nitrogen fixation (Kneip et al., 2007), and nitrogen fertilisation can make this capability redundant, thus leading to compositional changes in root communities, such as a lower proportion of nitrogen fixers.

Our results showed that soil nitrogen treatment did not affect microbiome richness and compositions in animal-associated niches. It was unlikely that animal-associated microbiomes were experiencing direct changes of nitrogen availability. In contrast, the soil nitrogen treatment would most likely influence ruminants by changing the plant materials grazed on the pasture (Henderson et al., 2015). Applications of nitrogen fertiliser was proven to favour pasture growth while the production of clover was unfavourable (Moir et al., 2003), leading to a higher ryegrass: white clover ratio. This vegetation change would then lead to changes in nutrient intake to a higher fibre and less protein diet for ruminants (Rattray and Joyce, 1974). Bowen et al. (2018) and Smith et al. (2020) showed that changes in sward types (i.e., ryegrass only vs. combined ryegrass and white clover) altered the rumen microbiome composition in dairy cows. In contrast, O’Callaghan et al. (2018) suggested that changes in feedings between ryegrass only and mixture of ryegrass and white clover had no impact on rumen microbiome composition, but altered the rumen metabolome. Unlike the previous microbiome studies on sward type alterations (i.e., with or without white clover) (Bowen et al., 2018; Smith et al., 2020), we did not expect to see a huge change in ryegrass and white clover ratio between control and fertilised ground [i.e., about 5–10% reduction in white clover on fertilised ground (Moir et al., 2003)]. This small change in vegetation might not lead to a detectable change in animal microbiomes.

Many N-responsive ASVs were present across niches (Supplementary Figures 10–12), but most of N-responsive ASVs only had strong correlations in animal-associated niches across amplicons, and few N-responsive ASVs in phyllosphere or rhizosphere niches (Supplementary Figures 5, 6). Zancarini et al. (2012) suggested that soil nitrogen treatments only affected rhizosphere microbiomes with certain host plant genotypes (Zancarini et al., 2012). Therefore, our results could mean that white clover and rye grass rhizosphere microbiomes were insensitive to soil nitrogen treatments. Beattie and Lindow (1999) stated that phyllosphere microorganisms were generally oligotrophs which can tolerate low-nutrient conditions or formed symbiosis with host plant to obtain more nutrients. When obtaining nutrient from soil, physical movement of nutrients, most likely from host roots to leaf surface, was inevitable for phyllosphere microbiomes. Hence, phyllosphere microbiomes might not experience much knock-on effects caused by soil nitrogen treatments due to physical barriers.

Among 16S N-responsive ASVs, over one-third (47 out of 126) were classified as Prevotella spp., all of which were enriched under medium- or high-nitrogen compared to no-nitrogen, suggesting that Prevotella spp. were potentially benefited from soil nitrogen treatment, or vegetation change (i.e., 5 – 10% reduction of white clover). Compared to white clover, ryegrass contained less protein and more soluble and structural carbohydrates (Rattray and Joyce, 1974). Prevotella was known to be an abundant and essential group of catabolic generalists in the rumen (Hungate, 2013; Emerson and Weimer, 2017), due to the ability of degrading a variety of carbohydrate substrates, including hemicellulose, xylan, and pectin (Dodd, 2010). Therefore, the enrichment of N-responsive Prevotella spp. was most likely due to increases in carbohydrates, or in other words higher ryegrass intake by cows. Unlike Prevotella spp., patterns were random for other N-responsive taxa. This suggested that the organisms under the same taxa groups could respond to the same disturbance differently. For example, Trichostomatia is a common group of rumen protozoa fermenters involved in degrading non-structural polysaccharides and soluble sugar (Wright, 2015). This aligned with our results as 696 out of 1199 N-responsive ASVs were Trichostomatia ASVs. Furthermore, because the range of carbohydrates metabolised was genus-dependent (Wright, 2015), it was unsurprising to see inconsistent change in N-responsive Trichostomatia ASVs.

Our results showed that for 16S and 18S communities, network connections in soil were completely eradicated under medium- and high-nitrogen, whereas ITS networks complexity and connectivity were drastically reduced but still detectable. Our findings aligned with a previous study showing agricultural fertilisation associated with decreases in microbial network complexity (Banerjee et al., 2019). The dramatic reduction in network linkages in response to nitrogen could be due to the redundancy of nitrogen-fixation functions. When nitrogen was limited in soil, microorganisms formed networks to transform inaccessible nitrogen into an accessible form for non-nitrogen fixers, supporting community growth (Kuypers et al., 2018). Under medium- or high-nitrogen, potentially non-nitrogen fixers would no longer require nitrogen fixation from nearby nitrogen fixers, leading to an alteration of microbiome networks. Instead, alternative nitrification pathways, such as aerobic ammonia oxidation, could potentially be induced and consequently led to population and activity increases of ammonia-oxidizing microbes (Di et al., 2009; Dai et al., 2013). However, further investigations on the transcriptome are required to corroborate the conjecture.

Other than reduction of soil network connectivity, ITS community networks were unresponsive to nitrogen treatment across niches, suggesting that fungal communities were less sensitive to nitrogen treatments compared to prokaryotes and other eukaryotes in general. However, since the current study was only based on amplicon sequencing, the presence of ‘relic fungal DNA’, such as dormant fungal spores, was potentially interfere with detection of relevant fungi and confound the results.