Arbuscular mycorrhizal (AM) fungi traditionally form symbioses with most plant species. Although AM fungi have critical effects on microbial communities, the pathways showing how AM fungi shape rhizosphere bacterial communities and their functions are rarely explored. Through three systematic experiments, AM fungi-bacteria interactions were first investigated in the rhizosphere of Lotus japonicus, then the interactions were confirmed by a second experiment with wild-type and a mycorrhiza-defective mutant ljcbx of L. japonicus. The mechanisms were presented by adding core bacteria and AM fungi to the plant rhizosphere with the third experiment. We found that AM fungi-bacteria interactions enhanced host plant growth and identified a core bacterial group that uniquely enhanced host plant growth. Adding core bacteria and AM fungi promoted host growth and nutrient acquisition compared to adding AM fungi or core bacteria independently. Allelopathic substances secreted by AM fungal colonizing host roots to recruit the rhizosphere bacteria were detected by the multi-omics joint analysis, showing that arachidonic acid was the main allelopathic substance that affected AM fungi–bacteria interactions. Our findings provide direct evidence that mycorrhizal infection simulated root exudation, such as arachidonic acid, recruited a beneficial microbiome to the host rhizosphere, increasing plant growth and soil nutrient turnover.

Greenhouse experiments

Experiment 1 was designed to determine whether rhizosphere bacteria affected the plant growth-promoting effects of AM fungi in the greenhouse. The factorial experimental design included the addition of AM fungi and one soil with two treatments (sterilized and unsterilized). To sterilize the soil, the soil was autoclaved at 121°C for 1 h and the procedure was repeated after 2 days. Wild-type (WT) L. japonicus plants were selected. In the dark, surface-sterilized seeds of L. japonicus were grown on 0.8% agar for 5 days (d) and were then transferred into pots containing AM fungal inoculum.  Plants were grown at 22°C to 24°C under 16 hours of light: 8 hours of dark photoperiodic lighting for 40 d. The soil samples for the first experiment were collected from a southern grassland in Yunnan Province, China (103°63′E, 27°66′N). Each sample was extracted from a depth of 0–20 cm. The mean soil organic carbon content (SOC) was 60.8 g kg^-1; the mean soil total nitrogen content (TN) was 4.7 g kg^-1, and the mean soil pH value was 4.7. AM fungi included Glomus mosseae, G. intraradices (alternatively R. intraradices), or G. versiforme. Each pot contained 500 g of either sterile or non-sterile soil and 250 spores of fungal inoculum. The experimental design as Fig. 1a shows. Experiment 1 had a total of 64 pots with 8 replicates for each treatment.

