Symbiosis with arbuscular mycorrhizal (AM) fungi improves plant nutrient capture from the soil, yet there is limited knowledge about the diversity, structure, functioning, and assembly processes of AM fungi-related microbial communities. Here, 16S rRNA gene sequencing and metagenomic sequencing were used to detect bacteria in the rhizosphere of Lotus japonicus inoculated with and without AM fungi, and the L. japonicus mutant ljcbx (defective in symbiosis) inoculated with AM fungi in southern grassland soil. Our results show that AM symbiosis significantly increased bacterial diversity and promoted deterministic processes of bacterial community construction, suggesting that mycorrhizal symbiosis resulted in the directional enrichment of bacterial communities and established a stable rhizosphere bacterial community. AM fungi promoted the enrichment of nine bacteria, including Ohtaekwangia, Niastella, Gemmatimonas, Devosia, Sphingomonas, Novosphingobium, Opitutus, Lysobacter, Brevundimonas, which are positively correlated with NPK-related parameters. Through a functional identification experiment, we found that six of these genera, including Brevundimonas, Lysobacter, Ohtaekwangia, Sphingomonas, Devosia, and Gemmatimonas, demonstrated the ability to mineralize organophosphate and dissolve inorganic phosphorus, nitrogen, and potassium. Our study revealed that AM fungi can regulate rhizosphere bacterial community assembly and attract specific rhizosphere bacteria to promote soil nutrient turnover in southern grasslands.

Greenhouse experiments

Experiment 1 was designed to explore whether AM fungi inoculation affected the bacterial community in the rhizosphere. The factorial experimental design included the addition of AM fungi and the non-addition of AM fungi (CK). Lotus japonicus (Miyakojima MG20) plants were selected and grown at 22°C to 24°C under a 16-h-day/8-h-night cycle. We collected soil samples from a southern grassland in Zhaotong, Yunnan Province, Southwest China (103°63′E, 27°66′N; Figure S1). The dominant native AM fungi in this soil were Acaulospora, as identified using special primers (Lee et al., 2008; Sato et al., 2005; Zhang et al., 2021) (Table S1). Surface soil (0–20 cm) was collected using a “five points” sampling strategy in a 25 m × 25 m field area. Soils were immediately transported for storage at 4°C. Soils were dried sufficiently to pass through a 2-mm sieve and plant tissues were removed. Soil samples were divided into two parts: one portion was used for chemical analysis, and another portion was used for the greenhouse experiment. Each pot contained 500 g of non-sterile soil. Three AM fungi, including Glomus mosseae, G. intraradices (alternatively Rhizophagus intraradices), and G. versiforme, were used in the experiment.

These AM fungi are widely used in scientific research and practical applications and have different symbiotic abilities (G. versiforme>G. intraradices>Glomus mosseae) with Lotus japonicus. AM fungi spores were propagated on maize (Zea mays L.), L. japonicus, and white clover (Trifolium repens L.) in a sand-ceramsite mixture (3:1w/w) for 4 months before the experiment. The matrix containing AM fungal propagation passed through sieves of 0.8mm, 0.25mm, and 0.055mm in turn. The inoculum in the bottom sieve including AM fungi spores and sand was used as AM fungi inoculum. Before use, we counted the spore density in the inoculum. Each gram of Glomus mosseae, G. intraradices, and G. versiforme inoculum contained approximately 60, 80, and 80 spores, respectively. A total of 250 AM fungal spores were added per pot in AM fungi addition experiments. AM fungi inoculum was mixed with soil before planting plants. Eight replicates for each treatment were set.

Experiment 2 was designed to determine: (1) if mycorrhizal symbiosis can promote the absorption and utilization of nitrogen, phosphorus, and potassium by L. japonicus; (2) if mycorrhizal symbiosis shapes microorganisms in the rhizosphere. Plants L. japonicus (wild-type, Gifu-129) and mycorrhiza-defective mutant ljcbx plants (Gifu-129 background) were selected to grow in non-sterile soil inoculated with G. versiforme (250 spores per pot), which have the strongest symbiotic ability with L. japonicus compared with Glomus mosseae and G. intraradices. 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 growth conditions were the same as in Experiment 1. There were eight replicates for WT and mutant ljcbx treatments in Experiment 2.

Plant and soil analyses

