Experimental design

We performed a greenhouse experiment with common bean (Phaseolus vulgaris var. Red Hawk) and switchgrass (Panicum virgatum var. Cave-in-Rock). The experimental design consisted of two crop levels (bean, switchgrass), two planting levels (planted, unplanted), two moisture levels (drought, watered), and five destructive plant and rhizosphere sampling points correlating with increased drought severity. The drought gradient included a pre-drought baseline that was watered, and then destructive sampling over the next six days, with the day 2 samples experiencing the lowest drought severity and day 6 samples experiencing the highest. Five replicates were collected at each time point for each combination of the factors tested (Fig. 1).

Soil sample collection

Bean soil was collected from Montcalm Research Center on September 6, 2018, at Stanton, MI (43.350885, −85.177044), and was most recently planted with common bean (Phaseolus vulgaris var. Red Hawk) that year. Switchgrass samples were collected from the Great Lakes Bioenergy Research Center switchgrass plots at Lux Arbor Reserve (42.475224, −85.444979) in Denton, MI, on August 27, 2018. This site has been under continuous switchgrass (Panicum virgatum var. Cave-in-rock) cultivation since 2011. Switchgrass soil was collected and homogenized from six sites randomly situated throughout the plots. All sampling materials were sterilized with ethanol before use and between samples. A shovel was used to collect bean and switchgrass soils to a depth of 10 cm. Large rocks and root fragments were removed from the soil. Within three hours (for bean soils) and five hours (for switchgrass soils), soils were transported to the lab and stored at 4 °C until sieving. Soil nutrient analysis was completed for the field soils at the Michigan State University Soil and Plant Nutrient Laboratory (Supplementary Table 4).

Soil samples were sieved (4 mm mesh) over three to five days immediately after the sample collection day and stored at 4 °C until needed for the experiment. The drought experiment commenced after eight months and six months of soil storage at 4 °C, for bean and switchgrass, respectively.

Seed germination

Approximately 250 switchgrass seeds (Sharp Bros. Seed of Mo., Inc) were placed in 50 ml conical tubes. 10 ml of 5% bleach was added to the tubes. The tubes were placed on a shaker at 200 rpm for 20 min. The bleach/water mix was decanted, and 10 ml of deionized water was added. The tubes were agitated again at 200 rpm for 5 min, then decanted. The deionized water rinse was repeated two more times. Approximately 60 switchgrass seeds were placed onto a wet filter paper (90 mm diameter) per sterile Petri plate. Seeds were moistened with deionized water, plates were wrapped in parafilm, covered with foil to block light, and incubated at 32 °C incubator for three days. After three days, 1 ml of deionized water was supplemented to each plate, which was re-wrapped and returned to 32 °C for an additional day.

Bean seeds were obtained from the Dry Bean Breeding and Genetics program at Michigan State University. Approximately 590 seeds were placed in a 1 L Erlenmeyer flask, and 1 L of sterilization solution was added (0.1% Tween-20, 10% bleach, and 90% DI water) to the seeds. The seeds were soaked for 15 min at room temperature and inverted every few minutes. The seeds were washed five times with deionized water, and 12 seeds were placed onto Petri dishes lined with moistened filter paper. Plates were wrapped in parafilm and incubated at room temperature in the dark. Over the next three days, 2 ml of Milli-q water was added to each plate. After four days, approximately 60% of the switchgrass and bean seeds germinated and had an emerged radicle, and these were advanced for use in the experiment.

Growth tube preparation and planting

Growth tubes were made by modifying 50 ml conical tubes to better control soil water content during the experiment. An ethanol-cleaned drill bit was used to drill a hole in the bottom of 50-ml conical tubes. A hole was also drilled in every cell of the tube racks used to hold the 50 ml tubes for switchgrass and bean plants. Small squares of Kimwipe® tissue papers (5 cm * 5 cm) were cut, each rolled, and one rolled sheet was used to plug the drain hole in each 50 ml conical tube (Supplementary Fig. 1). Autoclaved perlite was combined with sieved field soil (50:50 v/v) to fill the prepared tubes to the top. Soils were watered with 10 ml of deionized water, which reduced the total soil volume to approximately 35 ml. Tubes for the planted treatment received one germinated seed. Unplanted tubes were prepared the same but without a seed. Seeds were planted so that ~1 cm of the seedling was exposed above the soil and at a depth to ensure that the emerging root would not push the seedling out of the soil. After planting, additional soil: perlite mix was added to the tube to replace the volume lost after settling and supplemented with 3 ml of deionized water. Each tube was wrapped in foil to exclude light in the root zone. The greenhouse’s daily temperature range was 21-32 °C, with 14 h of supplemental daytime lighting (400-watt high-pressure sodium lamp).

