Soil collection and glasshouse experiment
Soil samples (0–20 cm) 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 soil is classified as a Dermosol (Isbell 2021) with the following physicochemical properties: 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 were selected for the pot experiment. 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. Wheat roots were gently rinsed with RO water to remove large soil particles, and subsequently shaken with 30 ml of phosphate-buffer solution (PBS) in a falcon tube at 250 rpm to remove any adhered soil. The soil and root samples were then stored at -20°C.
Plant and Soil Physicochemical Characterisation
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). Soil samples 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. 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, the 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
Root samples were ground into a fine powder with liquid nitrogen. The total DNA was extracted from soil and root samples using the DNeasy PowerSoil Pro Kits (QIAGEN), and quality was assessed using Nanodrop (NanoDrop One, Thermo Scientific). 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 (AGRF), Melbourne, Australia.
PCR primer sequences were removed from the raw sequences using cutadapt (v.4.2) (Martin 2011) and were quality-filtered to remove reads with an average quality score < 20 using DADA2 (v1.26.0) (Callahan et al. 2016). Filtered reads were further subject to DADA2 for demultiplexing and 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 (Bolyen et al. 2019) using the SILVA (v.138.1) (Quast et al. 2013) for the 16S rRNA gene. 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
Shotgun metagenomic sequencing was performed to investigate the functional differences in the rhizosphere microbiome under the N0 and N2 treatments. 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. 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). The Metaphor outputs were used for statistical analysis.
Quantitative PCR (qPCR) analysis
The abundance of root-associated bacteria was characterised with the 506F/816R primers (Caporaso et al. 2011). The 16S rRNA gene abundance was estimated with a standard curve generated by a 10-fold serial dilution of a plasmid containing a full-length copy of the Escherichia coli 16S rRNA gene. The 20 µl qPCR reactions reagent contained 10 µl SYBR MasterMix (Thermo Fisher), 0.8 µl of each 10 µM forward and reverse primers, 6.4 µl sterile DNA-free water, and 2 µl DNA template. The reaction was carried out using the real-time PCR system (CFX96 Touch Real-Time PCR Detection System, BIO-RAD).
Hydroponic experiment
The seeds of two wheat cultivars were surface sterilised with 80% ethanol and germinated in the 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 (2, 4 and 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). Roots were cut into ~ 3 cm pieces, homogenised and freeze-dried for 2 weeks. Freeze-dried root samples were ground to fine powder with liquid N and stored at -20°C for metabolite extraction.
Untargeted 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. Before the injection, the samples were rested for 60 minutes.
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).
Statistical analyses
All statistical analyses were performed using R (v.4.2.1) and all figures were generated with ggplot2 (v.3.4.4). 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. Adonis analyses with the Permutational Multivariate Analysis of Variance (PERMANOVA) using the ‘vegan’ package (Oksanen et al. 2022) were further performed to statistically compare the bacterial community compositions between treatments. 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 contig count matrix by COG functional code from Metaphor was subjected to the PCoA and PERMANOVA analyses. The contig count matrix by COG categories was further compared between the N0 and N2 treatments using the Kruskal-Wallis test. For the GC-MS untargeted 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 further compared with the PERMANOVA analysis. The fold-change of non-redundant known metabolites was analysed on the response values through multiple Student’s t-tests with False Discovery rate adjustment (adj. P < 0.05). A significance threshold of P < 0.05 was applied to all statistical results.