Soil collection and glasshouse experiment
Soil samples used in this study were collected in March 2022 from a wheat experimental farm on the Dookie Campus at The University of Melbourne, Australia (36.382° S, 145.711° E). The experimental field has a history of wheat-pasture rotation. We collected approximately 400 kg of soil (0–10 cm) using a shovel after removing the root debris for the pot experiment. The soil is classified as a Dermosol (Isbell 2024) with the following physicochemical properties (measured by the Nutrient Advantage lab, Victoria, Australia): total C, 0.021 g kg− 1 soil, total N, 0.002 g kg− 1 soil, organic C, 0.018 g kg− 1 soil, nitrate-N, 56 mg kg− 1 soil, ammonium-N, 2.9 mg kg− 1 soil, and soil pH of 6.3. The collected soil samples were sieved to 1 cm, air-dried, and homogenised before the glasshouse experiment. Pots with an inner diameter of 16.5 cm were filled with 4.6 kg of soil. Two wheat cultivars, high-NUE Mace and low-NUE Gladius, provided by the Australian Grain Technology (AGT) company, were selected for the pot experiment. Five uniform-sized seeds were sown at a depth of 3 cm and seedlings were thinned to three per pot at the two-leaf stage. To establish varying N levels, ammonium sulphate was applied at three rates: (1) N0, without N addition, (2) N1, 45 kg N ha− 1, and (3) N2, 90 kg N ha− 1, with eight replicates for each N treatment. The soil moisture content was maintained at around 80% field capacity throughout the experiment. The plant growth conditions were maintained at 24°C /10°C for 12 h/12 h day/night, 65% humidity and ~ 225 µmol m− 2 s− 1 during the daytime throughout the study.
The wheat rhizosphere and bulk soils, as well as plant roots, were destructively sampled at the stem-elongation and the mid-anthesis stages. The wheat root system was placed in a sampling bag and subjected to gentle shaking to obtain rhizosphere soil samples. Soils remained in the pots were randomly collected as the bulk soil. Wheat roots (5–20 cm) were separated from the plant and were gently rinsed with reverse osmosis (RO) water to remove large soil particles. Root samples were subsequently shaken with 30 ml of phosphate-buffer solution (PBS) in a falcon tube at 250 rpm to remove any adhered soil and briefly dried with filter paper. The soil and root samples were then stored at -20°C.
Plant and soil physicochemical characterisation
The fresh above ground plant materials were oven-dried at 60°C for 7 days for biomass measurement. The dried plant material was ground to a 0.2-mm fine powder using the ultra-centrifugal mill (ZM200 Ultra Centrifugal Mill, RETSCH). A set of 48 bulk soil samples (2 growth stages, 2 cultivars, 3 N levels and 4 replicates) were oven-dried at 40°C for 7 days, and then ball milled using a tissue grinder (TissueLyser II, QIAGEN). The resulting, finely ground plant and soil samples were subjected to the CN analyser (LECO TruMac Series) for total C and N measurements (Figure S1). Another set of 48 bulk soil samples were oven-dried at 105°C for 2 days for the estimation of soil moisture content. For soil pH and electrical conductivity (EC) measurements, an alternative set of fresh soil was suspended in Milli-Q water at a 1:5 ratio and shaken at 250 rpm for 60 min. Soil pH was measured with a professional benchtop pH meter (HANNA) and EC was measured with the conductivity meter (Multiparameter Laboratory Benchtop Conductivity meter, SmartCHEM-LAB).
DNA extraction and amplicon sequencing
Frozen root samples were ground into a fine powder with liquid N. The total DNA was extracted from 250 mg of the frozen rhizosphere soil and root samples using the DNeasy PowerSoil Pro Kits (QIAGEN), and quality was assessed using Nanodrop (NanoDrop One, Thermo Scientific). A total of 48 rhizosphere and 48 root endosphere DNA samples were subjected to polymerase chain reaction (PCR) amplification. The V3–V4 hypervariable region of the bacterial 16S rRNA gene was amplified with the 341F/806R primer pairs (Klindworth et al. 2013) and sequenced using the Illumina MiSeq System platform at the Australian Genome Research Facility, Melbourne, Australia.
Amplicon sequencing was processed using QIIME 2 (v.2022.11) (Bolyen et al. 2019). PCR primer sequences were removed from the raw sequences using cutadapt (v.4.2) (Martin 2011) and were filtered with an average quality score cutoff > 20 using DADA2 (v1.26.0) (Callahan et al. 2016). Filtered reads were further subject to DADA2 for denoising (Callahan et al. 2016). Non-chimeric amplicon sequence variants (ASVs) were obtained by merging read pairs that overlapped at least 12 bp and removing PCR chimeras. Taxonomic assignment of ASVs was performed with QIIME 2 (v.2022.11) (Bolyen et al. 2019) using the SILVA database (v.138.1) (Quast et al. 2013) for the 16S rRNA gene of the rhizosphere (Supplementary data 1) and root endosphere samples (Supplementary data 2). A minor proportion of archaeal ASVs were detected and their alpha and beta diversity had no significant differences across samples (Figure S2 & S3). Hence, ASVs belonging to mitochondria, chloroplast, archaea and Triticum aestivum were removed before the downstream analysis. The count matrix was rarefied using microeco (v.1.1.0) (Liu et al. 2021) to ensure even read depth for sample comparisons.
Metagenomic analysis
To further investigate the functional changes of rhizosphere bacterial communities in response to N deficiency, shotgun metagenomic sequencing was performed on the extracted DNA of the rhizosphere samples of Mace and Gladius under the two contrasting treatments, N0 and N2. A total of 16 metagenomic libraries of the rhizosphere soil DNA were generated with the NovaSeq 6000 sequencing platform using the Nextera Flex primer (Bruinsma et al. 2018), yielding 2 × 150 bp reads at the Australian Genome Research Facility, Australia. The resultant metagenome sequences (24 Gb per sample) were subjected to quality checks using fastqc (v.0.11.9) (Andrews 2010) and adapter trimming using trimmomatic (v.0.39) (Bolger et al. 2014). Samples were decontaminated using the Triticum aestivum genome (IWGSC) retrieved from the EnsemblPlants (Yates et al. 2022) using bowtie2 (v.2.4.2) (Langmead and Salzberg 2012). Clean reads were processed using the Metaphor (v.1.7.9) pipeline (Salazar et al. 2022). Two co-assemblies were generated by pooling the reads from samples of the same cultivar and assembled with MEGAHIT (v.1.2.9) (Li et al. 2015). Genes were predicted over the assembled contigs using Prodigal (v.2.6.3) (Hyatt et al. 2010) and functionally annotated with DIAMOND (v.2.1.0) (Buchfink et al. 2015) against the Clusters of Orthologous Genes (COGs) (v.1.0) database (Tatusov et al. 2000).
Metaphor outputs the relative abundance of contigs falling in different COGs functions (Supplementary data 3) and categories (Supplementary data 4). These files were used to compare the functional similarity between treatment samples. We integrated the contig count and COGs annotation output from Metaphor to explore the N metabolism functions in bacterial genera. Contig depth was estimated by the reads per kilobases (RPK) method (Supplementary data 5) and the contig annotation file contained COGs taxonomy and function information generated from DIAMOND (Supplementary data 6). We first filtered out contigs related to N-metabolism and then searched for the presence of the bacterial genera in Mace or Gladius under N deficiency, in their corresponding samples.
Hydroponic experiment
The seeds of two wheat cultivars were sterilised with 80% ethanol for 1 min, followed by 2.5% sodium hypochlorite (NaOCl) for 15 mins and rinsed with Milli-Q water for 5 times. Sterilised seeds were germinated on petri dishes with Milli-Q water inside a fume cupboard, covered in dark, at room temperature, until the two-leaf stage was reached (Lu et al. 2021). Two uniform seedlings of each cultivar were transferred into a 4.3-L hydroponic system. The hydroponic pots were filled with a modified Hoagland solution with the following composition: 1 mM KH2PO4; 1 mM MgSO4; 0.05 mM H3BO3; 0.01 mM Fe-EDTA; 0.009 mM MnSO4; 0.0007 mM ZnSO4; 0.0003 mM CuSO4; 0.0001 mM NaCl and 0.0001 mM H2MoO4. Three N levels, (1) N0, 2 mM, (2) N1, 4 mM and (3) N2, 8 mM were achieved by adding various amounts of Ca(NO3)2·4H2O, KNO3 and (NH4)2SO4, with a NH4+ to NO3− molar ratio of 1:4 across all N treatments (Thomas and Paparozzi 2013). Additional K2SO4 and CaCl2 were introduced to balance the levels of K+ and Ca2+ in the solution. The hydroponic solution was renewed with a half-strength Hoagland solution once every two weeks before the stem elongation stage and with a full-strength solution weekly thereafter. The solution pH was maintained within the range of 5.5–6.5 by adding 0.1 M KOH, while the water level was maintained using RO water. An aeration system was implemented in the hydroponic pots at the onset of tillering, and the arrangement of hydroponic pots was randomised on a weekly basis.
During the mid-anthesis stage, wheat root samples were collected and rinsed with 0.2 mM CaCl2, followed by immersion in 800 mL of 0.2 mM CaCl2 for 3 h with bubbling to remove adhered salt and microbes from the root surface (Liu et al. 2019). Root was separated from the plant and immediately stored in a freeze room at 20°C. Frozen root samples were freeze-dried for 2 weeks, cut into ~ 3 cm pieces, and homogenised (Saiman et al. 2012). Freeze-dried root samples were ground to fine powder with liquid N and stored at -20°C for metabolite extraction.
GC-MS analysis for polar metabolites
Organic compounds in root were extracted from 60 ± 2 mg of the stored lyophilised and ground materials using 500 µl of 100% MeOH containing 4% (v/v) of [13C6] sorbitol/[13C515N] valine (0.5 mg ml− 1). The supernatant was stored in a 2-ml reaction tube as the first extraction product. A second extraction was performed on the remaining pellet with 500 µl of Milli-Q water, vortexed and centrifuged. The two extracted supernatants were combined and 100 µl of the combined supernatant was dried using a speed vacuum dryer (John Morris Scientific Pty Ltd). Derivatization of the dried extracts was performed with 20 µl of methoxyamine hydrochloride in pyridine and bis-(trimethylsilyl)-trifluoroacetamide (BSTFA), as described by (Dias et al. 2015). Specifically, a solution of 30 mg ml− 1 methoxyamine hydrochloride in pyridine was added to all samples, followed by a derivatization process at 37°C for 120 minutes with agitation at 500 rpm (61.0 g). Subsequently, 20 µl of BSTFA was added to the samples and were further incubated on a thermomixer (500 rpm/61.0 g) at 37°C for an additional 30 minutes. Samples were settled for 60 minutes before the injection.
The samples were analysed with a GC-MS system comprised of an autosampler (Gerstel 2.5.2), a gas chromatograph (7890A Agilent), and a quadrupole mass spectrometer (5975C Agilent) (Agilent, Santa Clara, United States). An injection volume of 1 µl was used for each derivatized sample. The mass spectrometer was tuned using tris-(perfluorobutyl)-amine (CF43) following the manufacturer’s recommendations. Chromatography was carried out on a 30 m column (J&W VF-5MS Agilent) with a film thickness of 0.25 µm and an internal diameter of 0.25 mm. A 10 m Integra guard column was additionally equipped. The inlet temperature was set at 250°C; the transfer line of MS was at 280°C; the ion source was adjusted to 230°C; and the quadrupole was maintained at 150°C. Helium was served as the carrier gas at a flow rate of 1 ml min− 1.
Metabolite data pre-processing
Data were pre-processed in silico for mass spectral deconvolution, with peak picking and identification using the Automated Mass Spectral Deconvolution and Identification System (AMDIS) software (National Institute of Standards and Technology, Gaithersburg, MD, USA). The target component library was created with the PBQC file and peak identification was confirmed with the Agilent MassHunter Qualitative Analysis B.05.00 (Agilent Technologies, Inc., 2011). The peak area data was exported with the Agilent MassHunter Quantitative Analysis software version B.08.00/Build 7.0.457.0 (Agilent Technologies, Inc., 2008). The peak area data was normalised against the internal standard [13C6] sorbitol and the weight of root samples to generate response values, which were log-transformed for the downstream analysis (Gupta et al. 2022) (Supplementary data 7).
Statistical analyses
All statistical analyses were performed using R (v.4.2.1) and all figures were generated with ggplot2 (v.3.4.4). Wheat growth parameters were statistically compared using the analysis of variance (ANOVA). The number of samples used for statistical analysis is provided in Table S1. Shannon index was used to estimate the bacterial diversity and the differences between treatments were compared using the nonparametric Kruskal-Wallis test. Principal coordinates analysis (PCoA) based on the Bray-Curtis distance matrix was performed to assess the dissimilarity of the bacterial community compositions between treatments. Analysis of similarity (ANOSIM) with the ‘vegan’ package (Oksanen et al. 2022) was performed to statistically compare the similarity of bacterial community between treatment samples. Bacterial biomarkers were determined between N0 and N2 samples at the mid-anthesis stage. The linear discriminant analysis effect size (LEfSe) was quantified at the ASV level using the count per million (CPM) transformation with microeco (v.1.1.0). ASVs were considered significantly enriched with the criteria of p < 0.05 and a linear discriminant analysis (LDA) effect score > 2. For the metagenomic data, the relative abundance of COG functional code from Metaphor was subjected to the PCoA and ANOSIM analyses (Supplementary data 3). The relative abundance of COG categories was further compared between the N0 and N2 treatments using the Kruskal-Wallis test (Supplementary data 4). For the GC-MS metabolomics analysis, unsupervised learning principal component analysis (PCA) was performed to assess the similarity of metabolite profiles between treatments. The metabolite composition between treatment samples was statistically compared with the Permutational Multivariate Analysis of Variance (PERMANOVA) using the ‘vegan’ package (Oksanen et al. 2022). The fold-change of non-redundant known metabolites was analysed on the response values through multiple Student’s t-tests with false discovery rate (FDR) adjustment (adj. p < 0.05). A significance threshold of p < 0.05 was applied to all statistical results.