Cover crop root exudates impact soil microbiome functional trajectories in agricultural soils

Abstract

Background

Cover cropping is an agricultural practice that uses secondary crops to support the growth of primary crops through various mechanisms including erosion control, weed suppression, nutrient management, and enhanced biodiversity. Cover crops may elicit some of these ecosystem services through chemical interactions with the soil microbiome via root exudation, or the release of plant metabolites from roots. Phytohormones are one metabolite type exuded by plants that activate the rhizosphere microbiome, yet managing this chemical interaction remains an untapped mechanism for optimizing plant-soil-microbiome interactions. Currently, there is limited understanding on the diversity of cover crop phytohormone root exudation patterns and our aim was to understand how phytochemical signals selectively enrich specific microbial taxa and functionalities in agricultural soils.

Results

Here, we link variability in cover crop root exudate composition to changes in soil microbiome functionality. Exudate chemical profiles from 4 cover crop species (Sorghum bicolor, Vicia villosa, Brassica napus, and Secale cereal) were used as the chemical inputs to decipher microbial responses. These distinct exudate profiles, along with a no exudate control, were amended to agricultural soil microcosms with microbial responses tracked over time using metabolomes and genome-resolved metatranscriptomes. Our findings illustrated microbial metabolic patterns were unique in response to cover crop exudate inputs over time, particularly by sorghum and cereal rye amended microcosms. In these microcosms, we identify novel microbial members (at the genera and family level) who produced IAA and GA₄ over time. Additionally, we identified cover crop exudates exclusively enriched for bacterial nitrite oxidizers, while control microcosms were discriminated for nitrogen transport, mineralization, and assimilation, highlighting distinct changes in microbial nitrogen cycling in response to chemical inputs.

Conclusions

We highlight that root exudate amendments alter microbial community function (i.e., N cycling) and microbial phytohormone metabolisms, particularly in response to root exudates isolated from cereal rye and sorghum plants. Additionally, we constructed a soil microbial genomic catalog of microorganisms responding to commonly used cover crops, a public resource for agriculturally relevant microbes. Many of our exudate-stimulated microorganisms are representatives from poorly characterized or novel taxa, revealing the yet to be discovered metabolic reservoir harbored in agricultural soils. Our findings emphasize the tractability of high-resolution multi-omics approaches to investigate processes relevant for agricultural soils, opening the possibility of targeting specific soil biogeochemical outcomes through biological precision agricultural practices that use cover crops and the microbiome as levers for enhanced crop production.

Video Abstract

Background

More than 153,000 US farms currently integrate cover cropping practices on over 15 million acres of land and cover crop adoption has increased in recent years [1]. Cover cropping is the incorporation of a secondary, unharvested crop grown typically in rotation with a primary cash crop. Crop diversification through cover cropping influences aboveground plant growth responses and belowground soil characteristics, promoting soil and primary crop health alike [2,3,4,5,6,7]. For instance, aboveground, cover crops act as a physical barrier against wind and water erosion, protect against weeds, and when utilized to replace periods of fallow, cover crops can contribute to the maintenance of soil functionalities that can bolster primary crop yields [7,8,9,10]. Belowground, cover crops stimulate microbial populations that catalyze beneficial biogeochemical processes [4, 11, 12]. The belowground impact of cover cropping is influenced by both cover crop root exudation during plant growth as well as the cover crop litter as it decomposes, but rarely are these influences detangled [13,14,15]. Despite the adoption of cover cropping as a regenerative agriculture practice, there is a paucity of data on cover crop root exudate chemical diversity and its impacts on soil biogeochemistry.

Metabolites released by plants as root exudates include diverse sugars, organic and amino acids, enzymes, and secondary plant metabolites. These primarily water-soluble and low-molecular-weight compounds influence the rhizosphere, altering the chemical landscape adjacent to plant roots and shaping the microbial recruitment and associations that are beneficial for plant hosts [16,17,18,19]. For example, exuded sugars and amino acids are carbon and nitrogen growth substrates for rhizosphere microorganisms, stimulating microbial activity near roots that provides benefits to the host plant of enhanced nutrient acquisition and pathogen suppression [20,21,22]. Moreover, secondary metabolites in exudates, like indoles and derivatives, contribute to microbially-mediated plant defense responses, growth, and plant–microbe signaling, playing a crucial role in shaping the overall ecology of the rhizosphere [20, 23, 24]. However, harnessing exudates for targeted biological stimulation of the rhizosphere is complicated by the diversity in both quantity and quality of root exudate compounds resulting from different crop species, and cultivars within species, which is altered in response to growth phase, abiotic and biotic stress, and many edaphic conditions [14, 15, 25,26,27].

We previously used metabolomics to characterize variability in root exudation across 19 commonly used cover crop species [25]. We assessed primary and secondary root exudates, and importantly, quantified low abundance phytohormones released as exudates. The cover crop species selected here are a subset of this prior study and representative of four functional plant groups (legumes, brassicas, cool-season grasses, warm-season grasses) commonly used in US agriculture [1]. We previously demonstrated these cover crop species had distinct root exudate chemical profiles, a finding we hypothesized could yield distinct soil microbial outcomes.

To decipher the microbial responses to cover crop root exudate treatments, we utilized laboratory-scale soil microcosms and stimulated microbial communities with daily amendments of pure root exudates from each over crop species (n = 4) over 6 days, tracking the responses over a 21-day experiment using metabolomics and genome-resolved metatranscriptomics (Fig. 1A, B). We found metabolomics was more sensitive than metagenomics for detecting differences in microbial responses between cover crops treatments and over the timeseries. However, functional gene expression and chemical behavior were still tightly coupled, especially for microbial phytohormone metabolism. We created a microbial genomic resource for exudate-stimulated microorganisms from agricultural soils here referred to as the Agricultural exudate Responsive Metagenomic (ARM) database. This integrated research is a first step towards the development of biological precision agriculture, linking cover cropping microbiome management practices to support valuable ecosystem resources.

Fig. 1

Experimental design. A Root exudates were collected from 4 hydroponically grown cover crops [25]. The soil microcosms were amended with soil from the agricultural research station at Colorado State University in biological quadruplets from 5 treatments: cereal rye (dark blue), sorghum (green), rapeseed (light purple), hairy vetch (orange), and a water amended soil control (gray). B In this study, microcosms were amended with root exudates from cover crops as shown in A for 6 days (denoted as the exudate addition phase) and responses were surveyed for 21 days (post-exudate addition phase). Metagenomes (n = 26), metatranscriptomes (n = 40), and metabolomes (n = 199) were collected at specific timepoints, indicated by circles, to profile microbial responses. The number of samples collected for each ‘omic measurement is indicated by the bar charts on the right. C Schematic summary of the data collected for the metagenome-assembled genome (MAG) database, which includes content from a prior study [25] and findings reported here. These MAGs were used to construct the Agricultural exudate-Responsive Metagenomic (ARM) database, with contributions from each experiment highlighted by bullet points

Materials and methods

Cover crop hydroponics and root exudate collections

Mature seeds from the four cover crops (Sorghum bicolor (sorghum); Vicia villosa (hairy vetch); Brassica napus (rapeseed); and Secale cereal (cereal rye)) were first surface sterilized for hydroponics in a sterile tissue culture hood. All hydroponic methods were completed using sterile techniques (i.e., surface sterilizing all equipment, materials, and gloves, working near a 100% EtOH flame, and autoclaving all materials before use) to ensure microbes were not introduced to the plants. First, ~ 10 seeds were placed into 15-mL conical tubes with one of two seed sterilization methods (A or B, below) chosen based on optimization of germination rate for each cover crop species (see Additional File 1 for germination rate and sterilization method for each seed). Method A consisted of 2 mL of Clorox bleach solution (3% sodium hypochlorite; 1:1 bleach to water) + 1 drop of Tween-20 (to reduce surface tension and improve sterilization). Seed sterilization method B consisted of only 2 mL of Clorox bleach solution (3% sodium hypochlorite; 1:1 bleach to water). Seeds suspended in sterilization solution were shaken by hand for 1 min. In a sterile tissue culture hood, the sterilization solution was removed by pipetting and the seeds were rinsed five times with sterile deionized (DI) water. Seeds were either (A) plated to agar plates with MS basal salt mixture, 8 seeds per plate, (2.16 mg Murashige and Skoog (MS) media in 500 mL sterile DI water; MP Biomedicals, Santa Ana, CA) germinated in the dark for 3 days and then individual seedlings were transferred to sterile tubes with 3 mL liquid MS media after 3 days; or (B) a single seed was placed directly in sterile growth vessel with 3 mL MS media. Growth vessels were either glass tubes which were used for monocotyledon species (cereal rye and sorghum) or Magenta boxes which were used for dicotyledon species (rapeseed and hairy vetch) based on preliminary experiments to optimize growth of each seedling. For hydroponic growth, all seedlings (n = 20) were incubated in a growth chamber for 14 days with photoperiod 16 h light/8 h night at 25 °C for 14 days (Additional File 1). The 14-day-old seedlings utilized did not introduce any microbes to the system as no bacterial growth was observed on the MS agar plates or growth containers.