In each treatment, plants were randomly collected at the final harvest and then divided into roots and shoots to analyze for nutrient uptake. Plant organs were dried at 105°C for 30 min and then oven-dried at 75°C for 48 h to a constant weight (dry weight, DW). Samples were then ground to pass through a 0.25-mm sieve for chemical analysis. For plant N, P, and K measurements, samples were digested firstly in hot sulfuric acid with hydrogen peroxide as an additive (H2SO4–H2O2 method). In short, accurately weighed 0.1g ground sample in the digestion tube, added distilled water to fully wet it, then added 5 ml sulfuric acid to soak it overnight, and then digested it at 375℃ for 3h, cooled the digestion tube, added a few drops of hydrogen peroxide until it is completely digested. The volume of the digested solution to 50 ml. Took 10 ml supernatant to a volumetric flask, added water to 20 ml, added two drops of 0.3% 2,4-dinitrophenol, and then sodium hydroxide solution was added to adjust the solution to yellow (pH 6-7). Total N, P, and K were measured using a CleverChem 380 discontinuous water quality analyzer (DeChem-Tech, Hamburg, Germany). For soil nutrient analysis, the soil was first passed through a 2-mm sieve. Soil-available phosphorus (AP) was extracted with sodium bicarbonate and available potassium (AK) was extracted with ammonium acetate solution. Soil ammonium nitrogen (NH4+-N) and nitrate nitrogen (NO3-N) were extracted with 1 mol·l−1 KCl. Soil NH4+-N, NO3-N, AK, and AP were measured using a CleverChem 380 discontinuous water quality analyzer (DeChem-Tech, Hamburg, Germany).

Rhizosphere soil separation and 16S rRNA gene sequencing