Seedlings were provided sufficient water to achieve healthy growth before initiating the drought. Switchgrass tubes were watered with 10 ml of deionized water every other day or every two days for two months until the plants were large enough for the experiment, which was designated as when the roots had filled the tube. Though this was close to the flowering stage, no flowers were observed. Bean tubes were watered daily with 10 ml of deionized water for the first five days and then twice daily with 10 ml for four additional days. After ~9 days of watered conditions, approximately 55 bean plants were sorted by height as a proxy for total biomass, and similarly sized plants were selected to include in the experiment. At this vegetative stage, the first trifoliate leaf emerged and started to unfold in all bean plants, and roots were completely bound within the tube. Before the drought, the planted tubes of both switchgrass and bean had tube-bound root systems such that all their soil could be directly influenced by the plant and classified as the rhizosphere. Bean roots were tube-bound within a week, but switchgrass took about two months. The timing for the start of drought treatment was chosen based on the time needed for the roots to fill the tube so that the entire soil was plant-influenced rhizosphere. Thus, the growth stage at the time of drought commencement was different between the two crops.

Drought treatment and destructive sampling

The same experimental design and sampling strategy were applied to bean and switchgrass (Fig. 1). Twenty-five seedlings and 25 unplanted tubes were assigned to either the sufficiently watered (“watered”) or to the reduced water treatment (“drought”), totaling 100 samples. In addition, five planted and five unplanted tubes were collected before the drought (“pre-drought,” day 0). These samples were used as a baseline control to compare to the subsequent drought samples. The drought was initiated on the same day the pre-drought samples were collected (see Supplementary Tables 5 and 6 for the watering regimes for switchgrass and bean, respectively). The drought severity was increased over six days by adding progressively less water. Five replicate soil samples (individual tubes used for destructive sampling) were collected for each experimental condition from day 2 to day 6 for five post-drought sampling points.

To ensure that the drought bean plants remained viable (and that our assessment was of the rhizobiome response to live plants rather than to dead), an additional 36 plants were in parallel subjected to drought to determine their viability on the final days of the experiment. On day 5, when bean plants were browning and lacked turgor (day 5), half (18) of the extra plants were provided 15 ml of deionized water in the morning and the afternoon. All 18 plants were resuscitated, and we observed their leaves returned to full turgor by the next morning. Thus, the drought plants from day 5 were likely alive when sampled despite appearing well-desiccated. The same was repeated on day 6 with the remaining drought plants. Thirteen of the 18 plants resuscitated (leaves back to full turgor) within two days, affirming that most drought plants on day 6 were viable. As a relatively more drought-tolerant plant, switchgrass did not appear severely desiccated on days 5 and 6, and the resuscitation check was unnecessary to confirm its viability.

Soil, root, and shoot sample collection and gravimetric soil moisture

Materials and surfaces were sterilized and/or ethanol-cleaned before use and between samples. All tubes were destructively sampled. During each sampling, each tube’s contents were dumped onto ethanol-cleaned aluminum foil. The soil was brushed from the root system using ethanol-cleaned gloves and homogenized with an ethanol-cleaned metal spatula. Gloves were ethanol-cleaned again when handling new samples and new ethanol-cleaned gloves were used as needed. Roots were manually removed from the soil, and 0.5 g of soil was collected into a 1.5 ml microcentrifuge tube for DNA/RNA co-extractions. An additional ~6 g soil was retained as a backup. Perlite ( > 2 mm) was avoided. Soils were flash-frozen in liquid nitrogen and stored at −80 °C until nucleic acid extractions. An additional 4 g of soil was collected to determine gravimetric soil moisture content. After recording the mass of fresh soil collected, soils were dried in an oven at 50 °C for at least three days and then re-weighed. Percent gravimetric soil moisture was calculated as the percent water mass lost during drying using Eq. 1⁶².

After soil collection, the shoot system was cut at the base of the stem and put in pre-labeled envelopes. The shoots were dried at 50 °C for at least five days. Masses taken after drying were used to determine the biomass of the shoot. Treatment differences in shoot biomass and gravimetric soil moisture were analyzed with a three-way ANOVA type III test upon satisfaction of normality.