Root exudates were collected by first rinsing the roots in sterile DI water 3 times and transferring to a new, sterile vessel filled with 5 mL sterile DI water for a 24-h root exudate collection period. The suspensions containing the root exudates were then filtered through a 0.2 μm filter membrane to remove root detritus. Enough root exudates were collected and pooled together to achieve a treatment addition of 1.5 mg of pure exudates per day (36 mg total needed for the experiment), an amount representative of the average phytochemical exudation over a 24-h period [25]. Filtered root exudate suspensions were frozen at − 20 °C in 200-mL sterile glass jars. To aliquot into daily treatments, frozen root exudate samples were lyophilized completely (~ 72 h), then weighed before resuspension in sterile HPLC-grade water (Thermo Fisher Scientific, Waltham, MA). Resuspended exudates were vortexed thoroughly to remove all residue from the bottle (1 min of vortexing along all edges of the bottle). A volume equivalent to 1.5 mg was aliquoted into new sterile glass jars for daily addition during the experiment and frozen at − 20 °C until use.

Cover crop root exudate microcosm experimental setup

Microcosms were established and sampled as previously described [28, 29] and metadata collection of all samples follows the recommended standards for agricultural microbiome research [30]. The soil (microbial inoculum) was collected from an agricultural field, previously planted to corn the year prior, at the Colorado State University Agricultural Research, Development, and Education Center near Fort Collins, CO on June 21, 2022 (GPS coordinates: 40.65493226483286, -104.99860614985835). The climate at the site is semi-arid, with 408 mm mean annual precipitation and a mean annual temperature of 10.2 °C (1981–2010 average, https://usclimatedata.com/). The soil is classified as an Aridic Haplustalf. Three 2-cm diameter soil core samples (at approximately 12 cm depth) were collected from 5 random locations within the plot (n = 5 cores) to achieve a representative bulk soil sample across the plot. The soil was stored on ice during transportation and at 4 °C until microcosm construction the following day. Next, 20 g of soil from each of the soil cores taken from the 5 locations within the plot was pooled, homogenized, and sieved (2 mm) to create a representative bulk soil repository used in the microcosms experiment.

Homogenized soil (4 g) and 45 mL of sterilized water were added to quadruplet sterile 100-mL glass bottles to construct each microcosm (n = 20). Microcosm slurries were vortexed and allowed to settle for 5 min, with day 0 samples then taken by removing 3 mL of soil slurry for parallel microbiome analysis. The microcosm slurry samples were centrifuged at 10,000 rcf for 10 min, with the supernatant used as a metabolome sample and the pellet used in the RNA/DNA extraction. After baseline sampling (day 0), 4 mL of exudate treatment (1.5 mg pure, dried, root exudates in 4 mL sterilized LC–MS-grade water) was added to each microcosm replicate and briefly vortexed. At this point, for each microcosm, the bottle caps were removed and replaced with a sterilized foam stopper (model number: Outus-17731) for the rest of the experiment to maintain oxic conditions and prevent colonization by contaminating microbes. Microcosms were incubated in a benchtop orbital shaker set at 200 rpm at 24 °C for 21 days. Each exudate treatment (1.5 mg exudate/4 mL water) (Fig. 1A) was added to microcosms on days 0, 1, 2, 3, 4, and 5 at approximately the same time every day. After day 5, no additional exudate treatments were applied, but microcosms were maintained until day 21, which afforded additional samples taken at days 7, 10, 15, and 21 (Fig. 1B). Samples were collected (3 mL of soil slurry removed and centrifuged to obtain pellet for DNA/RNA and supernatant for metabolome) at roughly the same time each day and collected with aseptic techniques to ensure no additional microbial influence was introduced. Samples for metabolomics and DNA/RNA were immediately frozen at − 80 °C until processing.

Metabolomics: sample preparation

Metabolomic analyses (targeted UPLC-MS/MS and nontargeted LC–MS/MS) were performed on all microcosm treatments (n = 5) across each of the ten timepoints (days 0, 1, 2, 3, 4, 5, 7, 10, 15, 21, n = 10) in quadruplicate (total n = 200 samples per analysis type). Once ready for analysis, frozen metabolomic samples were lyophilized and weighed to record the sample mass (tubes were pre-weighed prior to the experiment in order to calculate mass). We then added 1 mL sterile HPLC-grade water to each sample, and a volume equivalent to 0.25 mg of sample was transferred to a new, 2 mL glass autosampler vial and re-lyophilized to create two 0.25-mg subsamples for targeted and non-targeted metabolomics.

Metabolomics: targeted UHPLC-MS/MS for phytohormone analysis

To measure low abundance phytohormones, the first 0.25 mg subsample was extracted in 75 µL of a spiked methanol solution containing 100% methanol with 65.2 ng/mL abscisic acid-d6, 62.5 ng/mL salicylic acid-d6, and 90.0 ng/mL jasmonic acid-d5 (Sigma). After solvent addition, samples were placed on a shaker plate for 1 h at the highest speed setting, centrifuged at 3500 × g at 4 °C for 5 min, and transferred to glass inserts. A final centrifuge step at 3500 × g for 15 min at 4 °C was completed to ensure any precipitate was in the bulb of the vial insert. To measure instrument function a pooled QC was created by combining 5uL of each sample into a separate vial. This pooled QC was run after every 6 injections. Five microliters of the extract was injected onto an LX50 UHPLC System, equipped with an LX50 Precision Sampling Module (20-μL sample loop, partial loop injection mode) (PerkinElmer, Waltham, MA, USA). An ACQUITY UPLC T3 column (1 × 100 mm, 1.8 μM; Waters Corporation) was used for chromatographic separation. Mobile phase A consisted of LC–MS grade water with 0.1% formic acid and mobile phase B was 100% acetonitrile. The elution gradient was initially 0.1% B for 1 min, which was increased to 55.0% B at 12 min and further increased to 97.0% B at 15 min, then decreased to 0.1% B at 15.5 min. The column was reequilibrated for 4.5 min for a total run time of 20 min. The flow rate was set to 200 μL/min and the column temperature was maintained at 45 °C. Samples were held at 4 °C in the autosampler. Detection was performed on a QSight™ 420 triple quadrupole mass spectrometer (MS) operated in selected reaction monitoring (SRM) mode. SRM transitions for each compound were optimized through analysis of authentic standards (Additional File 2). The MS was operated with the ESI voltage 5000 V in positive mode and − 5000 V in negative mode. Nebulizer gas flow was set at 350 arbitrary units and drying gas was set to 120 arbitrary units. The source temperature was 315 °C and hot-surface induced desolvation temperature was 200 °C.

Metabolomics: nontargeted LC–MS and fastDDA analysis

To broadly measure primary and secondary metabolites, the second 0.25-mg subsamples were extracted in 80 μL of 20% MeOH, sonicated for 1 h at 20 °C, centrifuged at max speed for 15 min, and transferred to glass inserts for analysis. From each sample, 5 μL was aliquoted into a separate vial to be used as a pooled QC to monitor proper instrument function and to detect any analytical variation. From each sample within a treatment, 5 μL was aliquoted into a 2-mL vial to be used as a treatment pool for fastDDA.

One microliter of each sample was injected onto a Waters Acquity UPLC system. Separation was achieved using a Waters ACQUITY UPLC Premier BEH Amide 1.7 μm Column (2.1 × 100 mm), using a gradient from solvent A (10 mM ammonium formate in water with 0.125% formic acid) to solvent B (95:5 acetonitrile:water, 10 mM ammonium formate in water with 0.125% formic acid) and a flow rate of 0.5 mL/min. Column eluent was infused into a Waters Xevo G2-XS Q-TOF–MS with an electrospray source in positive ion, sensitivity mode, with data-dependent acquisition (DDA). For individual samples, the following parameters were used for the MS1 survey scan: 50–1200 m/z at 0.2 s per scan, switching to MS/MS after individual ion intensity rises above 10,000. MS/MS acquisition occurred at a scan rate of 0.2 s or at an accumulated TIC threshold of 100,000, with 1 MS/MS event per MS1 scan. No inclusion list was used, but dynamic peak exclusion was used with a 30-s and 100 ppm mass difference exclusion window. A pooled QC was run after every 7 normal sample injections.

In addition, pooled cereal rye, rapeseed, sorghum, control soil, control water, and hairy vetch QC samples were run separately with 10 replicate injections per pool in iterative exclusion mode using AutoCat_V1 processing. In this case, the following parameters were used for the MS1 survey scan: 50–1200 m/z mass range at 0.1 s per scan, switching to MS/MS after individual ion intensity rises above 5000. MS/MS acquisition occurred at a scan rate of 0.5 s, or after an accumulated TIC threshold of 100,000, with 5 MS/MS events per MS1 scan. For all experiments, collision energy for MS/MS was ramped from 15 to 30 V. Calibration was performed using sodium formate with 1 ppm mass accuracy. The capillary voltage was held at 700 V in positive mode. The source temperature was held at 150 °C and the nitrogen desolvation temperature at 450 °C with a desolvation flow rate of 1,000 L/hr. Lockspray reference mass was used to correct for drift, with 40-s intervals between scans, 0.1 s/scan and signal averaged over 3 scans. LeuEnk was used for mass correction, with a reference mass of 556.2771 m/z. For pooled iterative exclusion samples, a lockspray signal was collected but correction was not applied until post-processing. These latter data were processed using Waters Tynebridge to produce mass-calibrated mzML files compatible with GNPS server. The column and samples were held at 30 °C and 6 °C, respectively.

Metabolomics: data analysis

For the analysis of low abundance phytohormones, Simplicity 3Q software (Version 3.0.2, PerkinElmer, Waltham, MA) was used to detect and integrate peak areas and to calculate linear regression of analytical standards used for quantification. Each peak was normalized to an appropriate internal standard. The corresponding linear regression equation was used for quantification (ng/mL) for each analyte, which was then adjusted for the precise volume of the sample (ng/mL). The limit of detection was calculated as 3 times the standard deviation of the blank divided by the slope of the calibration curve. Likewise, the limit of quantitation was calculated as 10 times the standard deviation of the blank divided by the slope of the calibration curve.

For non-targeted LC–MS data, .mzML files were processed through the following workflow: (1) XCMS software was used for preprocessing to identify molecular features [31]; (2) features were further normalized to total ion current; (3) the package RAMClustR (v1.3.1) [32] was used for clustering features into spectra; (4) the package InterpretMSSpectrum (v1.4.5) [33] was used to produce clusters which were then input into MSFinder [34] for spectral matching, formula inference, and tentative structure assignment.

For LC–MS/MS fastDDA data, .mzML files for each iterative exclusion across all treatments (n = 50) were uploaded to GNPS for molecular networking and annotation. A molecular network was created using the online workflow (https://ccms-ucsd.github.io/GNPSDocumentation/) on the GNPS website (http://gnps.ucsd.edu). The data was filtered by removing all MS/MS fragment ions within ± 17 Da of the precursor m/z. The precursor ion mass tolerance was set to 0.5 Da and a MS/MS fragment ion tolerance of 0.02 Da. A network was then created where edges were filtered to have a cosine score above 0.5 and more than 5 matched peaks. Further, edges between two nodes were kept in the network if and only if each of the nodes appeared in each other’s respective top 20 most similar nodes. Finally, the maximum size of a molecular family was set to 100, and the lowest scoring edges were removed from molecular families until the molecular family size was below this threshold. The spectra in the network were then searched against GNPS’ spectral libraries. The library spectra were filtered in the same manner as the input data. All matches kept between network spectra and library spectra were required to have a score above 0.5 and at least 4 matched peaks.

Annotation of compounds to the level 2 classification was completed through manual matching between the MS1 data from all samples and the MS2 data collected from the pooled QCs. Specifically, annotation results from GNPS molecular networking were exported and precursor masses were matched to MS1 spectra within an 11 s retention time window (to account for analytical drift). Quantitative data (untargeted LC–MS and targeted UHPLC-MS/MS for phytohormone analysis) were z-scored and combined for the final quantitative dataset (Additional File 2).

Metabolomics: chemical taxonomic assignments

For non-targeted LC–MS/MS molecular features, MSFinder InChIKey results were exported. For targeted LC–MS/MS, InChIKeys were gathered from PubChem. All InChIKeys were uploaded to ClassyFire Batch (https://cfb.fiehnlab.ucdavis.edu/) for batch annotation of chemical taxonomy [35]. Metabolites quantified via targeted UHPLC-MS/MS were assigned a level 1 annotation. A level 2 annotation was assigned to molecular features in which a feature in MS1 was also present in MS2. A level 3 annotation was assigned to molecular features that did not have corresponding MS1 and MS2 data, but could be assigned chemical taxonomy according to ClassyFire.

Metagenomics: sample preparation

For days 3, 5, 7, 10, and 21 (n = 5 timepoints), we obtained samples for metagenomics analysis from each treatment (n = 5) from a single microcosm (total n = 25 metagenomes; Fig. 1B). For day 0, the control microcosm was used for a total of 26 metagenomes sequenced. For these samples, genomic DNA was prepared for metagenomic sequencing using the ZymoBiomics DNA/RNA Miniprep kit and the Tecan Ovation Ultralow System V2 library prep kit and was sequenced at the University of Colorado Anschutz Sequencing Core on the Illumina NovaSeq 6000 with 2 × 150 bp chemistry at 20 million read pairs per sample. This resulted in an average of 9.12 Gbp per sample, for a total of 237.15 Gbp of metagenomic sequencing.

Metagenomics: assembly, binning, and annotation

Fastq files were trimmed using Sickle (v1.33) [36], and trimmed reads were assembled using MEGAHIT (v1.2.9) [37]. To maximize genome recovery, we performed 3 assemblies: (1) individual assembly on each fastq (n = 26), (2) co-assembly by combining reads from each treatment fastq metagenomes to increase assembly coverage (n = 5), and (3) a single iterative assembly was performed on reads that did not map to binned assembled scaffolds ≥ 2.5 kb at 97% identity on all metagenomes. Information for metagenome statistics, including assembly information, is found in Additional File 3. For each assembly, scaffolds ≥ 2.5 kb were binned using MetaBAT2 (v2.12.1) [38], and MAG quality was assessed using checkM (v1.1.2) [39]. This resulted in 737 total MAGs, but only MAGs ≥ 50% complete and ≤ 10% contaminated were kept in the database. This resulted in 441 medium- and high-quality MAGs which were then dereplicated at 99% identity using dRep (v2.6.2) [40] into 326 MAGs. MAG taxonomy was assigned using GTDB-tk (v2.3.0 Release 08-RS214) [41]. MAGs and assemblies were annotated using DRAM (v1.4.4) [42]. CAZymes were inferred from the DRAM hits. See Additional File 4 for DRAM distillate. See Zenodo for raw annotations (https://doiorg.publicaciones.saludcastillayleon.es/10.5281/zenodo.11390843) and for the ARM database fasta (https://doiorg.publicaciones.saludcastillayleon.es/10.5281/zenodo.11390808).

To quantify MAG relative abundance in each metagenome, trimmed metagenomic reads were mapped to the dereplicated MAG database using Bowtie2 (v2.4.5) [43], output sam files were converted to sorted bam files using samtools (v1.9) [44] and entries in the bam file were filtered to 95% identity using the reformat.sh script from the BBTools package (v38.90) at minidfilter = 95 [45]. We had two requirements for a MAG to be present in a sample: first, we required reads to map to at least 75% of a MAG in a given sample, and second, the MAG had to have at least 3X coverage in that sample. To determine MAGs that had reads mapped to at least 75% of the MAG, we used coverM (v0.6.0) [46] in genome mode to output MAGs that passed this threshold (–min-covered-fraction 75). To obtain MAG coverage, we used coverM (v0.6.0) in genome mode to output reads_per_base (reads mapped/genome length), and from this calculated MAG coverage as reads_per_base × 151 bp. A bin was present in a treatment or in control if it was found in any of the timepoints.

Metagenomics: gene annotation curation

Genes for cysteine desulfhydrase are known to misannotate as 1-aminocyclopropane carboxylic acid (ACC) deaminase (acdS) [47,48,49]. To curate these, we took genes that were annotated as “K01505” from DRAM and aligned them against a known and validated acdS (UniProt: Q5PWZ8) using MUSCLE (default parameters) [50]. We then looked for the key acdS active site residues E295 and L322 and confirmed the presence of K51, Y268, and Y294 (Additional File 3). Genes for gibberellic acid (GA) biosynthesis in KEGG were annotated for fungi and plants and do not annotate model microorganisms that have GA biosynthetic capacity. Instead, we used a BLAST analysis with individual protein sequences of each gene as the query to survey MAGs containing genes located in the GA operon [51]. MAGs that contained sequence matches were considered if the BLAST E value was 10⁻¹⁰ or less and a bitscore > 150. Because CYP115 is the key enzyme converting GA₉ to GA₄, this was used as the determinate for a MAG that produces bioactive GA₄. Preliminary nitrogen oxidation and reduction DRAM annotations were confirmed using a combination of phylogeny and active site analysis where appropriate.

Metatranscriptomics: sample preparation

For days 0, 5, and 21, we obtained triplicate metatranscriptomes for each treatment/timepoint (n = 45). However, library preparation failed for replicate 2 of each treatment at day 5 leading to n = 40 metatranscriptomes (Fig. 1B). RNA was prepared for sequencing using the ZymoBiomics DNA/RNA Miniprep kit and cleaned using ZymoBiomics RNA Clean & Concentrate Kit and sent to the Department of Energy’s Joint Genome Institute for sequencing on the Illumina NovaSeq 6000 with 2 × 151 bp chemistry at a target depth of 100 M reads per sample. rRNA was depleted from an input of 10 ng of total RNA using QIAseq FastSelect™ – 5S/16S/23S, rRNA Plant, and rRNA Yeast Kits (Qiagen). Using TruSeq stranded mRNA kit (Illumina), the 300 bp—400 bp heat fragmented RNA was reverse transcribed to create the first strand of cDNA with random hexamers and SuperScript™ II Reverse Transcriptase (Thermo Fisher Scientific) followed by second strand synthesis. The double-stranded cDNA fragments were treated with A-tailing, ligation with JGI’s unique dual-indexed adapters, and enriched using 13 cycles of PCR.

Metatranscriptomics: mapping and analysis

Fastq files were trimmed, and adapters were removed using bbduk (v38.90) [45] with the parameters ktrim = r, k = 23, mink = 11, hdist = 1, and filtered using rqcfilter2 (v38.90). Trimmed, filtered reads were then mapped via Bowtie2 (v2.4.5) [43] to the MAG database (dereplicated to 99% ID). Sam files were transformed to bam files using samtools, (v1.9) filtered to 95% id using reformat.sh and name sorted using samtools [44]. Transcripts were counted for each gene with feature-counts (v1.5.3) [52]. Genes were considered if they had counts greater than or equal to 3 in at least 3 samples. The counts were then transformed to geTMM (gene length corrected trimmed mean of M-values) in R (v4.3.1) using edgeR (v3.42.4) package [53]. MAGs were considered present in a sample if they had greater than or equal to 10 genes in a sample (Additional File 5). We also analyzed the metatranscriptome at the function level, where KEGG, MEROPS, and dbCAN2 assigned annotations were used to infer gene function. Gene geTMMs within a function were summed to obtain the overall functional activity of the microbial community within a metatranscriptome (Additional File 5).

Statistics

Peak area files for each metabolomic analysis were combined into one file (Additional File 2) and normalized by z-scoring. Univariate statistics were conducted in GraphPad Prism (v9.4.1). Line graphs of targeted and nontargeted data were visualized in Prism using absolute concentrations and TIC normalized values, respectively. Two-way ANOVAs (multiple comparisons were corrected using the Dunnett’s multiple comparisons test) were calculated using log₁₀ transformed data to satisfy the assumption of data normality in order to determine if a compound was significantly different between a treatment and the control at each timepoint. Multivariate metabolomic statistics were conducted in SIMCA (v17.0.1) to generate partial least square discriminate analysis (PLS-DA) models and gather variable importance in projection (VIP) scores (here, we used a cutoff of VIP > 1.8). PLS-DA was performed using annotated and unannotated metabolites using z-scored and UV-scaled data in SIMCA. The list of scores and loadings coordinates were then plotted in Prism. The significance between metabolic changes across time and treatment was quantified using adonis2 from the vegan package (v2.6–4) in R (v4.3.1) with a PERMANOVA model [54]. Pairwise comparisons were completed using the pairwise_permanova wrapper function using adonis2 from the vegan package (v2.6–4) in R (v4.3.1) [54, 55].

Metatranscriptomic non-metric multidimensional scaling (NMDS) plots were generated in R (v4.3.1) to estimate beta diversity across treatments. We utilized the Bray–Curtis dissimilarity matrix visualized by NMDS in R with the ggplot2 (v.3.4.4) package with stress of the non-parametric fit for the goodness of fit for the model [56]. NMDS scores were exported from R and imported for visualization into Prism. The significance of compositional differences across treatments and the interaction of treatment and time was quantified using a multi-response permutation procedure, and the betadisper commands from the vegan package (v2.6–4) with an ANOVA model in R. Metatranscriptomic hierarchical clustering was completed in R using the hclust function from the stats package (v4.3.1). Visualization of the clustering used the as.dendrogram function in the stats package [55].

Discriminant metatranscriptome MAGs and functions were identified using MaAsLin2 (v1.14.1) [57] using default parameters with various tests considered: (1) day 5, MAG-level expression, control versus each individual treatment; (2) day 5, MAG-level expression, control versus all exudate treatments; (3) day 5, function-level expression, control versus all exudate treatments. All MaAsLin2 analyses were ran in R using the treatment (control or exudate) as a fixed effect. See Additional File 6.

Results

Agricultural exudate-Responsive Metagenomic (ARM) database expands microbiome knowledge in Colorado soils

To understand how microbial genomic content is influenced in response to exudate stimulation, we established the Agricultural exudate-Responsive Metagenomic (ARM) microbial genomic catalog (see the “Availability of data and materials” section). In total, we collected 41 metagenomes from two exudate stimulation soil microcosm experiments (Fig. 1C), resulting in 621 Gbp of total metagenomic sequencing. The prior proof of concept study evaluated sorghum genotype exudate profiles using microcosms operated similarly except exudate profiles guided the formulation of synthetic root exudates [28], while here we directly added the cover crop-produced exudates to the soil microcosms. From both experiments we reconstructed 441 medium- and high-quality metagenome-assembled genomes (MAGs) that were dereplicated at 99% identity into 326 MAG clusters, representative of distinct microorganisms (Fig. 2A). ARM contains MAGs assigned to 21 bacterial phyla and 1 archaeal phylum, with the database composed mainly of members of Pseudomonadota (29.1%), Bacteroidota (14.1%), Actinomycetota (13.1%), and Bacillota_A (12.2%). The MAGs of ARM span 43 different genera (Fig. 1; Additional File 7 Fig. 1).

Fig. 2

Curation of an exudate-responsive metagenome-assembled genome (MAG) database. A The taxonomy of the dereplicated MAGs is shown by sequentially colored rings ordered from domain (D, inner ring; gray = bacteria, purple = archaea), phylum (P), class (C), order (O), family (F), genus (G), to species (S, outer ring) assignment. Ring color corresponds to phylum, with the taxonomic assignment denoted in the legend to the left. Gaps at each level represent MAGs that were unclassified at that level of taxonomy (according to GTDB v2.3.0 08-RS214). B Stacked bar chart shows novelty distributions of ARM MAGs. Bars indicate the number of dereplicated MAGs recovered that represent unassigned families and genera as well as MAGs assigned only an alphanumeric identifier. Here, novelty is defined as the first unnamed level (i.e., unnamed family or genus) and the level where an alphanumeric identifier is used (i.e., alphanumeric family). Coloring corresponds to MAG phylum

Highlighting the genomic novelty of these soils, the ARM database contains MAGs that represent previously unidentified (i.e., lacking a taxonomic assignment) families (n = 10 MAGs, 3%) and genera (n = 39 MAGs, 11.9%) in the Genome Taxonomy Database (GTDB) (Fig. 2B, Additional File 3). Additionally, a large proportion of our MAG database belonged to lineages only recognized by alphanumeric identifiers (e.g., Draft Genome Sequence) in GTDB at the class (n = 7), order (n = 46), and family (n = 104) levels (Fig. 2B). Together these findings emphasize the phylogenetic novelty yet to be genomically captured from agricultural soils. Ultimately, accessible ARM resources provide new genomic information for ecologically relevant taxa, with the goal towards enabling taxonomic analyses and metabolic reconstruction of microorganisms from agricultural soils.

Cover crop exudates result in distinct microbial chemical landscapes

We first examined how soil microbial metabolite pools changed in response to the addition of cover crop exudates over time. Combining non-targeted and targeted LC–MS/MS data resulted in the detection of 641 molecular features to resolve community metabolic changes across the enrichment timeseries in each treatment (Additional File 2; Additional File 7 Fig. 2). The metabolite data showed restructuring by cover crop exudate treatment and time (Fig. 3 and Additional File 7 Fig 3 and Additional File 7 Table 1). Additionally, the cereal rye and hairy vetch treatments exhibited significantly different metabolomes from the control at each matched timepoint following the first day of exudate addition, indicating adding exudates not only rapidly influenced microbial outputs during the time they were added, but had an impact on restructuring of metabolic outputs after the amendment period. In the sorghum and rapeseed exudate treatment, the exudate amendments only differed from control on days 3 and 5, respectively, of the amendment period, yet despite these slight chemical changes early on, they did result in a temporally significant trajectory different from the control later in the timeseries. These results suggest that that cereal rye and hairy vetch root exudates may have a more distinct chemical influence over the soil microbiome than sorghum and rapeseed treatments.

Fig. 3

Cover crop exudate treatments influence soil microbial metabolomes. Two-dimensional scores plot of partial least squares discriminate analysis (PLS-DA) (R₂X = 0.101, R₂Y = 0.0597, Q₂ = 0.784, PERMANOVA, p < 0.05) between all treatment metabolomes. Each point represents a metabolome, with colors representing time. Shapes denote treatment: cereal rye (circles), hairy vetch (triangles), sorghum (diamonds), octagons (rapeseed) and control (square) metabolomes. Ellipses represent 95% confidence intervals and are colored by treatment. Corresponding loadings biplot can be found in Additional File 7 Fig 2

The 641 detected metabolites were then classified into 14 superclasses across all samples: lipids and lipid-like molecules (n = 179); organoheterocyclic compounds (n = 106); organic acids and derivatives (n = 97); phenylpropanoids and polyketides (n = 64); organic oxygen compounds (n = 52); benzenoids (n = 38); nucleosides, nucleotides, and analogs (n = 20); alkaloids and derivatives (n = 18); organic nitrogen compounds (n = 15); lignans, neolignans. and related compounds (n = 7); organosulfur compounds (n = 3); organic 1,3-dipolar compounds (n = 1); organic polymers (n = 1) (Additional File 7 Fig. 4; Additional File 2). We found a significant difference between treatment and control metabolomes at this metabolite classification level (Additional File 7 Fig. 5 and Table 2). Specifically, all exudate-treated microcosms were significantly enriched for lipids, benzenoids, alkaloids and derivatives, and organoheterocyclic compounds. Water-amended soil controls were enriched for organic acids and derivatives, nucleosides, nucleotides, analogs, and organic polymers (Additional File 7 Table 3). Interestingly, cereal rye was significantly different from the control and other exudate treatments being enriched for metabolites from the phenylpropanoids and polyketides superclass (Additional File 7 Table 3). Taken together our broad metabolite profiling indicates that cover crop exudates can result in distinct chemical trajectories, possibly through altered microbial metabolism.

Fig. 4

Overall expressed gene content, but not transcriptionally active genera, are altered by cover crop exudates. A Average metatranscriptome profiles of the MAGs from control and exudate microcosms. MAG metatranscriptome abundance was summed at the genus level and then averaged across replicates. Sorted by phyla, colors correspond to MAG genus and MAGs with an abundance less than 1% are shown in black. B Nitrogen cycling functions differed between control and all cover crop exudate metatranscriptomes at day 5. Arrow thickness corresponds to significant MaAsLin2 coefficients (roughly effect size), and arrow color indicates the treatment with which the feature associated (gray for the control, purple for the exudates). Dashed black lines correspond to functions that were detected in the metatranscriptomes but were not discriminant. For organisms discriminate to each function, see Additional File 7 Fig. 12. For plots of nitrogen function expression, see Additional File 7, Fig. 11. To see all discriminant functions and their statistics, see Additional File 6

Fig. 5

Cover crop treatment influences microbial phytohormone biosynthesis across time. A Partial least squares discriminate analysis (PLS-DA) (R₂X = 0.101, R₂Y = 0.0597) biplot shows the relative contribution of each phytohormone abundance found within a treatment metabolome. Colored stars indicate a treatment and colored circles indicate a phytohormone. Colored circles with text labels indicate phytohormones discussed in the text. B–D Line graphs show temporal dynamics of 4 selected phytohormones colored by treatment. Circles represent the average concentration at that timepoint, and error bars represent one standard deviation. B Indole-3-acetic acid (IAA), C 1-aminocyopropane carboxylic acid (ACC), D gibberellic acid 4 (GA₄), and E indole-3-butyric acid (IBA)

Microbial community function is altered by exudate addition

Next, we wanted to examine the microbial community response to exudates at the genome membership and gene expression levels. We used the ARM database to resolve metatranscriptome expression profiles across the enrichment timeseries (Additional File 5). Unlike the metabolomes that had a strong individual cover crop effect, the overall composition of metatranscriptionally active MAGs did not change by cover crop treatment, but rather time (Additional File 7 Fig. 6 and Tables 4–5). Yet we did observe that the active genera changed over the course of the experiment. For example, genera in the family Nitrososphaeraceae (phylum Thermoproteota, previously Thaumarchaeota) contributed roughly 50% of exudate-addition phase metatranscriptome expression (Fig. 4A, green shading). But by day 21 of enrichment, control and cover crop microcosm metatranscriptomes were dominated by MAGs from the Bacteroidota phyla, including unclassified genera in the Crocinitomicaceae (genus M2408) and JADKCL01 (genus JADKCL01) families (Fig. 4A). Comparing these membership and metabolite findings, the chemical data showed much stronger restructuring by treatment and time than the genus level active members.

Fig. 6

Gene potential and expression across phytohormone metabolic routes. A Heatmap shows the expression (in red) or the gene potential (in pink) for each gene in ACC degradation (acdS; purple), GA₄ biosynthesis (CYP115; blue), or IAA biosynthesis (shades of green correspond to different IAA pathways; darkest green for the indole-3-pyruvate (IPA) pathway, dark green for the indole-3-acetamide (IAM) pathway, green for the tryptamine (TAM) pathway, and light green for the indole-3-acetonitrile (IAN) pathway, the final lightest shade of green represent the terminal oxidation of indole acetaldehyde (IAAld) to IAA)). MAGs listed are those only detected in exudate amended microcosms and not controls as indicated in the main text. When the genus name was undescribed, the family name was used. Colored boxes to the right indicate MAG phylum with taxonomy in the legend to the right. Colored boxes at the bottom indicate the pathway. B ACC deamination (acdS) yields ammonium and α-ketobutyric acid, a compound that can be further transformed into propionyl-CoA or isoleucine

To better identify specific organismal responses to individual cover crop treatments we next focused on the genome (not genus) resolved responses within the microbial community. During exudate amendment, only a single MAG from an undescribed species in the Actinomycetota increased in activity in response to hairy vetch amendment (JAJTCL01 sp021323255), a single MAG from a novel genus in the Patescibacteria (family UBA1547) was enriched in response to cereal rye exudates, while a single MAG from an undescribed species in the Pseudomonadota (genus Hydrogenophaga) increased in response to cereal rye and sorghum exudates (Additional File 6, Additional File 7 Fig. 7). Next, we identified MAGs that broadly responded to the cover crop exudates compared to the control. We found 23 MAGs were enriched in the control metatranscriptomes and 9 MAGs were enriched across the cover crop exudate treatment metatranscriptomes (Additional File 7 Figs. 8 and 9. The latter group included three MAGs in the Nitrospiraceae (including one in the Nitrospira_C and classified as Nitrospira_A sp021323395), two MAGs in the Gammaproteobacterial genus Hydrogenophaga, one MAG classified as Nitrobacter vulgaris, and MAGs representing novel species of Arthrobacter and Lacipirellula. Collectively, this suggests cover crops produced root exudate compounds that stimulated specific microorganisms in the broader community.

Fig. 7

Indole-3-acetic acid and gibberellic acid 4 biosynthesis routes and associated MAGs. A IAA biosynthetic pathways. Four routes for IAA biosynthesis were detected by metagenomics and/or metatranscriptomics where detection is indicated by arrow type: hollow arrows indicate the gene was not detected, dotted arrows indicate the gene was only encoded (metaG only), solid arrows indicate the gene was encoded and expressed. Enzyme names are in gray next to the corresponding reaction arrow and given an enzyme number (“enzyme #”). Pathways are given a corresponding ID # in circles. Colors correspond to the IAA pathway. B Heatmaps show the z-scored geTMM value for each summed gene abundance across all MAGs expressing the gene (corresponding to the enzyme number in A) within a treatment and timepoint. C Downstream GA₄ biosynthetic route for production of bioactive GA₄. D The GA operon highlighting a potentially new bacterium of the Thermomicrobiales, which expressed 5 of 9 GA biosynthetic proteins required to produce GA₄, including the final enzyme (CYP115) converting GA₉ to GA₄. Gray arrows indicate the gene was neither expressed nor encoded in the Thermomicrobiales MAG, blue arrows show an encoded gene in the Thermomicrobiales MAG, and a blue arrow with a black border shows a gene was encoded and expressed in the Thermomicrobiales MAG at day 21. E Averaged geTMM abundance of MAGs producing CYP115 in cereal rye and control microcosms. F Expression of CYP115 from the potentially novel Thermomicrobiales MAG in cereal rye and control microcosms

Next, the impact of exudate addition on gene expression was evaluated by aggregating genes at the functional level (by annotation protein type). We identified 114 microbial functions that were enriched in the control metatranscriptome, and 145 functions that were enriched in the cover crop exudate amended microcosms, at day 5 (Additional File 7 Fig. 10). We observed distinct nitrogen cycling gene expression between control and exudate metatranscriptomes at day 5 (Fig. 4B; Additional File 7 Fig. 11). Functions associated with nitrogen transport, mineralization, and assimilation were discriminant to the control microcosm metatranscriptomes. In support of the microorganisms in the control microcosms being nitrogen limited, we confirmed the exudate metabolites contained organic nitrogen in the form of urea and several amino and organic acids while the controls received only water amendments. Additionally, ammonia oxidation gene expression (amoB) was more enriched in the control, with contributions from four genera of Thaumarchaeota (JAJNBK01, Nitrososphaera, TH1177, TH5893; Additional File 6; Additional File 7 Fig. 12). Nitrite reduction (nirK) was also enriched in the control, with an average of 97% control gene expression from the same four Thaumarchaeota genera. The increased expression of amoB further supports the idea that the control soil microcosms are nitrogen limited, as genera from archaeal ammonia oxidizers were detected and active in all cover crop and control treatments (Fig. 4A), but this increased gene expression in controls may be due to the low ammonia environment [58].

The only nitrogen function enriched in the exudate microcosms was nitrite oxidation (nxrB), or the second step of nitrification resulting in nitrite being oxidized to nitrate. Gene expression was derived from Nitrospira_C MAGs and a MAG from an undescribed genus (WHTF01) in the Binatia class (Additional File 6; Additional File 7 Fig. 12). Concomitant with this gene result, as mentioned above, a third of the discriminant cover crop transcriptionally active MAGs were bacterial nitrite oxidizers (Nitrospira_C, Nitrobacter). The enrichment of bacterial nitrite oxidizers was unique to the cover crop treatments, suggesting these nitrifiers selectively responded to exudate components. Despite similar redox conditions between exudate and control reactors, in the control soils gene expression data suggest that nitrite is not oxidized but reduced in a closed loop (Fig. 4B). Instead in the controls, nitrite reduction was mediated by Thaumarchaeota (nirK) possibly to detoxify ammonification produced nitrite due to the absence of active nitrite oxidizers or was assimilated to ammonia (nirBD). Finally, 8 genera contributed to nirD expression in the control, but 83% of the expression comes from 2 MAGs from the Cellvibrio and Pseudomonas_E genera. An additional 7 genera contribute to nirB expression in the control, with 38% of the expression from the same Pseudomonas_E MAG (Additional File 7 Fig. 12). These findings reflect an increase in control microorganisms to utilize specific nitrogen cycling functions due to the lack of available organic nitrogen (Fig. 4B). While we acknowledge limitations in translating these laboratory findings to the field scale, our multi-omics results do highlight how small-scale nutrient landscape changes caused by exudation could alter microbially expressed metabolic regimes (e.g., nitrification versus nitrite reduction) to result in disparate ecosystem manifested outcomes.

Metabolite evidence for cover crop and microbially produced phytohormones

Given the recognized importance of microbial phytohormones in plant growth promotion, we next considered whether cover crop exudation could stimulate soil microbial phytohormone metabolism. We performed targeted metabolomics to quantify phytohormones in the microcosms, identifying 8 phytohormones that fluctuated over time: indole-3-acetic acid (IAA), indole-3-butyric acid (IBA), gibberellic acid 4 (GA₄), 1-aminocyclopropane carboxylic acid (ACC), salicylic acid (SA), methyl salicylate (mSA), benzoic acid (BA), and jasmonic acid (JA) (Fig. 5). Of these phytohormones, IAA, IBA, BA, GA₄, and ACC had the strongest effect across the treatments (Fig. 5A; Additional File 7 Tables 6–9). Notably, only IAA and ACC were detected in cover crop exudate additions and were thus not detected in the water controls (Fig. 5B, C). Of the other phytohormones, IBA and GA₄ were stimulated to be produced by the microbial community in response to added root exudates (Fig. 5D–E), while BA and mSA production was unique to the non-exudate stimulated controls (Additional File 7 Fig. 13). Our findings show that cover crop roots exude unique phytohormone profiles, but also that the soil microbiome likely encodes diversity of phytohormone producing and consuming metabolisms.

Of the different phytohormones, we were particularly interested in ACC, IAA, IBA, and GA₄ dynamics because they had a clear effect on the microbial metabolic responses from different cover crop treatments (Fig. 5A). Our exudate metabolomics data confirmed ACC was detected as a component of the endogenous root exudates across all cover crops. As a result, ACC increased during the exudate addition phase (up to day 5). After day 7, ACC abundance stabilized (e.g., hairy vetch) or decreased (e.g., rapeseed, sorghum, cereal rye), the latter indicating a possible microbial consumption (Fig. 5C). Like ACC, IAA was a product of cover crop direct addition in the hairy vetch but was not a component of other exudate treatments. However, the IAA accumulation signal was not as clear as ACC during the exudate addition phase (up to day 5) (Fig. 5B). After exudate stimulation stopped, IAA generally decreased in 3 cover crop treatments (hairy vetch, cereal rye, rapeseed), except for sorghum where IAA increased overtime (Fig. 5B; Additional File 7 Table 10). Similarly, IBA, an indole related to IAA, accumulated significantly only in the cereal rye at the post exudation phase (Fig. 5E; Additional File 7 Table 6). These latter two findings suggest that the exudate chemistry stimulates members of the microbiome for indole production.

Based on these interesting indole cover crop responses, we mined our untargeted LC–MS data for additional indoles, identifying 2 unknown compounds that shared structural similarity to IAA (see Methods). Unknown indoles were from the chemical subclasses (i) indoles and (ii) indolyl carboxylic acids and derivatives and could be further categorized as a 3-alkylindole compound and an IAA derivative, respectively (Additional File 7 Table 11). The unknown IAA derivative was significantly more abundant in the rapeseed metabolomes compared to the control at each timepoint following day 1 (Additional File 7 Fig. 14A and Table 12). Exhibiting a broader response, the unknown 3-alkylindole was significantly enriched at each timepoint following the first day of exudate addition in rapeseed, sorghum, and hairy vetch compared to the control (Additional File 7 Fig. 14B and Table 13).

While there are more than 100 gibberellic acids (GA) known, only four (GA₁, GA₃, GA₄, and GA₇) have characterized biological activity [59]. Of these, GA₄, a C20-gibberellin, was detected exclusively, and increasingly over time, in the cereal rye-treated microcosms (Fig. 5D; Additional File 7 Table 7). Further supporting broader microbial GA production, we found two GA features that shared structural similarity with GA₄ that were further characterized as type C20- and C19-GA (Additional File 7 Table 11). Given that GA compounds were not detected in the exudates, and they collectively increased after exudate stimulation in the cereal rye (Additional File 7 Fig. 15A–B and Tables 14–15), it is likely that cereal rye exudates contained compounds that stimulated microbial production of these GA. These findings show the power of combining untargeted metabolomics with targeted phytohormone quantification to enable tractability of microbial phytohormone biosynthesis. Our metabolite analyses highlight how much yet to be discovered biochemistry is likely residing in agricultural soils.

Metatranscriptome evidence for a microbial role in phytohormone transformations

Based on the possible ACC microbial consumption observed in our metabolite dynamics (Fig. 5C) and the critical role this compound plays as a recruitment cue for plant growth-promoting bacteria [60], we mined ARM MAGs for the ACC deaminase (acdS) gene. Many bacteria use acdS to metabolize ACC via deamination, curbing the abiotic stress-induced ethylene production and its adverse effect on plants [60]. In total, 66 genes from 63 MAGs were annotated as acdS (Additional File 3). As ACC deaminase genes are commonly misannotated as D-cysteine sulfhydrase [47, 48], we curated these putative hits for key active site residues to identify 13 high-confidence acdS genes from distinct MAGs in the Actinomycetota (n = 9) and the Pseudomonadota (formerly Proteobacteriota, n = 4) (Fig. 6A).

Two of these Actinomycetota MAGs (Pseudonocardia, Ornithinibacter) expressed ACC deaminase across treatments and timepoints. These strains represent the most likely biological culprits contributing to ACC removal in the treatment microcosms. Analysis of the broader gene expression of these two MAGs suggests the α-ketobutyric acid generated from ACC deaminase activity could be used to support isoleucine synthesis for these microorganisms (Fig. 6B; Additional File 3). Additionally, the Ornithinibacter MAG expressed genes for participation in arabinan utilization, a plant polysaccharide found in high concentrations in roots [61], as well as the genes for nitrate reduction (Nap). These paired metabolite and genome-resolved expression data uncover how cover crop ACC exudation can enrich microorganisms with the capacity to modulate soil carbon and nitrogen cycles that benefit plant health.

Microbial biosynthesis of indole derivatives, like auxins, can impact plant health directing proper plant growth and development [62, 63]. Given that IAA was only a component of the hairy vetch exudate pool, the IAA detected in sorghum, cereal rye, and rapeseed was likely microbially generated. This compound and other closely related indole intermediates can be synthesized via five microbial pathways derived from tryptophan (Trp), all converging on the final oxidation step converting indole acetaldehyde (IAAld) to IAA by aldehyde dehydrogenase (ALDH) [62, 64]. We profiled genes for IAA pathways in the ARM database and found 4 pathways were represented in our MAGs (Fig. 7A).

Of the 4 pathways, the indole-3-pyruvate (IPA) pathway was likely not contributing to the metabolite signal. We found 10 MAGs carrying IPA decarboxylase, yet we found no expression of either gene in the pathway (Tam1 and ipdC; Fig. 6A). For the second pathway, the indole-3-acetamide (IAM) pathway, we failed to identify any MAGs that encoded both genes in this pathway. However, it has been suggested that this could be community metabolism [65], thus we still examined the gene expression patterns in our data. One MAG from the unnamed genera belonging to the Geminicoccaceae family expressed iaaM and a different MAG from the Burkholderiales order expressed iaaH. Expression of these two genes was highest in the rapeseed, hairy vetch, and cereal rye exudate treatments at day 5 (Fig. 7B). Our combined metabolite and metatranscriptome data suggest that this pathway was likely contributing to the IAA-like metabolite signal in the cereal rye and rapeseed potentially by community metabolism.

The third pathway, the tryptamine (TAM) pathway, best explains IAA production in later time points of the sorghum treatment microcosms (Fig. 7B). This pathway is composed of two genes that were encoded across 116 MAGs, but 1 MAG from the Woeseiaceae family in the Pseudomonadota expressed both genes in the sorghum microcosms making this organism the most likely source for the increased IAA detected in sorghum metabolites at later time points. The fourth pathway, the indole-3-acetonitrile (IAN), is identified by the nitrilase gene. This gene was encoded by 55 MAGs, with expression detected in 5 MAGs from three genera in Actinomycetota, one in Pseudomonadota, and one in Chloroflexota. Expression was greatest in rapeseed at timepoints 5 and 21 (Fig. 7B). Holistically, our multi-omics indicates higher functional diversity expressed in the rapeseed treatment, where IAA synthesis genes were expressed by diverse organisms using two different pathways (IAN, IAM). Lastly, 66% of the MAGs in ARM had the potential for catalyzing the final step in IAA synthesis (ALDH), 15 of which expressed this gene only in exudate microcosms during (day 5) or after exudate stimulation (day 21; Fig. 7B). Our findings clearly validate a role for cover crops in producing a chemical environment that promotes microbial IAA production.

Gibberellic acids (GA) are important phytohormones that promote plant root and stem elongation [66]. Because the cereal rye metabolites showed a clear microbial GA₄ production signal, we next wanted to identify microbial taxa with this capacity. We found 176 MAGs that encoded at least one gene in the GA operon (Additional File 7 Fig. 15C), with 63 encoding CYP115, the gene responsible for converting GA₉ into bioactive GA₄ (Fig. 7C, D; Additional File 3). Supporting our GA₄ metabolite data, there were 5 MAGs that expressed CYP115 only in cereal rye metatranscriptomes at days 5 and 21 (Fig. 6A and E). A novel MAG in the Thermomicrobiales order (Chloroflexota) encoded 6 of 9 GA operon genes and actively expressed 5 (CYP112, CYP114, CYP115, GGPS, SDR) of the transcripts at day 21 in the cereal rye treatment and not in the control making this MAG the most likely culprit for GA₄ production (Fig. 7D). CYP115 expression by this MAG increased over time showing a significant positive correlation with GA₄ metabolite production (Fig. 7F; Additional File 7 Table 16). Our findings show the Thermomicrobiales order may be important in the production of GA₄, demonstrating the power of these laboratory experiments for uncovering new microbial physiologies with relevance to plant growth promotion.

Discussion

Cover crop root exudate stimulation expands the microbial genomic cataloging of the agricultural soil microbiome

Cultivation-independent, community-level genomic techniques are increasingly necessary to untangle the mechanistic relationship between soil microbial community metabolism and plant root exudates. However, obtaining genome-resolved metagenomic data from agricultural field soils is challenging due to microorganismal richness and high strain diversity, the latter especially confounds existing metagenomic recovery methods [67]. Additionally, the physical and chemical heterogeneity that changes across millimeter scales, also confounds the linkages of plant chemical inputs to microbial responses. Furthermore, the complexity of soil management systems (e.g., tillage and fertilization practices) that impact soil microbial function also obscures the effects of cover crop-microbiome interactions in field scale studies. Thus, using intermediate complexity laboratory-scale ecosystems (microcosms) allows for the temporal resolution of microbial dynamics and linkages of microbial gene expression patterns to changes in soil chemistry.

These microcosms were designed to illuminate the immediate impact of cover crop root exudation on agricultural soil microbiomes. As a result, this experimental design does not capture the myriad ways cover crops can shape soil microbiomes (e.g., through root architecture and plant litter decomposition), nor were these batch-operated experiments conducive for tracking the legacy effect (months, years) of cover crops exudates on soil microbial communities. A second limitation of this study was the constraints of using pure root exudates, which is tedious and low yield. Using a higher quantity of root exudate, although not as physiologically relevant considering root exudate mass per g of soil used here could have illuminated a greater metatranscriptome responsiveness of exudate-responsive taxa and contributed to a higher degree of microbial compositional changes in response to treatments as compared to time (Additional File 7 Fig. 6). We do, however, see that the effects of cover crop treatments on the slurry metabolome are indicative of metabolic adaptations in response to different root exudates which is possibly due in part of metabolomics as a higher sensitivity platform. While not a perfect replica of the complexity of cover cropping systems in a field environment, this research highlights how microcosms provided a controlled and replicable environment to construct a metabolic blueprint of agriculturally relevant taxa, many of which have not been genome sequenced or cultivated. Our approach and results provide mechanistic linkages between root exudates and soil microbial communities that can inform future research questions and the design of field-scale experiments.

Here, we establish ARM, the first microbial genome catalog of plant root exudate-responsive agricultural soil community members. Agricultural soils were used as an inoculum and stimulated over 6 days of various exudate additions with the metabolite and microbial community responses tracked for 15 days post-stimulation in two different experiments that collectively make up ARM (Fig. 1, https://doiorg.publicaciones.saludcastillayleon.es/10.5281/zenodo.11390808). Coupled to ARM is an extensive catalog of metabolomic and metatranscriptomic data from two cover crop exudate experiments, enabling us to decode the microbial membership, metabolic, and expression responses to exudate amendments from commonly used cover crops. Apart from providing multi-omic resources for ecologically relevant taxa, ARM MAGs offer access to genomic information from novel taxa as 70% of MAGs in the database lack a species assignment, and 87% are only recognized by an alphanumeric identifier (Fig. 2). Our intention was to establish a publicly accessible repository to support fellow researchers engaged in agricultural soil microbiome studies. This repository aims to serve as a practical resource for enhancing the retrieval of microbial genomes in agricultural fields, promoting emerging knowledge of relevant taxa, specifically those involved in microbial phytohormone biosynthesis.

This experiment highlighted new expressed capabilities of soil microbiomes. For example, in agricultural soils, Thaumarchaeota are well-known ammonia oxidizers that govern the fate of soil nitrogen. Here, we contribute to the growing genomic catalog of these ecologically relevant taxa [68, 69], adding 11 total MAGs that belong to 6 genera, 5 of which were only identified by alphanumeric identifiers. Our metatranscriptome data shows these microbes are active in exudate-treated and unstimulated soils, especially during the initial phases of the experiment, thus providing insights that can guide further cultivation efforts. Interestingly, however, the second phase of ammonia oxidation, nitrite oxidation, was clearly enriched only in exudate microcosms by non-comammox members of the Nitrospira_C genus, resulting in different nitrogen outcomes with and without exudation (Fig. 4). These findings highlighted how small chemical changes produced by cover crop exudation could result in differential ecosystem outputs, as others have shown at both the laboratory and field scales [70, 71].

Additionally, our phytohormone expression data revealed new genera that could be plant growth-promoting taxa. For example, our data showed that members of currently undefined genera Saprospiraceae (Bacteroidota), GWC2-71–9 (Gemmatimonadota), Geminicoccaceae, and Woeseiaceae (Pseudomonadota) expressed genes for IAA or GA₄ production (Fig. 6). Only two exudate-discriminant MAGs expressed the capacity for producing multiple phytohormones, one of these was the novel Woeseiaceae. Together our gene expression data show how cover crops through exudation can stimulate changes in nutrient cycling and phytohormone production.

Root exudate amendments alter microbial community function and metabolisms

Plant root exudates can influence the rhizosphere microbiome, an area of dynamic microbial activity around plant roots, where the abundance of nutrients stimulates diverse organisms to colonize the rhizosphere [15, 19, 28, 72, 73]. This hotspot of biological activity has the potential to influence plant and soil functional outcomes through the manipulation of root exudates. However, the intentional manipulation of the rhizosphere microbiome is limited in practice and has historically involved bacterial or fungal soil inoculation, such as inoculating seeds with plant growth-promoting microorganisms, to encourage proliferation without the need for external chemicals [74,75,76]. The primary pitfall to inoculant-focused strategies is the lack of inoculant resilience to competition or lack of proliferation due to unfavorable soil conditions which decreases their impact. Alternatively, successful experiments using the complementary approach of altering root exudate chemistry to induce changes in the microbiome are sparse but suggest great potential [77]. Thus, understanding how variability in root exudation across plant species and genotypes can impact plant-microbial associations will provide important insights for exudate-focused intentional microbiome manipulation to enhance primary crop performance and sustainability. Here, our combined metabolite and metatranscript data from cover crop exudate and non-exudate microcosms revealed important microbial roles in the degradation of ACC and the production of IAA and GA₄.

The phytohormone ACC is a precursor to the synthesis of ethylene in plants and has been shown to enhance stress resistance and relieve ethylene inhibition of plant growth. Our metabolite and metatranscript data suggested ACC consumption was possible in all exudate-treated micrcosms. ARM contained 13 MAGs from the phyla Actinomycetota (n = 9) and Pseudomonadota (n = 4) that encoded the acdS gene (Fig. 6A), deaminating ACC into α-ketobutyric acid and ammonia (Fig. 6B). Again, highlighting metabolic novelty uncovered in these reactors, for all 9 of the Actinomycetota we could not find prior accounts for ACC deamination. This has important ramifications for designing microbiome-informed cover cropping strategies, as there is a growing awareness that ACC not only acts in ethylene regulation, but also attracts ACC-degrading bacteria to establish a positive rhizosphere interaction [78,79,80,81]. In fact, this gene has been cited as an important trait for plant growth-promoting rhizobacteria in the Bacillus and Pseudomonadota phyla, though our results suggest members of the Actinomycetota phylum maybe play a similarly important role in ACC deamination [60].

The important endogenous plant auxin, IAA, is necessary for directing proper plant growth and development with production capacity being widespread among plant-associated microorganisms [62]. Within ARM, we detected genomic evidence for four of the five known Trp-dependent IAA biosynthetic pathways [63]. We found 215 MAGs encoded for at least one gene in any IAA pathway and 29 MAGs expressed one gene within any IAA pathway. Further, previous studies have demonstrated that bacteria are primed to produce IAA via the IPA or IAM pathways [63, 82], but our results highlight that the genomes analyzed here more dominantly encoded and expressed genes from the TAM and IAN pathways (Additional File 7 Fig. 14C–D). Out of the 15 MAGs expressing genes for IAA synthesis, 20% and 33% expressed DDC and MAO, respectively (TAM pathway), while only 7% expressed iaaM and iaaH (IAM pathway), and 33% expressed nitrilase (IAN pathway) (Additional File 7 Fig. 14C–D).

MAGs expressing DDC include 3 novel genera from diverse phyla (Gemmatimonadota, Bacteroidota, and Pseudomonadota) and one novel species of Actinomycetota. Similarly, two novel genera of Pseudomonadota, and novel species of Nitrospirota, Bacteroidota, and Actinomycetota expressed MAO. Of these MAGs, only members of Pseudomonadota have been found to produce IAA via the TAM pathway, leaving open the possibility that more diverse members may play a role in this pathway in agricultural soils [62, 83,84,85]. The two MAGs expressing iaaM and iaaH (IAM pathway) are both members of the Pseudomonadota phylum (Gammaproteobacteria and Alphaproteobacteria). Many species of Gammaproteobacteria, namely Pseudomonas species, are known to produce IAA through the IAM pathway suggesting the phylogenetically related Alphaproteobacterial members could have similar genomic potential [62, 86, 87]. Finally, our results suggest that MAGs encoding the IPA or IAM pathway were primarily from the Pseudomonadota, Acidobacteroita, or Bacillota phyla, as others have found in agricultural soils, while MAGs encoding the TAM or IAN were much more diverse (Additional File 7 Fig. 14C–D) [62, 84, 88, 89]. Our findings show that IAA production in soils is functionally redundant and that multiple cover crop exudation strategies can upregulate indole biosynthetic routes across diverse phyla which may serve to benefit plant growth [24, 89]. Future research to identify the exudate compounds that directly stimulate microbial phytohormone production could lead to new targets for cover crop breeding and their integration into biological precision agriculture approaches.

The GA phytohormones act endogenously in plants to promote seed germination, seedling growth, flowering, and leaf expansion and influence plant–microbe interactions [90, 91]. Plant-associated bacteria can produce GA to stimulate primary root elongation and lateral root extension [92,93,94]. Here, we found cereal rye exudate amendments increased the microbial production of bioactive GA₄ over the timeseries, with two other GA intermediates produced during the exudate-addition phase (Fig. 4E; Additional File 7 Fig. 15A–B). In terms of GA genetic potential, the core GA biosynthetic gene operon contains 8 genes and is widely distributed in rhizobia whereas CYP115 has only been found in certain species and is required for GA₄ synthesis (Additional File 7 Fig. 15C) [51, 95]. While 5 MAGs encoded and expressed CYP115 exclusive to the cereal rye treatment, the MAG most likely responsible for GA₄ production in cereal rye was a member from the Thermomicrobiales order (Chloroflexota). There are limited reports about the Thermomicrobiales in agricultural soils, yet 16S rRNA gene surveys have suggested members of this phylum are enriched in arid agricultural soils much like the soil inoculum used here [96]. Our findings uncovered bacteria that could be explored as plant growth-promoting targets for enhancing bioactive GA₄ production, potentially creating a favorable legacy environment for cash crop growth following a period of cereal rye cover cropping.

ARM is an important first step to begin illuminating the complex relationship between cover crop exudation and associated microbial membership and functional responses. We identified how different cover crop root exudates can alter the metabolic outputs in soils, and key pathways in which these metabolic changes have the potential to impact plant health. We have also resolved key microbial genomes and their transcriptome-based responses. A comparison of the genomes in ARM revealed similarity to prior laboratory and field-based work, highlighting that some of these key taxa are present in Colorado agricultural soils [71]. However, we recognize additional research is needed to ascertain the broader value of ARM to field relevant scenarios, spanning geographical regions and agroecosystem management practices. Future research can include examination of key responsive ARM genomes in cover crop trials via linkages of 16S rRNA sequences, or via q-RT-PCR responses for phytohormone genes uncovered here. This would provide additional evidence and/or validation of these lab scale responses at a field scale and also demonstrate how high-resolution techniques used here can generate hypotheses that scale to high throughput, targeted field-relevant responses.

Conclusions

Integrating a simplified and manageable laboratory soil system with advanced high-resolution multi-omics techniques afforded new insights into the intricate interplay between plant root exudates, the soil microbiome, and the capacity for soil microbial communities to contribute to phytohormone pools. We first provide ARM, a public agricultural genome database, creating a community resource for use in related microbiome-cover crop workflows that capture over 300 microbial genomes. Of these genomes, 215 encode hallmark genes of phytohormone metabolism, advancing the knowledge of the diversity of the physiology. Using ARM to contextualize our metabolite and metatranscriptome data, we then demonstrate cover crop exudate amendments applied at a physiologically relevant concentration impact microbial metabolism and membership function. Future studies can scale these hypotheses to the field to analyze the impact of plant genotype, root exudates, and associated microbial communities on overall plant performance. Additional research can detangle the precise chemical cues in these exudate cocktails that elicit direct microbial metabolic outcomes beneficial for soil and plant health. In summary, the results of this study represent an important step toward decoding the complex chemical language between plants and their rhizosphere microbial communities, a translation required to optimize these interactions for enhancing the environmental sustainability and productivity of agricultural systems.