Rhizosphere soil was collected according to established protocols (McPherson et al., 2018) after plants had grown in soil for 5 weeks. In brief, soil loosely attached to roots was removed with gentle shaking. Steps to collect rhizosphere soil were performed on ice. Roots with adhered soil (rhizosphere 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. The total DNA of rhizosphere samples was extracted using a PowerSoil® DNA isolation kit according to the manufacturer’s instructions. Miseq sequencing targeting the bacterial 16S rRNA gene V3-V4 region was performed on an Illumina MiSeq platform by the company Sangon Biotech (Shanghai, China).

The 16S rRNA gene V3-V4 region was amplified using the16S rRNA gene primers 341F (5′-CCTACGGGNGGCWGCAG-3′) and 805R (5′-GACTACHVGGGTATCTAATCC-3′). The amplification procedure is as follows: 94 °C for 3 min, followed by 5 cycles consisting of 94 °C for 30 s, 45 °C for 20 s, 65℃ for 30 s, then 20 cycles consisting of 94℃ for 20 s, 55℃ for 20 s, 72℃ for 30 s, finally, 72 °C for 5 min. There was an average of 40000 reads per sample via 16S rRNA sequencing. Amplicon sequencing reads were quality-filtered and demultiplexed using QIIME (Caporaso et al., 2010) and USEARCH (Edgar, 2010) pipelines. Demultiplexed files from each library were processed separately with the “DADA2” (Callahan et al., 2016) pipeline in R (version 3.6.0). Amplicon sequence variants (ASVs) in each sample were inferred followed by trimming and filtering. Trim parameters were as follows: trim left 0 bases and truncate at base 250 for both forward and reverse reads. ASVs were then aligned to the SILVA database (Quast et al., 2013) for the taxonomic assignment. Sequences were clustered into OTUs at 97% identity using UPARSE algorithm.

To detect genes related to metabolism by amplicon sequencing, the Kyoto Encyclopedia of Genes and Genomes (KEGG, http://www.kegg.jp) database and Pfam (Finn et al., 2013) database were used. DIAMOND (Buchfink et al., 2015) was used to compare the protein sequences of the gene set with KEGG database to obtain the KEGG orthology (KO), and the annotation of pathway and module of the sequence was obtained subsequently. The abundance of each functional level of KEGG in each sample was counted. Using KEGG and Pfam, functional genes were searched by protein domain. The functional gene set involved in N and P cycle and related KO information was constructed. Finally, according to the KO result of functional prediction of the amplicon through PICRUSt (version 1.1.4) using the KEGG database, the N and P functional genes were extracted and the heat map was drawn.

Metagenome sequencing

Metagenomic shotgun sequencing was performed on an Illumina HiSeq 2500 platform at Sangon Biotech (Shanghai, China). A total of six samples were conducted metagenomic sequencing, including samples not inoculated with AM fungi (CK, three replicates) and samples inoculated with G. versiforme (three replicates). Finally, an average of 6 GB sequencing data were obtained from each sample. Raw sequencing data were evaluated with FastQC (version 0.11.8) software (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) and filtered using Trimmomatic (version 0.36) (Bolger et al., 2014) to obtain relatively accurate and effective data. IDBA_UD was used to splice clean reads of each sample to assemble long sequence contigs (Peng et al., 2012). Contigs were obtained according to overlapping relations between reads (minimum overlap=100). Assembly results of multiple K-mers were comprehensively evaluated, and the best K-mer assembly results (k-mer=31) were selected. Prodigal (version 2.60) (Liu et al., 2013) was used to predict ORFs (Open Reading Frames) of splicing results, and sequences with length ≥ 100 bp were selected and translated into amino acid sequences. We used CD-HIT software for de-redundancy to obtain a nonredundant gene set for the gene prediction results of each sample (Li & Godzik, 2006). 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). All reads sequences were mapped back to assembled contigs to obtain contig coverage, and then the contigs were clustered by integrating the information of GC content, nucleic acid composition, single copy gene, and other genomic characteristics (Zeng et al., 2019) to separate contig sequences belonging to different genomes. Finally, high-quality Binners (MAG, metagenome assembled genomes) were obtained. To further improve the quality and integrity of bins, we re-mapped the original sequence based on the obtained bins and optimized the assembly. Using SPAdes for fine assemble to remove contigs less than 1500 bp, and then evaluate the quality by checkM (version 1.1.3) (Parks et al., 2015). Each bin was annotated against the National Center for Biotechnology Information (NCBI) nonredundant Nucleotide Sequence Database (NT) using BLASTn (You et al., 2022).

Species annotation was conducted by searching against NCBI non-redundant protein sequences (NR, http://ncbi.nlm.nih.gov/) database via DIAMOND blastp (e-value<1e-5) (Buchfink et al., 2015). The translated proteins from all detected coding regions of each metagenome were annotated by searching against three custom databases. Kyoto Encyclopedia of Genes and Genomes (KEGG, http://www.kegg.jp) database was searched using GhostKOALA (e-value<1e-5) (Kanehisa et al., 2016) to predict functional pathways and KEGG Orthology (KO) was obtained. Functional assignments were also conducted through the SEED database (http://www.theseed.org/wiki/Main_Page) using DIAMOND (e-value<1e-5). The evolutionary genealogy of genes: Non-supervised Orthologous Groups (eggNOG, http://eggnogdb.embl.de/) database was searched via DIAMOND (e-value<1e-5), and the Clusters of Orthologous Groups of proteins (COG) corresponding to the genes were obtained. On this basis, the gene sets were annotated and classified, and the abundance of each functional level of COG in each sample was counted. Then, using perl to retrieve COG annotation information corresponding to the core species. Gene abundances (gene counts per KO) were normalized to the number of amino acids detected in the metagenome of each sample.

Isolation and identification of rhizosphere bacteria

To isolate rhizosphere bacteria, we resuspended 1 g of rhizosphere soil in sterile phosphate buffer and incubated the mixture for 30 min with shaking at 150 rpm. Strains harvested from rhizospheres were cultured on three media: nutrient agar (NA), Reasoner’s 2A agar (R2A), and trypticase soy agar (TSA). A 1-ml suspension was taken from each mixture for serial dilution on NA, R2A, and TSA. One representative colony (determined by its morphology) from each plate was selected, purified, and stored at −80°C in R2A medium containing 20% glycerol. The complete 16S rDNA was amplified via 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 by 16S rRNA gene sequencing, which was performed by Sangon Biotech. We analyzed resulting high-quality sequences with BLASTn (NCBI; http://blast.ncbi.nlm.nih.gov) to confirm the identity of each bacterium. Gram staining was performed with a three-step Gram stain procedure kit (Solarbio Co., Ltd., Beijing, China), in strict accordance with the standard procedure. Stained cells were observed with a compound microscope.

Determination of rhizosphere bacterial functions

Fresh soil samples (10 g) were suspended in 90 ml of sterile saline solution (NaCl, 0.9% w/v), and mixtures were shaken (150 rpm) for 30 min at room temperature. Ten-fold dilutions were made in sterile saline solution, and 0.1 ml was spread on the appropriate medium. Bacteria were cultured in standard NA (Cultimed, Spain) for 48 h at 30°C. Cultivable bacteria and organophosphorus-degrading bacteria (OPDB) were estimated using a standard dilution-plating procedure. The OPDB were cultivated in modified Pikovskaya’s agar medium (pH 7.0–7.5) containing the following (l⁻¹): glucose, 10.0 g; (NH4)2SO4, 0.5 g; NaCl, 0.3 g; KCl, 0.3 g; yeast extract, 0.5 g; MgSO4·7H2O, 0.3 g; FeSO4·7H2O, 0.03 g; MnSO4·4H2O, 0.03 g; CaCO3, 5.0 g; soybean lecithin (P-lecithin), 15.0 g; and Bacto-agar, 15.0 g. The phosphate-solubilizing ability of each isolate was tested using insoluble tricalcium phosphate (Ca3(PO4)2) as the sole P source in Pikovskaya’s medium. Phosphate solubilization was assessed by spot-inoculating individual bacterial isolates on Pikovskaya’s medium containing the following (l⁻¹): yeast extract, 0.50 g; dextrose, 10.0 g; calcium phosphate, 5.00 g; ammonium sulfate, 0.50 g; potassium chloride, 0.20 g; magnesium sulfate, 0.10 g; manganese sulfate, 0.0001 g; ferrous sulfate, 0.0001 g; and Bacto-agar, 15.00 g. Media pH was adjusted to 7.0 before autoclaving. Mixtures were poured into sterile petri plates (25 ml of agar per plate). Plates were incubated for 4 to 5 d at 28 ± 1°C. A clear zone around a bacterial colony was considered a positive indication of phosphate solubilization (Pande et al., 2017). To test the bacterial ability to solubilize K, all colonies growing on plates were screened on modified Aleksandrov agar medium containing the following (l−1): glucose, 3.50 g; MgSO4.7H2O, 0.50 g; CaCO3, 0.10 g; FeCl3, 0.0005 g; Ca3PO4, 2.00 g; insoluble mica powder, 1.0 g, as the K source; and agar, 15.00 g. Colonies selected from plates were spot-inoculated on new plates and incubated at 28°C for 7 d. Nitrogen-fixing bacteria were isolated using solid N-free medium (NFM) containing the following (l⁻¹): K2HPO4, 0.20; KH2PO4, 0.50 g; MgSO4.7H2O, 0.20 g; FeSO4·7H2O, 0.10 g; Na2MoO4·2H2O, 0.005 g; NaCl, 0.20 g; glucose, 10.00 g; and noble agar, 15.00 g. Serial dilutions were spread on freshly prepared solid NFM in triplicate. Plates were incubated for 5 to 7 d at 28°C. Colonies were selected based on morphology and were streaked three times on NFM to ensure purity before they were isolated.

Data analysis

All analyses and data visualization were performed in R version 3.6.0 (Hu et al., 2022; Team, 2019). Differences in the alpha diversity richness were estimated using one-way and two-way ANOVAs. Differences in the taxonomic composition of bacterial communities were visualized using non-metric multidimensional (NMDS) analyses. Alpha diversity and NMDS were calculated using the ‘vegan’ package (version 2.0-10) (Oksanen et al., 2013) in R. To describe potential co-occurrence between all OTUs in the rhizosphere, co-occurrence networks of each sample were constructed using the ‘cooccur’ package in R. Using psych package in R, the genus level abundance matrix was transformed into species correlation matrix. Then, the Gephi (https://gephi.org/) was used to calculate the information of nodes and edges of the network. At last, network characteristic data was exported to draw figures. For Sloan neutral model (Rezki et al., 2017; Sloan et al., 2006) 95% binomial confidence intervals for the neutral model based on the Wilson method using the ‘Hmisc’ package in R were constructed (Morris et al., 2013). The deterministic processes (selection and environmental filtering) and stochastic processes (dispersal and drift) were assessed by comparing the βNTI values: deterministic processes (|βNTI|≥2), stochastic processes (−2 < βNTI < +2) (Dini-Andreote et al., 2015). Redundancy analysis was performed using ‘vegan’ package (version 2.0-10). In addition, linear discriminant analysis effect size (LEfSe) analysis was performed using LEfSe software to identify specific bacterial genera that are associated with AM fungi.

Data availability

Raw transcriptome data have been deposited in the National Center for Biotechnology Information (NCBI) SRA database with the BioProject accession number PRJNA791779 (SRR17335095-SRR17335100). The 16S rRNA gene sequencing raw data were deposited in the NCBI SRA database under accession number PRJNA792500 (SRR17344331-SRR17335100). Metagenomic sequences were deposited in the SRA database, which is available at NCBI with accession number PRJNA792495 (SRR17333937-SRR17333942).

Code availability

Scripts used for computational analyses described in this study are available at https://github.com/xuyunjian1992/R-code-and-software, to ensure replicability and reproducibility of these results. No unpublished algorithms or methods were used.