An additional 30 switchgrass plants (15 watered and 15 drought) were grown for metabolite analysis on day 6. Soils were homogenized, shoots were sampled as previously described, and root systems were cleaned in deionized water. Then, samples were flash-frozen in liquid nitrogen in 15 ml tubes. Soils of the unplanted treatments (watered and drought) from day 6, stored at −80 °C, were used to compare with the planted treatments. We also analyzed metabolites from pre-drought soils (four replicates of planted and five unplanted treatments) stored at −80 °C. In addition, ten unplanted soils from day 6 (five drought and five watered) and nine pre-drought samples (four planted, five unplanted) were analyzed for metabolites. In summary, 15 switchgrass plants from the day 6 watered and drought treatments each yielded 30 samples for soil, shoot, and root metabolites.

Soil DNA/RNA co-extractions

DNA/RNA co-extractions were performed using the protocol specified in ref. ⁶³ with minor modifications. Briefly, 0.5 g of flash-frozen soil was added to 0.7 mm PowerBead® garnet bead tubes (Qiagen, Germantown, MD, purchased early 2021 before being discontinued). 0.5 ml of CTAB-phosphate buffer (120 mM, pH 8) and 0.5 ml of phenol:chloroform: isoamyl alcohol solution (25:24:1, nuclease-free, Invitrogen®) were added to the tube and placed on a bead beater to beat for 30 s. Tubes were then centrifuged at 12000 g for 10 min at 4 °C. The top aqueous layer was extracted and placed into a new 1.5 ml microcentrifuge tube. Next, 0.5 ml of chloroform: isoamyl alcohol (24:1, nuclease-free, Sigma®) was added to the tube and inverted several times to form an emulsion to remove residual phenol. Tubes were centrifuged at 16000 g for 5 min at 4 °C, the top aqueous layer extracted and placed into a new 1.5 ml microcentrifuge tube. Nucleic acids were precipitated by adding two volumes of 30% polyethylene glycol solution (PEG6000, 1.6 M NaCl) and mixing a few times. Tubes were incubated on ice for two hours. Tubes were then centrifuged at 16000 g for 20 min at 4 °C. The supernatant was removed, and 1 ml of 70% ice-cold ethanol was added. Tubes were centrifuged at 16000 g for 15 min at 4 °C. Ethanol wash was pipetted out, being careful not to remove the pellet, and placed back in the centrifuge for a final spin for 10 s to collect residual ethanol. The remaining ethanol was pipetted out of tubes, and the pellet was allowed to air dry. The nucleic acid pellet was resuspended in 30 μl of nuclease-free water. Negative controls for extractions included tubes to which no soil sample was added, and only included reagents and the garnet beads to check contamination in the extraction reagents. Two negative controls were processed alongside the experimental samples for each extraction day. All DNA and RNA samples were quantified using Qubit® dsDNA BR assay kit and RNA HS assay kit on a qubit 2.0 fluorometer (Invitrogen, Carlsbad, CA, USA). All nucleic acid raw coextracts (containing DNA and RNA) were visualized using agarose gel electrophoresis and validated with a band for DNA and RNA⁶³.

DNase treatment and cDNA synthesis

The DNA/RNA coextract was used to prepare purified RNA using the Invitrogen TURBO® DNA-free Kit with minor modifications. 1 μl of 10X TURBO® DNase Buffer and 3 μl of TURBO® DNase enzyme were added to a 6 μl aliquot of the DNA/RNA coextract. The mixture was then incubated for 30 mins at 37 °C, after which 2 μl of a DNAse inactivation reagent was added to each tube. The resulting solution was mixed by flicking the tubes by hand and then incubated at room temperature for 5 mins. The samples were centrifuged at 2000g for 5 min at room temperature. The purified, DNAse-treated RNA samples were transferred to a clean, sterile tube and immediately processed to make complementary DNA (cDNA). cDNA was prepared using the Invitrogen® SuperScript III First Strand Synthesis System using random hexamers per the manufacturer’s instructions. cDNA samples were stored at −20 °C until further processing. Negative controls were prepared for the cDNA synthesis step to check for reagent contamination.

PCR and RT-PCR