Experiment 2 was designed to determine if mycorrhizal symbiosis shapes plant growth and enrichment of N and P cycling microorganisms in the rhizosphere. Wild-type L. japonicus and mutant ljcbx plants were grown in nonsterilized soil inoculated with G. versiforme (250 spores per pot). The growth conditions were the same to Experiment 1. Mycorrhiza-defective L. japonicus mutant ljcbx (30057597) plants were kindly provided by Dr. Stig U. Andersen and Jens Stougaard (Lotus Base, https://lotus.au.dk/) (Małolepszy et al. 2016). The ljcbx is a mutant of LjCBX1, which has a specific function in the symbiotic exchange of nutrients in mycorrhizal L. japonicus (Xue et al. 2018). Trypan blue staining used for the detection of mycorrhizal symbiosis was performed according to the methodology previously described by Liu et al. 2020. Mycorrhizal colonization (%) was quantified using the gridline intersect method previously described (Wang et al. 2014). Experiment 2 had a total of 16 pots with 8 replicates of WT and mutant ljcbx, respectively.

Experiment 3 included two steps to find a simplified synthetic community (SynCom) and determine if the combined effects of AM fungi and SynCom can promote plant growth. To isolate rhizospheric bacteria, rhizosphere soil suspensions were serially diluted on nutrient agar (NA), Reasoner’s 2A agar (R2A), and trypticase soy agar (TSA). The complete 16S rDNA was amplified using PCR with universal bacterial primers 27F (5′-AGA GTT TGA TCC TGG CTC AG-3′) and 1492R (5′-ACG GCT ACC TTG TTA CGA CTT-3′) (Defez et al. 2017). Genus-level classifications of isolates were confirmed through 16S rRNA amplicon sequencing performed by Sangon Biotech (Shanghai, China). High-quality sequences were analyzed using the NCBI BLASTn (NCBI; http://blast.ncbi.nlm.nih.gov) database to confirm the identity of each bacterium. Subsequently, a SynCom consisting of seven core genera of bacteria was constructed. The standards for core bacterial genera used to construct the SynCom included (1) the relative abundances of bacteria that enriched in rhizospheres of wild-type L. japonicus colonized by three AM fungi, but not in rhizospheres of AM-colonized ljcbx1; (2) enriched bacteria in the rhizosphere of mycorrhizal plants containing functions associated with N, P, and K metabolisms. To explore whether the combined effects of AM fungi and SynCom promote plant growth, the SynCom was mixed with equal volumes of each bacterium (~10⁸ CFU·ml^-1). The SynCom inoculum was adjusted to an optical density (OD₆₀₀) of 0.02, and each pot was inoculated with 200 µl of bacterial solution or 10 mM MgCl2 (Pfeilmeier et al. 2021). The plants were planted in pots with sterilized soil. Experiment 3 had 64 pots with 8 replicates for each treatment.

Plant and rhizosphere soil analyses

Plant shoots and roots were collected in each pot and oven-dried to constant mass at 70°C. Rhizosphere soil was determined to be the 1 mm width root-adhering soil and was collected using the method previously described (McPherson et al. 2018). In brief, plants were removed from each pot, and soil loosely attached on roots was removed with gentle shaking. The remaining root-adhering soil was the rhizosphere, which was approximately 1 mm thick. Roots with adhering soil were placed in a sterile 50-ml falcon tube containing 30 ml of sterile, pre-cooled PBS (phosphate-buffered saline) buffer (pH 7.5). The tube was mixed in a vortex apparatus for 15 s, and the turbid mixture was passed through a 100-μm aseptic nylon mesh strainer (to remove root fragments and large soil particles) into a new 50-ml tube, which was centrifuged for 5 min at 12,000 ×g at 4°C. The supernatant was discarded, and the collected soil, defined as rhizosphere soil, was frozen with liquid nitrogen and stored at -80°C.

Available nutrients and total NPK were measured using a CleverChem 380 discontinuous water quality analyzer (DeChem-Tech, Hamburg, Germany). Soil enzymatic activities (urease, ammonia monooxygenase, nitrogenase, pyrroloquinoline-quinone synthase, and phytase) were determined using an enzyme-linked immunosorbent assay (ELISA) kit (Mlbio, China). Total genomic soil DNA (gDNA) was extracted using a FastDNA SPIN kit (MP Biomedicals, France) following the manufacturer’s instructions. Ammonia monooxygenase (AMO, which catalyzes the oxidation of ammonia to hydroxylamine), nitrogenase (nifH, which catalyzes nitrogen fixation), urease (ureC, which catalyzes urea to ammonia and carbon dioxide), pyrroloquinoline-quinone synthase (pqqC, which catalyzes the final step of the pyrroloquinoline quinone biosynthesis), and phytase (BPP, associated with the hydrolysis of phytate) were functional genes assessed through quantitative PCR (qPCR) due to their involvement in nitrification, N fixation, organic N mineralization, organic P mineralization (phytic acid mineralization), and inorganic P dissolution. Primer sequences for target genes were designed according to previous studies (Table S1) (Cotta et al. 2016; Zheng et al. 2017). Quantification was performed with TB Green Premix Ex Taq II (Tli RNaseH Plus) on ABI QuantStudio 7 (Applied Biosystems, Germany). The conditions for qPCR were summarized in Table S1.

For 16S rRNA sequencing, sample DNA was extracted, and dual bar-coded MiSeq libraries of 16S rRNA V3-V4 gene amplicons were prepared and sequenced on an Illumina MiSeq platform. Demultiplexed files from each library were processed separately with the “DADA2” (Callahan et al. 2016) to construct a sequence count table (Table S2, S3). The count table, phylogenetic tree, taxonomy table, and sample metadata were combined in a “phyloseq” object for further analysis (McMurdie & Holmes 2013). Amplicon sequence variants (ASVs) were normalized to account for variable sequencing depth between samples using DESeq2 (Pfeilmeier et al. 2021). Relative abundances of each strain or tax on were calculated by proportional normalization of each sample by its sequencing depth.  A random forest model was conducted by using the "randomForest" R package to identify the biomarker taxa correlated with AM fungi symbiosis. Cross-validation was conducted using the rfcv() function to choose appropriate characteristics, and the varImpPlot function was used to determine the importance of characteristics in the classification (Li et al. 2022). The importance of characteristics and the cross-validation curve were visualized using the ggplot2 package (Wickham 2016). Both model training and testing used parameters as follows: Model, Random Forest; Error type, out-of-bag; Estimated error, 0.10000; Baseline error (for random guessing), 0.50000; Ratio baseline error to observed error, 5.00000; Number of trees, 500.

For metagenomic sequencing and downstream analyses, metagenomic shotgun sequencing of soil samples was performed on an Illumina HiSeq 2500 platform by Sangon Biotech (Shanghai, China). High-quality, raw sequencing data were assembled using IDBA_UD (Peng et al. 2012). CD-HIT was used for de-redundancy to obtain a nonredundant gene set (Li & Godzik 2006), and Prodigal was used to translate into protein sequences (Liu et al. 2013). Bowtie2 was used to align clean reads of each sample to nonredundant gene set sequences (Langmead & Salzberg 2012), and SAMtools was used to obtain aligned reads (Li et al. 2009). Functional analysis of metagenomic data was conducted using the SEED protein database (http://www.theseed.org/wiki/Main_Page) and using DIAMOND (Buchfink et al. 2015) with e-value cutoff of 1e^-5. The corresponding metagenomic data of these bacteria were extracted using perl software. The obtained species and functional gene annotations were shown in Table S4.

The metabolites were extracted as follows(Yang et al. 2017): 100 μL of each sample was placed into a 1.5 mL centrifuge tube and 300 μL of methanol containing 20 μL of the internal standard was added, vortex for 30 s; then sonicated for 10 min in the ice water bath, incubated at -20℃ for 1 h, and centrifuged 15min at 13000 rpm; transferred 200 μL of supernatant to a 2 mL injection bottle, mixed 20 μL of each sample into QC sample, and then took 200 μL for testing on the computer.  The metabolic substances were analyzed using UHPLC/MS system consisting of an Agilent 1290 Infinity II Series UHPLC (Agilent Technologies, Santa Clara, CA, USA), an automated multisampler module, a high-speed binary pump, and an AB Sciex TripleTOF® 5600 Mass Spectrometer (AB Sciex, CA, USA). Briefly, soil and root metabolic samples were injected onto an ACQUITY UPLC BEH Amide HPLC column (1.7 μm, 100 × 2.1 mm, Waters) with a flow rate of 0.5 ml·min^-1 at 55°C. In this mode, acquisition software (Analyst TF 1.7, AB Sciex) continuously evaluated the full scan of MS data as it was collected and triggered the acquisition of MS/MS spectra depending on preselected criteria (Table S5, S6). For RNA-sequencing, the total RNA of roots was isolated using Trizol reagent (Takara). Libraries were prepared using NEBNext UltraTM RNA Library Prep Kit for Illumina (NEB, USA) and were sequenced on the Illumina Nova. Clean reads were mapped to the L. japonicus Miyakojima MG20 reference genome (https://lotus.au.dk/data/download) using bowtie2 (Du et al. 2021). Gene expression levels were calculated in FPKM (fragment per kilobase of gene per million reads mapped). Differentially expressed genes were identified using the P package in DEGseq (Table S7). For each sample, three biological replicates were used. Data sets were analyzed in RStudio1.3.959 with R 3.6.0 (https://r-project.org/).