DNA (PCR) and cDNA samples (RT-PCR) were amplified using a standard protocol specified by the Research Technology Support Facility at Michigan State University. A 15 μl reaction volume was prepared for each sample to amplify the V4 hypervariable region of the 16S rRNA gene with 2X GoTaq Green Mastermix, primers 515 F (5′- GTGCCAGCMGCCGCGGTAA- 3′) and 806 R (5′-GGACTACHVGGGTWTCTAAT-3′)⁶⁴ and 1 μl template. The final concentration of the GoTaq Green Mastermix was 1X, and the final concentration of each primer was 0.1 μM. The PCR cycle was run using the following cycling parameters: initial denaturation at 94 °C for 3 min; 30 cycles of denaturation at 94 °C for 45 s, annealing at 50 °C for 60 s, and extension at 72 °C for 90 s, followed by a final extension at 72 °C for 10 min. Samples were kept at 4 °C and immediately visualized on a 1% agarose gel using 100-bp ladder. All RNA samples purified from DNA/RNA coextracts and used to prepare cDNA were also used to run PCR and check on 1% agarose gel to ensure there were no contaminating DNA bands in RNA samples. All PCR reactions included a no-template negative control and an E. coli DNA template as a positive control. All negative controls used for extractions were also included as samples in the PCRs. Results indicated no contamination from extraction or PCR reagents and consistent amplification performance with the cycling parameters.

Illumina sequencing of the 16S rRNA gene and 16S rRNA

Amplicon sequencing was performed at the Genomics Research Technology Support Facility at Michigan State University. The V4 hypervariable region of the 16S rRNA gene was amplified using dual-indexed Illumina compatible primers 515F and 806R as described in ref. ⁶⁵. PCR products were batch normalized using an Invitrogen SequalPrep DNA Normalization Plate, and the normalized products recovered from each of the six plates submitted were pooled. The pools were cleaned up and concentrated using a QIAquick PCR Purification column followed by AMPureXP magnetic beads; it was quality controlled and quantified using a combination of Qubit dsDNA HS, Agilent 4200 TapeStation HS DNA1000, and Invitrogen Collibri Library Quantification qPCR assays.

Each plate submitted for sequencing had approximately 85 samples (6 plates); each plate was used for one MiSeq run. The pools from each plate were each loaded onto an Illumina MiSeq v2 standard flow cell, and sequencing was performed in a 2 x 250bp paired-end format using a MiSeq v2 500 cycle reagent cartridge. Custom Sequencing and index primers were added to appropriate wells of the reagent cartridge. Base calling was done by Illumina Real Time Analysis (RTA) v1.18.54, and the output of RTA was demultiplexed and converted to FastQ format with Illumina Bcl2fastq v2.20.0. In addition to experimental samples, all negative controls from each extraction day were sequenced, and these samples were evenly distributed throughout all the MiSeq runs. Two positive controls (mock communities) were sequenced with each MiSeq run. One positive control was an in-house Mock community prepared in the Shade Lab⁶⁶, and the other positive control was provided by the RTSF at MSU for library preparation.

Metabolite extraction from plant tissue and soil

The switchgrass shoot, root, and rhizosphere soil samples collected from 15 individual plants in the drought experiment were pooled into seven replicates (Supplementary Fig. 7A, B). We pooled plant tissue samples to achieve the 0.1-0. 5 g of dry plant tissue required to extract the total metabolites for metabolomics analysis. To ensure enough mass for extraction and appropriate replication for statistical analysis, we divided the 15 individual plants into 7 groups: 6 groups included 2 plants and the last group included 3 plants. After grouping, the tissues were ground into powders and metabolites were then extracted from an equivalent mass per pool. Though the soil samples did not require pooling to achieve sufficient mass, we performed the same grouping and extraction protocols to the soils to maintain consistency.

The plant tissues were lyophilized and ground into powders using an automated tissue homogenizer (SPEX SamplePrep, Metuchen, NJ). The plant metabolites were extracted from the powders using 80% methanol containing 1 µM telmisartan internal standard and normalized by tissue weight to achieve equal concentration as described in ref. ⁸. The soil metabolites were extracted following the same protocol with differences in the soil-to-solvent ratio (1:2, v/v) and incubation method (mixing on a laboratory tube rocker). Extracts were centrifuged at 4000 g for 20 min at room temperature to remove solids. The supernatant was completely dried using a SpeedVac vacuum concentrator (ThermoFisher, Waltham, MA) and reconstituted using 1/10 the original solvent volume.

Liquid chromatography–mass spectrometry (LC-MS) based untargeted metabolomics

The chromatographic separation and MS analysis for the switchgrass metabolites were performed using a reversed-phase, UPLC BEH C18, 2.1 mm × 150 mm, 1.7 μm column (Waters, Milford, MA) and an Electrospray Ionization – Quadrupole Time-of-Flight MS (ESI-QToF-MS, Waters). Mass spectra, under positive ionization, were acquired by data-independent acquisition (DIA)as described in Li, et al.⁸. The untargeted metabolomics data processing, including retention time (RT) alignment, lock mass correction, peak detection, adduct grouping and deconvolution, and metabolite annotation, were done using the Progenesis QI software package (v3.0, Waters) following the protocol in Li, et al.⁸. The identified analytical signals were defined by the RT and mass-to-charge ratio (m/z) information and referred to as the features. Measurement of each feature across the sample panel was filtered by interquartile range, log-transformed, and scaled for multivariate analyses using R Studio v3.1.1 and R package MetabolAnalyze (scaling function set to type “pareto”). The log-transformed and pareto scaled (normalized) abundance of features for soil were used for principal component analysis (PCA). A Partial Least-Squares Discriminant Analysis (PLS-DA) model was generated using the MetaboAnalyst 4.0 online tool platform⁶⁷ to assess each metabolite feature’s variable of importance (VIP) coefficient. The top 50 features with the highest VIP coefficient were used to visualize in a heatmap. For this, the features were log-transformed and scaled from 0-1 (by maximum) across all the samples. This was completed using R software and then visualized using the heatmap function in the R package ComplexHeatmap. Distance and clustering methods were set to Euclidean and Ward.D to generate hierarchical clustering for the heatmap. The data dependent acquisition (DDA) was done only for one pooled sample in order to obtain MS/MS spectra for annotating the top 50 PLS-DA features using CANOPUS (class assignment and ontology prediction using mass spectrometry)⁶⁸ machine learning function built in the SIRIUS 4 (https://bio.informatik.uni-jena.de/sirius/), a computational tool for systematic compound class annotation. The identification level was denoted for each annotated metabolite based on the criteria for metabolite identification as per Sumner, et al.⁶⁹.

16S rRNA and 16S rRNA gene sequence analysis

Sequence data were analyzed using QIIME2⁷⁰. All paired-end sequences with quality scores were compressed and denoised using the DADA2 plugin⁷¹. The denoising step dereplicated sequences, filtered chimeras, and merged paired-end reads. The truncation parameters to use with the DADA2 plugin were determined using FIGARO⁷². FIGARO analyzes error rates in a directory of FASTQ files to determine the optimal trimming parameters for sequencing pipelines that utilize DADA2. The truncation length was set to 123F and 162R for all the data, with minimum overlap set to 30 base pairs, which resulted in 93% merging success. All truncation was performed from the 3’ end for consistent final read lengths. The DNA and cDNA datasets were separately denoised. The resulting DNA and cDNA count tables were merged into a single QIIME2 artifact using the feature-table merge command.

Similarly, the DNA and cDNA representative sequences were merged into a single QIIME2 artifact using the feature-table merge-seqs command. The representative sequences from the combined count tables were clustered at 99% identity de-novo, and the clustered representative sequences were classified using SILVA v138⁷³ to generate the taxonomy file. Ninety-nine percent sequence identity was used to define OTUs to conservatively account for any potential amplification errors that may have occurred during the cDNA synthesis from the RNA. The resulting OTU (operational taxonomic unit) table and taxonomy files were exported to R for ecological analysis.

Designating the active community members

All downstream analyses were performed in R version 4.0. The R package decontam⁷⁴ was used to determine the number and identity of contaminants in the dataset (Supplementary Fig. 8A and 8B) and remove them using the prevalence method. Contaminating taxa, mitochondria, and chloroplast sequences were filtered from the datasets. Based on rarefaction curves, a subsampling depth of 15,000 reads per sample was selected (Supplementary Fig. 8C). After subsampling, 16S rRNA to rRNA gene ratios (hereafter, 16S rRNA:rRNA gene) were computed from the DNA and cDNA datasets as described in ref. ²¹. While we compared a few methods therein for this dataset (please see Supplemental Information for details), we ultimately chose the method that applied a 16S rRNA:rRNA gene ratio threshold > =1. The chosen method was statistically robust in overarching patterns of beta-diversity given the exclusion or inclusion of phantom taxa (taxa with detected cDNA but not DNA counts). For phantom taxa that were detected in greater than 5% of samples, the DNA counts = 0 were changed to DNA = = 1 (as in “method 2” in ref. ²¹) (see Supplemental Materials). All other phantom taxa were excluded. The DNA OTU table was filtered to include only sequence counts of active taxa in the samples determined to meet our ratio threshold. Consequently, while every DNA and cDNA sample of sequence counts was initially rarefied to 15,000 reads, each sample’s active community varied in their total reads (2000–6000). Relativized abundances were used for ecological statistics.

Microbiome data analysis for active community

The OTU table of the DNA counts of active taxa, taxonomy table, and metadata files was merged using the phyloseq package⁷⁵ in R. We used Bray-Curtis dissimilarity to determine beta diversity but also tested weighted UniFrac distance, and results were comparable. Permutational analysis of variance (PERMANOVA) was conducted using the adonis function in vegan package⁷⁶ to assess differences in community structure by treatment and interactions: drought treatment (watered or drought), drought severity/sampling (days 2-6), and plant treatment (planted or unplanted). For post hoc tests, pairwise comparisons between drought levels were computed using the pairwise.adonis2 function in vegan package. PERMANOVA tests using the adonis function in R vegan package were also done on metabolite feature abundances used for multivariate analysis to understand the effects of planting, drought, and sampling day factors. To understand drought dynamics, we analyzed the Bray Curtis similarity of the microbiome to the pre-drought samples over the covariate of time for planted and unplanted treatments in bean and switchgrass. This was visualized using a smoothed conditional means (geom_smooth function) with a linear model. This same approach was also used to assess general, relative fluctuations in rhizobiome size across samples, as proxied by DNA concentration.

Richness was the observed number of OTUs, using the estimate_richness function from 100 re-samples of the community. The normality of alpha diversity metrics was assessed using a Shapiro-Wilk test with a cut-off of W > 0.9 for normality assumptions. Since data were normally distributed, a parametric three-way ANOVA test was used to assess the main effects and interaction effects between factors on richness. Contrasts were set to “sum” before running the ANOVA model to ensure that type III ANOVA tests were valid. We used non-parametric Kruskal-Wallis tests (Kruskal.test in R) to calculate differences in class-level relative abundances between treatments.

We used the indicspecies package in R to determine indicator species associated with experimental conditions. We used the abundance-based counterpart of Pearson’s phi coefficient of association within the multipatt function⁷⁷. To correct the phi coefficient for unequal group sizes the “func” parameter within multipatt was set to “r.g.” P values were adjusted for false discovery rates. Indicator species were calculated at the OTU level for each treatment combination and the OTU and family level for the crop-specific indicators.

We created a heatmap visualization of the 50 most abundant and active OTUs in the bean dataset for indicator species. We used a maximum standardization approach using the decostand function in the vegan package in R. We distinguished taxa that changed relative abundance over the drought gradient from taxa that changed in their detection. For samples without detection of an OTU’s activity or DNA (e.g., DNA = 0, cDNA = 0), that OTU’s abundance was coded as NA. For samples that had no detection of an OTU’s activity but had detection of its DNA (e.g., DNA > 0, cDNA = 0), that OTU’s abundance was coded as 0. For samples that detected an OTU’s activity and DNA (e.g., DNA > 0, cDNA >0), that OTU’s abundance was the value of its DNA sequence count. We also included a special consideration of “phantom” OTUs (e.g., DNA = 0, cDNA>0), which were assigned a sequence count of 1.

Statistics and Reproducibility

In summary, we conducted a greenhouse experiment to understand active bacterial community dynamics during short-term drought. We used field soil with a legacy of growing two crops, switchgrass and common bean and conducted a drought treatment exposure over 6 days with incremental drought severity over time. We used watered controls and a pre-drought baseline for all comparisons. Appropriate multivariate statistical tests were conducted on the rhizobiome and metabolomics data, along with statistical tests for univariate data such as gravimetric moisture, shoot biomass and DNA biomass yield. Details are provided in the Methods sections above. All greenhouse experiments were conducted at least twice, which included a pilot study, which dictated the water addition levels appropriate to inflict drought in the respective crops without affecting the survival of those crops. Treatments were assigned randomly to crops during the greenhouse experiment. Even though no statistical method was used to predetermine sample size, we opted for 5 replicates as they accurately captured any differences observed in plant physiology during the experiment. A few data points were excluded from the rhizobiome analyses if they did not pass the rarefaction threshold.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.