Effect of strigolactones on recruitment of the rice root-associated microbiome

Abstract

Strigolactones are endogenous plant hormones regulating plant development and are exuded into the rhizosphere when plants experience nutrient deficiency. There, they promote the mutualistic association of plants with arbuscular mycorrhizal fungi that help the plant with the uptake of nutrients from the soil. This shows that plants actively establish—through the exudation of strigolactones—mutualistic interactions with microbes to overcome inadequate nutrition. The signaling function of strigolactones could possibly extend to other microbial partners, but the effect of strigolactones on the global root and rhizosphere microbiome remains poorly understood. Therefore, we analyzed the bacterial and fungal microbial communities of 16 rice genotypes differing in their root strigolactone exudation. Using multivariate analyses, distinctive differences in the microbiome composition were uncovered depending on strigolactone exudation. Moreover, the results of regression modeling showed that structural differences in the exuded strigolactones affected different sets of microbes. In particular, orobanchol was linked to the relative abundance of Burkholderia–Caballeronia–Paraburkholderia and Acidobacteria that potentially solubilize phosphate, while 4-deoxyorobanchol was associated with the genera Dyella and Umbelopsis. With this research, we provide new insight into the role of strigolactones in the interplay between plants and microbes in the rhizosphere.

rice, strigolactones, rhizosphere, bacterial communities, fungal communities

Introduction

Plants and microbes are key players in belowground ecosystems and their existence is closely intertwined. Moreover, microbial communities greatly affect plant growth and health. Certain microbial species can solubilize nutrients making them accessible to plants, produce phytohormones that promote plant growth and/or increase resistance against biotic and abiotic stresses (Bruto et al. 2014, Trivedi et al. 2020). Plants, in turn, provide an energy source for the microbes through the exudation of a blend of hundreds, if not thousands, of compounds into the soil surrounding the roots, called the rhizosphere (Sasse et al. 2018). These compounds not only are used by soil microbes for their growth and reproduction but are also a means of chemical communication between plants and microbes (van Dam and Bouwmeester 2016). Several authors have reported that plants actively exude signaling molecules in order to recruit beneficial microbes that can help plants to overcome hostile environmental conditions (Planchamp et al. 2015, Abedini et al. 2021, Rizaludin et al. 2021).

A prime example of plant signaling molecules in the interaction with microorganisms is the strigolactones (SLs). SLs are a carotenoid-derived plant hormone that regulates plant development, including traits such as shoot branching and root architecture, in interplay with other hormones, such as auxin, abscisic acid and cytokinin (Waters et al. 2017). SL biosynthesis begins with conversion of all-trans-β-carotene to 9-cis-β-carotene by carotenoid isomerase DWARF27 (D27; Alder et al. 2012). Carotenoid cleavage dioxygenase 7 (CCD7) (rice D17, Arabidopsis thaliana AtMAX3) and CCD8 (D10, AtMAX4) catalyze intermediate conversions together resulting in the formation of carlactone, the last common precursor for all SLs. The formation of the large variety of SLs, identified to date, is subsequently mediated by MORE AXILLARY GROWTH 1 (MAX1) homologs, the oxidase LBO and a number of other cytochromes P450 (Wakabayashi et al. 2019, Bouwmeester et al. 2021).

Besides their role as endogenous plant hormone, SLs also act as a rhizosphere signaling molecule. SLs are secreted—especially under phosphorus (P) and nitrogen (N) deficiency—from roots of many plant species, such as tomato, rice, sorghum, alfalfa, lettuce and wheat (Yoneyama et al. 2007, López-Ráez et al. 2008, Jamil et al. 2012, Yoneyama et al. 2012). These SLs in the rhizosphere mediate the interaction between plants and symbiotic arbuscular mycorrhizal fungi (AMF) that colonize up to 90% of land plants and help them to uptake P from soil, by inducing hyphal branching, a process that precedes, and is needed for, root colonization (Akiyama et al. 2005, Begum et al. 2019). So far, over 35 different SLs have been identified in the root exudates of a range of plant species that usually produce blends of up to 7 different SLs by each plant species (Bouwmeester et al. 2021). Why plants exude so many different SLs is, however, unknown (Wang and Bouwmeester 2018). Akiyama et al. 2010 applied 44 different, synthetic and natural, SLs to Gigaspora margarita and showed that there is substantial variation in the hyphal branching activity induced by these SLs. It is unclear whether these activity differences play a role under field conditions.

Recently, a number of authors provided evidence that SLs may have an impact on the composition of the rhizosphere microbial community, using knockout mutants and overexpression lines for SL biosynthesis and/or perception genes in Arabidopsis, rice and soybean (Carvalhais et al. 2019, Nasir et al. 2019, Liu et al. 2020). Here, we want to further prove the hypothesis that SLs affect the plant microbiome and investigate whether structurally different SLs have different effects. Therefore, we analyzed the SLs of 16 rice genotypes, 11 varieties and 6 mutant lines, cultivated in two different soils, and assessed the bacterial and fungal communities recruited to the rhizosphere and root using 16S rRNA gene and Internal Transcribed Spacer (ITS) region amplicon sequencing, respectively. A significant effect of individual SLs on the composition of both bacterial and fungal communities was observed. Moreover, using constrained multivariate analysis and negative binomial regression modeling, correlations were pinpointed between the level of specific SLs and individual microbial taxa with putative functions that could potentially mitigate abiotic stress.

Materials and methods

Soil collection and their properties

To compare nutritionally rich and poor soil, we used two soil types—agricultural soil and forest soil—in this study. The agricultural soil was provided by the Netherlands Institute for Ecology (NIOO-KNAW), Wageningen, The Netherlands. Its chemical properties were described in detail in a previous study (Schlemper et al. 2017). The second soil was a nutrient-poor, acidic, forest soil, consisting of the top 20 cm layer of a natural heather field collected at Hoge Veluwe National Park, The Netherlands (52°10′24″N and 5°47′52″E) and sieved through a 2-mm mesh. The chemical and physical properties of the forest soil were analyzed at Eurofins Agro, Wageningen, The Netherlands (Table S1, Supporting Information). One day before transplanting seedlings, a 2.5-L size pot (17 cm Ø × 16 cm H) was filled with 2.4 kg of soil with a moisture content of 80% of its maximum soil water holding capacity.

Plant growth

All rice genotypes used in this study are listed in Table S2 (Supporting Information). Ten natural rice varieties differing in the amount of SLs produced under phosphate starvation were selected based on an earlier study (Jamil et al. 2012). Seeds of these varieties were provided by the International Rice Research Institute, Philippines. We also included mutant lines impaired in SL biosynthesis (d10, d17 and d27) and perception (d3 and d14) and its wild type (WT; Shiokari; Haider et al. 2018). Rice seeds were surface-sterilized in half strength of commercial bleach (approximate final concentration: 2%) for 30 min and washed five times with sterile demineralized water. The seeds were then germinated on half-strength Murashige and Skoog medium plates with 1.5% sugar and 0.9% agar in dark at 24°C for 4 days. Five seedlings of each genotype were finally transplanted per pot of each soil type and grown under controlled greenhouse conditions for 28 days at 25–28°C and 12/12 h day/night together with a pot without plant (bulk soil). Each treatment had five replicates, and all pots were watered with demineralized water every day to keep 80% water holding capacity until harvest (Fig. S1, Supporting Information).

Plant harvest

Two compartments, roots and rhizosphere soil of 32-day-old rice plants (vegetative stage), were harvested as previously described with small modifications (Lundberg et al. 2012). Briefly, plant roots were first cut from shoot and gently tapped until 1–2 mm soil layer was left on the root surface. Then, the root system was shaken in 35 mL sterile phosphate buffer and transferred to a sterile Petri dish filled with the sterile phosphate buffer. Roots were gently washed using sterile tweezers until soil particles were removed. After that, the fresh weight of roots and shoots was measured and they were immediately frozen in liquid N. Soil taken from roots in the phosphate buffer was centrifuged at 3000 rpm for 20 min. Finally, supernatant was removed, and soil was kept as a rhizosphere soil sample and frozen in liquid N. All the samples were stored in −80°C freezer until the next step.

Three out of five replicates were used for each analysis, because the amount of root material was sometimes insufficient for both SL and DNA extraction on the same sample. Therefore, some replicates were used for either SL or DNA extraction.

Strigolactone analysis

Three types of SLs produced by rice, orobanchol, 4-deoxyorobanchol (4DO, also known as 2′-epi-5-deoxystrigol) and methoxy-5-deoxystrigol isomers 1 (MeO5DS), were extracted and identified as previously described from root tissue with small modifications (Jamil et al. 2012); the other two reported isomers of methoxy-5-deoxystrigol were not detected. The frozen roots were ground in a cold mortar with liquid N. About 100–300 mg of ground roots were extracted with precooled 1 mL ethyl acetate. After vortexing, 10 µL of 5 × 10^–7 M GR24 and 10^–7 M [²H₆]-5-deoxystrigol in 50% acetonitrile were added as internal standard (IS). GR24 was used for the semi-quantification of orobanchol and MeO5D, and [²H₆]-5-deoxystrigol was used for the semi-quantification of 4DO. The samples were sonicated for 10 min in ultrasonic bath (Elmasonic S30, Singen, Germany) and centrifuged for 10 min at 2000 rpm. Supernatant was kept temporarily, and the same root sample went through the same extraction procedure again as described earlier. The supernatant from the second extraction was pooled with the supernatant from the first extraction and dried in a vacuum concentrator (ScanSpeed 40; LaboGene, Lynge, Denmark) at 2000 rpm for 2 h. The dried residue was dissolved in 50 µL of ethyl acetate and 4 mL of hexane and mixed thoroughly by vortexing. Solid-phase extraction (SPE; Phenomenex Strata^® SI-1 Silica, 200 mg/3 mL; Phenomenex, Torrance, United States) cartridges were conditioned on manifolds with 2 mL of ethyl acetate and 4 mL of hexane in order, and the samples were loaded through SPE cartridges. After that, SPE were washed with 20 mL of hexane and eluted with 3 mL of ethyl acetate/hexane (9:1) mixture. The eluted samples were dried in a vacuum concentrator at 2000 rpm for 2 h. The residue of the sample was reconstituted in 100 µL of 25% acetonitrile and filtered through micro-centrifugal nylon filter (0.2 µm; Thermo Fisher Scientific, Waltham, United States) in a centrifuge for 1 min at 12 000 rpm.

The Acquity Ultra High Performance Liquid Chromatography (UPLC) I-Class System (Waters, Milford, United States) coupled to Xevo^® TQ-XS tandem quadrupole mass spectrometer (Waters, Milford, United States) with electrospray ionization interface was employed to SL analysis, as previously described (Floková et al. 2020). The separation was performed on Acquity UPLC BEH C18 column (2.1 mm × 100 mm,  1.7 μm particle size; Waters, Milford, United States), maintained at temperature 45°C. Analytes were eluted with 12 min binary gradient, consisted of water (A) and acetonitrile (B) with 15 mM formic acid as an eluent additive. The elution profile was as follows: 0–0.4 min isocratic elution at 15% B; linear gradient increase to 27%, 40% and 65% B in 0.65, 5 and 8 min, respectively; and isocratic elution at 65% until 8.7 min. The column was washed with 95% B for 1.6 min and equilibrated to initial conditions for 1.7 min. The eluate was introduced in ionization source of mass spectrometer, operating in positive mode at following conditions: capillary voltage 1.2 kV; ion source/desolvation temperature 120/550°C; desolvation/cone gas flow 1000/150L h^–1; and collision gas flow 0.15 mL min^–1. Compounds were analyzed in multiple reaction monitoring mode using characteristic diagnostic and confirming transitions (Table S3, Supporting Information) with optimized cone voltage (20–22 V) and collision energy (18–25 eV). MassLynx (v4.2; Waters, Milford, United States) was used to control the instrument, and acquire and process the MS data. SLs were semi-quantified using peak area and its relation to the peak area of the ISs and expressed per gram of fresh root weight.

Microbiome analysis

DNA from the bulk soil, rhizosphere and roots were extracted following protocol given by Mo Bio PowerSoil DNA isolation kit (Qiagen, Hilden, Germany). Additionally, soil was centrifuged in order to remove extra water before extraction as recommended by the supplier for wet soils. We used 300 mg of grinded frozen roots and 250 mg of soil for DNA extraction.

Primers were selected to amplify V3–V4 region of the 16S rRNA gene (16S: 341F 3′-CCTACGGGNBGCASCAG-5′ and 806R 3′-GGACTACNVGGGTWTCTAAT-5′). For fungi, we used the universal ITS primers ITS1 3′-CTTGGTCATTTAGAGGAAGTAA-5′ and ITS2 3′-GCTGCGTTCTTCATCGATGC-5′. Sequencing was performed by StarSEQ, Mainz, Germany, using Ilumina MiSeq. We obtained sequences reads from 198 demultiplexed samples from the company and trimmed the primer and index sequences using the package cutadapt (v1.9.1) and preprocessed the sequencing data using R package DADA2 (v3.11). According to the sequence quality, the 16S rRNA gene reads were filtered to keep 290 bp for the forward read, 170 bp for the reverse read, less than expected errors = 2, greater than Q score = 2 and minimum length = 50 to avoid poor quality and ambiguous sequences. Similarly, ITS reads were filtered to have less than expected errors = 1, greater than Q score = 2 and greater than length = 50. Chimeras were removed after merging denoised pair-end sequence. Each unique amplicon sequence variant (ASV) was assigned to taxa according to SILVA database (v132, https://zenodo.org/record/1172783#.X7PlK2hKibg) and UNITE (04.02.2020 release, https://unite.ut.ee/repository.php) for 16S rRNA gene and ITS, respectively, and an ASV relative abundance table was created for each data set. Afterward, ASVs belonging to plant host were removed (e.g. mitochondria, chloroplasts).

Data analysis and statistics

For microbiome diversity and composition analyses, we used the R packages phyloseq and vegan. First, alpha-diversity indices, including the number of observed species, Chao1, evenness and Shannon were calculated on raw ASV counts. Afterward, the count table was rarefied according to the smallest sample size (903 for bacteria and 5452 for fungi). The sample named ‘R_J9_e’ was removed before rarefying due to its low sequencing depth. Rarefied ASV counts were used to assess beta diversity. Principal coordinates analysis (PCoA) based on Bray–Curtis distance was performed and permutational multivariate analysis of variance (PERMANOVA) was used to test the effect of the soil type, compartment and plant genotype on the beta diversity of the microbial community. Finally, rarefied ASV counts were summed according to their assigned taxa to compare the phylum/order composition using an in-house script in R (see github repository reference in the ‘Data availability’ section). The relationship between plant total biomass and relative abundance of orders belonging to the Glomeromycota was analyzed by linear modeling with P-value validation by 1000 permutations using the R package lmPerm.

To find correlations between SLs level and microbiome composition, we used the data set obtained from the forest soil because only in the plants grown on that soil we could detect SLs. For following analyses, we used the data from the samples for which we had both SL and microbiome data. We performed three analyses to look for the effect of SLs on (i) alpha diversity (i.e. richness, Chao1, evenness and Shannon) using linear regression modeling, (ii) combinations of ASVs (as a community trait) by constrained PCoA (CAP) based on Bray–Curtis distance together with validation by an ANOVA-like permutation test and (iii) the relative abundance of individual both bacterial and fungal ASVs/genera and individual functions predicted for the bacterial community by negative binomial generalized linear model (nbGLM).

For the first analysis, a linear regression model was created with 1000 permutations using the function ‘lmp’ in R package lmPerm. For the second analysis, CAP analysis based on Bray–Curtis distance was conducted together with ANOVA-like permutation with 1000 permutations using functions ‘capscale’ and ‘anova.cca’ in the R package vegan. Further visualization of CAP analysis was performed using the R package phyloseq. For the last analysis, low abundant counts in all input data were first filtered out prior to nbGLM. More specifically, bacterial/fungal ASVs and enzymes/pathways within a bacterial community present more than two times in at least 40% of the total samples were kept. In case of genus counts, bacteria genera that were less abundant than 1% and 0.5% in the roots and rhizosphere, respectively, were removed. Similarly, fungal genera with less than 0.3% and 0.1% abundance in the roots and rhizosphere, respectively, were eliminated. Then, nbGLM was performed to find all correlations between SL level and ASV/genus count with significance level P < 0.05 and between SL level and predicted bacterial pathway/enzyme count at significance level P < 0.01. To ensure goodness of model fit, only models with Shapiro P > 0.05 on residual, less than three outliers based on Cook's distance score (cutoff level < 0.3) and validated P-value by 1000 permutations were selected using an in-house script in R. Bacterial pathway and enzyme count tables were created using PICRUST2 (v2.4.1) that infers putative functions for the bacterial community using the 16S marker gene based on a genome database obtained from environmental samples including soils (Douglas et al. 2020). Subsequently, the count tables were rarefied using the R package phyloseq.

All comparison tests of alpha-diversity indices, biomass and SL levels between experimental conditions were performed using the nonparametric Kruskal–Wallis test and post-hoc Dunn test using R package FSA, rcompanion and multcompView and false discovery rate was corrected by Benjamini–Hochberg (adj. P). Principal component analysis (PCA) was performed on three SLs using the R function ‘prcomp’ with the parameter ‘scale = T’. In general, data was handled using R packages dplyr, tibble, tidyr, stringr and reshape2. Box plot, stack bar plot, regression plot and PCA plot were made using R package ggplot2. The Venn diagram was made by the R package ggvenn.

Data availability

The raw sequence data have been submitted to public data repository Zenodo (DOI: 10.5281/zenodo.4 604 914). All data sets, scripts and information of packages for the data analysis performed in this study can be found in: https://github.com/BoraKim2018/2021_Effect_of_strigolactones_on_recruitment_of_the_rice_root_associated_microbiome.

Results

Plant growth and strigolactone production

The 16 rice genotypes were grown on two different soils, agricultural and forest soil, for a period of 32 days after which plant biomass and SL production were measured (Figs S1 and S2, Supporting Information). Total fresh biomass was significantly affected by soil type (Kruskal–Wallis test, adj. P < 0.001; Fig. S2A, Supporting Information). Plants grown on the agricultural soil had about 1.25 to 5.3 times higher total biomass than plants grown on forest soil. There was a significant genotype effect on biomass for both soils (Dunn test, adj. P < 0.001; Fig. S2B and C, Supporting Information), with d3 and Kinko having the lowest and IAC165, IAC1246, GWD and SC a relatively higher biomass than others, on the agricultural soil. The variation in biomass among genotypes grown on the forest soil was smaller than for the agricultural soil.

As a result, probably, of the higher phosphate level of the agricultural soil (Table S1, Supporting Information), SL levels in rice roots grown on this soil were below the detection limit. In the roots of plants grown on the forest soil, however, we detected three endogenous SLs that have been reported in rice before: orobanchol, 4-deoxyorobanchol (4DO) and methoxy-5-deoxystrigol isomer 1 (MeO5DS), a (putative) SL with unknown structure (Jamil et al. 2012; Fig. 1). Of the 16 genotypes studied, SLs were detected in root tissue of nine, while they were under the detection limit in four low producing varieties (i.e. Kinko, TN1, Bhas and SC) and in the three mutants impaired in SL biosynthesis (i.e.d10, d17 and d27). The amount of SLs significantly varied among rice genotypes (Dunn test, adj. P < 0.05; Fig. 1) and the highest levels of SLs were observed in IAC165, IAC1246, Shiokari and the SL perception mutant d14. PCA showed that the level of the three SLs correlates quite well, as 78% of the total variance was explained by PC1, but the loading of the SL (arrows) also showed that there is variation in the ratio between the SLs, and 4DO and MeO5DS were closer related to each other than to orobanchol as visible on PC2 (Fig. 1D).

Figure 1.

Strigolactone concentration in the roots of rice plants grown on forest soil. SLs—orobanchol (A), 4-deoxyorobanchol (4DO) (B) and methoxy-5-deoxystrigol isomer 1 (Me5DS1) (C)—were analyzed using liquid chromatography–mass spectrometry. Each dot represents a single sample and letters indicate significant differences of means between samples (corrected P < 0.05 for multiple testing, Dunn's test). (D) PCA of SLs. Each dot represents a sample point and the arrows represent the loadings of the SLs. n.d, not detected; GWD, Gangweondo; Bina, Binagimbing; Sonk, Sonkanoir; Bhas, Bhasmanik; SC, Shuang-Chiang.

Bacterial and fungal community diversity and composition

The bacterial and fungal communities in the rhizosphere and roots of the 16 rice genotypes were analyzed using 16S rRNA gene and ITS amplicon sequencing. All samples had sufficient sequencing depth as the number of ASVs was saturated for each sample in the rarefaction curve (Fig. S3, Supporting Information). Bacterial alpha-diversity indices were higher in agricultural soil than forest soil (Kruskal–Wallis test, adj. P < 0.001; Fig. 2A; Table S4, Supporting Information). Furthermore, bacterial diversity in the roots was lower than in the rhizosphere and in the bulk soil, and similar in the rhizosphere and bulk soil (Dunn test, adj. P < 0.05; Fig. 2A; Table S5, Supporting Information). For fungi, the evenness index was not different between the two soils, but the richness, Shannon and Chao1 indices were significantly higher in the forest than in the agricultural soil (Fig. 2B; Table S4, Supporting Information). Additionally, a significant lower fungal alpha diversity was observed in the root compartment than in the rhizosphere and the bulk soil in both soils (Dunn test, adj. P < 0.001, Fig. 2B; Table S5, Supporting Information).

Figure 2.

Diversity and composition of the rice rhizosphere and roots bacterial (A, CandE) and fungal (B, D and F) communities. (A and B) Shannon index. Each dot represents a sample point. Different letters indicate significant differences between samples (corrected P < 0.05 for multiple testing, Dunn's test). (C and D) Principal coordinate ordination based on Bray–Curtis distance. Each dot represents a sample point. (E and F) Microbial composition at phylum level. Ag, agricultural soil; Fo, forest soil; BK, bulk soil; RS, rhizosphere; RT, root.

Beta diversity based on Bray–Curtis distance matrix was analyzed using PCoA (Fig. 2C and D). As seen on the first and second component of the PCoA, the soil type and compartment had an effect on the composition of both bacterial and fungal communities. Subsequently, PERMANOVA showed that the total variation in the microbial community was explained first by the soil type (43% and 45% for the bacterial and fungal communities, respectively; P < 0.001 for both), second by the compartment (11% and 15%, respectively, for the bacterial and fungal communities; P < 0.001 for both) and lastly by plant genotype (4%, for both; P < 0.01 for the bacterial community; P < 0.001 for the fungal community; Table S6, Supporting Information). To assess the influence of the genotype on the composition of the microbial community without the soil and compartment effect, PERMANOVA was employed in each subset of soil and compartment (i.e. rhizosphere in agricultural soil, roots in agricultural soil, rhizosphere in forest soil and root in forest soil). For each of the four data subsets, the genotype explained 35 to 48% of the total variance for most of the models at P < 0.05, except for the fungal community in the rhizosphere of the forest soil (Table S7, Supporting Information).

Also, at the phylum level, there was considerable variation in the bacterial and fungal community composition among soil types and compartments (Fig. 2E and F). In both soils the communities were dominated by Proteobacteria (42.7% and 60.6% in the agricultural and the forest soil, respectively), Bacteroidetes (26.5% and 5.2%, respectively), Verrucomicrobia (12% and 18.7%, respectively) and Acidobacteria (8.5% and 7.1%, respectively) for the bacterial community, and likewise, for the fungal community by Ascomycota (63.7% and 41.9%, respectively),Basidiomycota (8.3% and 32.2%, respectively) and Mortierellomycota (11.4% and 7.3%, respectively). Moreover, Glomeromycota represented 3.2 and 2.5% of the fungal community in the agricultural and the forest soil, respectively. In the Proteobacteria, the Gammaproteobacteria, Alphaproteobacteria and Deltaproteobacteria represented 62%, 28.5% and 9.4% of the sequences, respectively, in the agricultural soil, and 73.7%, 24.2% and 2.1%, respectively, in the forest soil.

Furthermore, we observed an effect of the compartment on the relative abundances of the Glomeromycota (i.e. AMF) that were 9.1 and 11.9 times higher in the root compartment than in the rhizosphere of the agricultural soil and forest soil, respectively, and 51.2 and 42.4 times higher than in the bulk soil, respectively (Dunn test, adj. P < 0.05; Fig. 3A). Although the Glomeromycota were observed to have similar relative abundance in both soils, their composition differed at the order level (Fig. 3B). Glomerales dominated the Glomeromycota community in the agricultural soil, whereas Archaeosporales dominated in the forest soil. In addition, the relative amount of Diversisporales decreased in the presence of the host in the agricultural soil but increased in the forest soil. As AMF are known to help plants to uptake nutrients and thereby increase their biomass on poor soils, we performed linear regression modeling to check correlations between biomass and the abundance of orders belonging to the Glomeromycota. Interestingly, the relative abundances of Glomerales and the Disversisporales did not have any significant correlation with the total biomass of plants grown on the agricultural soil, while the relative abundance of Archaeosporales in both the roots and the rhizosphere and the relative abundance of Diversisporales in the roots displayed a significant positive correlation with the plant biomass on the forest soil (Fig. 3C–E).

Figure 3.

Glomeromycota in the rhizosphere and roots of rice grown on two natural soils. (A) Relative abundance of Glomeromycota. Each dot represents a sample point. Different letters indicate significant differences between samples (corrected P < 0.05 for multiple testing, Dunn's test). (B) Composition of orders belonging to the Glomeromycota. Significant correlation between the plant total fresh biomass and the relative abundance of Glomeromycota orders in (C) the rhizosphere and (D and E) the roots. All correlation studies were performed using linear regression modeling (red line) as indicated by the P-value with a validation by permutation test. Each dot represents a sample point. Ag, agricultural soil; Fo, forest soil; BK, bulk soil; RS, rhizosphere; RT, root.

Effect of SLs on the microbial community in rhizosphere and root

Besides AMF, there might be other microbes that respond to SLs as a recruitment signal. We applied three strategies to identify relationships between the level of SLs and the microbial community diversity and composition by testing the effect of SLs levels on: (i) alpha-diversity indices (i.e. richness, Chao1, evenness and Shannon) using linear regression modeling, (ii) combinations of ASVs, as a community level trait, using CAP analysis based on Bray–Curtis distance and (iii) the relative abundance of individual ASVs/genera using nbGLM. As SLs were under our detection threshold in plant roots grown in agricultural soil, we performed these analyses with the data obtained from the forest soil only.

• Linear regression showed a positive correlation of the bacterial evenness in the root compartment (Fig. S4A, Supporting Information) and a negative correlation of fungal Shannon index in the rhizosphere (Fig. S4B, Supporting Information) with the amount of orobanchol, but not with the other SLs.

• CAP showed that the level of orobanchol significantly correlated with a number of ASVs, but not the other two SLs (Fig. 4). For the bacterial community, the CAP model showed there is a significant orobanchol-associated community in both the rhizosphere and root compartment (Fig. 4; Table S8, Supporting Information). More specifically, in the rhizosphere several ASVs—annotated as Burkholderia–Caballeronia–Paraburkholderia (BCP; bASV25, bASV3), Granulicella (bASV46, bASV31) and Dyella (bASV47, bASV10)—positively associated with the orobanchol level, while Chitinophagaceae (bASV9) and Saccharimonadales (bASV18) negatively associated. Similarly, in the root compartment several ASVs belonging to BCP (bASV3, bASV76, bASV85, bASV25, bASV30) positively associated with the orobanchol level, while Rhodanobacter (bASV24, bASV6, bASV78) both positively and negatively associated as observed on the first axis of the CAP constrained by orobanchol. In addition, we detected a significant effect of orobanchol on the fungal community in both the rhizosphere and root compartment (Fig. 4; Table S9, Supporting Information). In the rhizosphere, Basidiomycota (fASV2, fASV16), Fusarium acutatum (fASV1) and Ascomycota (fASV9) displayed a positive association with orobanchol, while Mortierella (fASV8, fASV20) was negatively associated. In the root compartment, orobanchol was positively associated with Basidiomycota (fASV2, fASV16), Archaeosporales (fASV34) and Chrysozymaceae (fASV75), and negatively associated with Mortierella gemmifera (fASV8), Sordariomycetes (fASV31) and Clonostachys candelabrum (fASV6).

• nbGLM analysis was used to find potential individual candidate ASVs in both the rhizosphere and root compartment correlating with the level of each SL. Including bacteria and fungi, in total five and three ASVs positively correlated with the level of SLs, and 13 and 8 ASVs negatively correlated in the rhizosphere and roots, respectively (Fig. 5A and B; Table S10, Supporting Information). Intriguingly, we found quite different relationships for the individual SLs and there were only few overlapping ASVs between the SLs. For example, in the rhizosphere three significantly correlating ASVs were in common between orobanchol and MeO5DS, while no significantly correlating ASVs were shared between orobanchol and 4DO (Fig. 5A). Just as observed for the CAP described above, two ASVs belonging to Mortierella (fASV8 and fASV20) were negatively associated with orobanchol, and the one ASV belonging to the AMF, fASV34 (Archaeosporales), displayed a positive correlation with orobanchol in both rhizosphere and root compartment (Fig. 5A and B; Table S10, Supporting Information).

Even though we found that a number of ASVs belonging to the genus BCP positively correlated with orobanchol level in the results of CAP (Table S8, Supporting Information), they did not show a significant correlation in nbGLM analysis. We hypothesized that this may be due to the fact that BCP responds to SLs at the genus level but not necessarily at the individual ASV level. Therefore, we created a sum count for each genus—except for the Archaeoporales as most of them were annotated at the order level—and performed nbGLM analysis again. This showed that BCP did display a significant positive correlation with orobanchol, as well as Granulicella and the Archaeoporales, while Mortierella, Solicoccozyma, Saitozyma and Discosia displayed a negative correlation with orobanchol and MeO5DS in the rhizosphere (Fig. 5C; Table S11, Supporting Information). In the roots, the order Archaeoporales and the genus Acaulospora, which both belong to the Glomeromycota, positively correlated with orobanchol (Fig. 5D). 4DO displayed a relationship with different microbes than orobanchol, for example, a negative correlation with Clonostachys in the rhizosphere and a positive correlation with Dyella and Umbelopsis in the root compartment (Table S11, Supporting Information).

Combining the results of CAP and nbGLM analysis suggests that a number of microbial taxa, in both the bacterial and fungal community, are significantly associated with the level of specific SLs. Microbiome components that are consistently positively associated with SLs are BCP, Granulicella and Archaeoporales for orobanchol and Umbelopsis for 4DO. Consistent negative associations of orobanchol were detected with Mortierella, Solicoccozyma, Saitozyma, and Discosia and those of 4DO with Clonostachys.

Orobanchol-associated bacterial functions

As SLs are released under nutrient deficiency as a cry for help, we decided to analyze the putative functions of the bacteria attracted by SLs. Hereto, we predicted bacterial functions based on their 16S rRNA gene sequence using PICRUSt2 (Douglas et al. 2020) and obtained the relative abundance of pathways/enzymes as an output. Then, we chose orobanchol—as the strongest predictor of the presence of AMF and a number of other ASVs/bacteria—to perform nbGLM on its level and the relative abundance of pathways/enzymes.

In the rhizosphere, we found 23 pathways that positively associated with the level of orobanchol, 12 of which are involved in aromatic compound degradation: eight catechol, two toluene, one cinnamate and one 4-hydroxyphenylacetate degradation pathways (Table S12, Supporting Information). Similarly, five out of the eight high-orobanchol enriched pathways in the root compartment were two catechol, one toluene, one aminophenol and one nicotinate degradation pathway (Table S12, Supporting Information).

For a more targeted approach we focused on enzymes that have been reported in the literature to play a role in plant-microbe interaction and are potentially beneficial for plants in a poor soil environment (Table S13, Supporting Information). In the rhizosphere, several enzymes involved in solubilizing phosphate (Pi)/cycling P displayed a positive correlation with the level of orobanchol that are, involved in the production of gluconic acid, pyrroloquinoline-quinone (PQQ) synthase (EC 1.3.3.11), PQQ-dependent glucose 1-dehydrogenase (EC 1.1.5.2), and cellobionic acid phosphorylase (EC 2.4.1.321). However, enzymes producing other organic acids such as oxalic acid, citric acid and tartaric acid did not display a correlation with orobanchol, but fumarate hydratase (EC 4.2.1.2) and malate dehydrogenase (EC 1.1.1.37), involved in the production of malic acid, displayed a negative correlation with orobanchol. Furthermore, a number of phosphatases that are involved in the solubilization of inorganic Pi such as inorganic diphosphatase (EC 3.6.1.1) and exopolyphosphatase (EC 3.6.1.11) also showed a negative correlation with orobanchol, while phytase (EC 3.1.3.8 and EC 3.1.3.26) that mineralizes organic P did not display a correlation with orobanchol. Also, enzymes in N cycling and production of volatiles and siderophores did not display a correlation with orobanchol. 1-Aminocyclopropane-1-carboxylate deaminase (ACC deaminase, EC 3.5.99.7)—that can degrade the ethylene precursor ACC—and several enzymes involved in quorum sensing/biofilm/secretion system—acyl-homoserine-lactone (AHL) synthase (EC 2.3.1.184), extracellular serine protease (EC 3.4.21), mannose-6-phosphate isomerase (EC 5.3.1.8), UDP-forming cellulose synthase (EC 2.4.1.12), phosphogluconate dehydratase (EC 4.3.2.2), and peptidoglycan muramidase (EC 3.2.1.17)—displayed a positive correlation with orobanchol in the rhizosphere. In the same category in the root compartment, only cyclic-guanylate-specific phosphodiesterase (EC 3.1.4.52) displayed a positive association with orobanchol, but several enzymes were negatively correlated; diaminohydroxyphosphoribosylaminopyrimidine deaminase (EC 3.5.4.26), 5-amino-6-(5-phosphoribosylamino) uracil reductase (EC 1.1.1.193), serine/threonine-protein kinase (EC 2.7.11.1), glycogen phosphorylase (EC 2.4.1.1) and chemotaxis protein methyltransferase (EC 2.1.1.80).

Discussion

So far, more than 35 different SLs have been reported in nature. This structural diversity most likely is the result of selective pressure on specificity in signaling, possibly in multiple relationships with microorganisms including as yet undiscovered ones. Here, we provide new information about candidate microorganisms for such relations, as their occurrence in the rhizosphere and/or roots displayed a relationship with the level of SLs produced by the rice host. Moreover, we show that structural differences affect this relationship, with orobanchol displaying the strongest correlations with the relative abundance of soil microbes. Indeed, there was a significant relationship between the level of orobanchol and the diversity and composition of the microbial community in both the rhizosphere and root compartments (Fig. S4, Supporting Information; Fig. 4). Furthermore, orobanchol associated with a larger number of candidates than the other SLs and displayed a significant correlation with AMF (Fig. 5; Tables S10 and S11, Supporting Information). These results suggest that orobanchol may have a stronger impact on the microbial community composition than the other two SLs, even though its level in the roots was about 10-fold lower than that of the other two.

Figure 4.

Principal coordinate ordination based on Bray–Curtis distance constrained by the level of orobanchol in rice grown on forest soil. (AandB) Bacterial community. (CandD) Fungal community. (A and C) Rhizosphere. (B and D) Roots. Arrows present ASVs having high species score extracted from first axis of CAP1. The models were validated using an ANOVA-like permutation test (1000 permutations) as indicated by the P-value. Each dot represents individual sample points and colors represent amount of orobanchol.

Figure 5.

Correlation between the level of SLs and the relative abundance of ASVs/genera in the rhizosphere and roots of rice grown on forest soil. Venn diagram of bacterial and fungal ASVs positively (red label) and negatively (blue label) correlated with the level of SLs (A) in the rhizosphere and (B) in the roots. Family and genus belonging to AMF (Archaeosporales and Acaulospora) and bacterial genera (BCP and Granulicella) positively correlated with the level of orobanchol (C) in the rhizosphere and (D) in the roots. All correlation studies were performed using negative binomial generalized linear modeling (red line) as indicated by the P-value with a validation by permutation test (1000). Each dot represents a sample point. 4DO, 4-deoxyorobanchol; MeO5DS, methoxy-5-deoxystrigol isomer 1; BCP, Burkholderia–Caballeronia–Paraburkholderia.

SLs affect the root-associated microbiome

In this study, we show that the levels of SL have an effect on the diversity and composition of the rice root-associated bacterial and fungal communities. This effect will also depend on the initial microbial inoculum (in the soil), the host species/genotype and the stresses the host is exposed to. So far, three studies have reported that changing the SL pathway alters the microbiome composition in the rhizosphere. When the A. thaliana SL mutant, max4, was grown on a low-P sand–peat moss mixture, its fungal community composition differed from WT, but not the bacterial community (Carvalhais et al. 2019). In rice, however, there were clear differences in both fungal and bacterial communities when comparing WT to d17 and d14 mutants in biosynthesis and perception, respectively (Nasir et al. 2019), despite the fact that plants were grown on a natural soil with sufficient P, which is known to suppress SL exudation. Finally, overexpression of MAX1, D14 and MAX2 in soybean greatly affected the composition of the bacterial community but not the fungal community including AMF Glomeraceae (Liu et al. 2020). In the latter study, plants were grown on natural soil, but physiochemical properties of the soil were not provided. None of these three studies measured SL production and therefore it is unclear whether SLs were exuded and differed between the genotypes investigated. While we showed that the ‘genotype’ factor explained 34–48% of the total variance in the microbial composition, in accordance with these three previous studies (16–40% of the total variance), measuring SL levels allowed us to show that the amount of orobanchol explained 4.3–6.9% of the total variance in the microbial composition (Fig. 4). This may be an underestimation due to inaccuracies in the analysis of SL levels but does suggest that SLs are not the only host trait affecting the microbiome. Studying the microbiome using mutants not completely devoid of SLs but with different SL composition (that are not available yet) is the next step that needs to be made to further demonstrate the importance of individual SLs as a rhizosphere signaling molecule.

Structural specificity of SLs in the relationship with the root-associated microbiome

In a study on the specificity in the perception of SLs by AMF, Akiyama et al. (2010) showed that orobanchol was one of the most active inducers of hyphal branching in G. margarita. They assessed the hyphal branching-inducing activity of 44 SLs and showed that the minimum effective concentration of orobanchol for hyphal branching activity is one of the lowest and about 30-fold lower than that for 4DO. This fits with the results of the present study in which we could barely detect a correlation between 4DO (and MeO5DS) and AMF colonization, while we did find a significant correlation with orobanchol, for both the rhizosphere and the root compartment (Fig. 5C and D).

Although SLs were undetectable in the roots grown on agricultural soil, native AMF did colonize rice roots regardless of the soil type (Fig. 4). This also suggests that plants grown on the agricultural soil might produce SLs but in a concentration under the detection limit of our Liquid Chromatography Triple Quadrupole Mass Spectrometry (1.25 fmol for orobanchol and 0.125 fmol for 4DO; Floková et al. 2020). Nevertheless, plant biomass correlated positively with the relative amount of AMF only in the forest soil (Fig. 4 C–E). Indeed, several studies showed that the positive effect of AMF depends on the host species/genotype, the AMF species and environmental factors, such as nutrient availability (Nagy et al. 2009, Hoeksema et al. 2010, Pozo et al. 2015, Campos et al. 2018). A lack of nutrients in soil (especially P) is the most important factor explaining a positive effect of AMF on biomass (Hoeksema et al. 2010), suggesting that plants reject help from AMF when they have enough nutrients. This likely explains that in the present study rice plants grown on forest soil benefited more from AMF than those on the agricultural soil. To further underpin the effects on AMF, it would be good to use primers more suitable for these fungi (Andreo-Jimenez et al. 2019) as the ITS primers used in our study potentially underestimate the AMF and hence their relationship with SLs.

An intriguing question is why rice is also exuding 4DO and MeO5DS, if orobanchol is such an efficient rhizosphere signaling molecule? Interestingly, we found quite different correlations, with 4DO correlating with the genera Clonostachys, Dyella and Umbelopsis, and orobanchol correlating with BCP,Granuciella and Archaeoporales (Fig. 5; Table S11, Supporting Information). Possibly, the production of both 4DO and orobanchol evolved to communicate with different microbes. For MeO5DS, this is less clear as almost all genus/ASV candidates that correlate with MeO5DS overlapped with those of orobanchol and 4DO. More rigorous tests may be needed to uncover the role of MeO5DS such as with mutants not making this SL anymore. At present, that is not feasible as the biosynthesis and structure of MeO5DS have not yet been elucidated.

Do SLs recruit specific microbial functions?

Rather than recruiting specific microbial species, SLs could also be recruiting specific microbial functions, which could be present in different microbial species. Indeed, in the rhizosphere compartment, we found a positive correlation between orobanchol and several predicted enzymes (i.e. AHL synthase, cellulose synthase and ACC deaminase) known to function in plant–microbe interaction and potentially beneficial for plants (Table S13, Supporting Information). AHL synthase and cellulose synthase are key enzymes in bacterial quorum sensing, which regulates population density, biofilm formation and mobility, and are required for colonization of the plant roots (Zúñiga et al. 2013). Bacterial ACC deaminase reduces plant stress by lowering the level of ethylene in plants (Bulgarelli et al. 2013) and could help rice plants growing on the poor soil in our study. As SLs are actively released from the roots under P deficiency, we anticipated that we may find a relationship with enzymes responsible for Pi solubilization. Indeed, orobanchol positively correlated with gluconic acid-producing enzymes, while it displayed a negative or no correlation with phosphatase, phytase and other organic acid-producing enzymes. This could mean that orobanchol specifically recruits a set of Pi-solubilizing bacteria that produce gluconic acid, rather than recruiting Pi-solubilizing bacteria in a general manner. However, since these putative genes were inferred based on bacterial taxa through bioinformatic analysis, further validation will be needed to confirm the association between orobanchol and these genes.

Orobanchol displays correlations with plant growth-promoting bacteria

Two bacterial genera, BCP and Granuciella, belonging to Acidobacteria subdivision 1 showed a positive correlation with orobanchol (Fig. 5C; Table S11, Supporting Information). Plant-associated Burkholderia (Suárez-Moreno et al. 2012, Kuramae et al. 2020) and Acidobacteria subdivision 1 are mostly plant growth-promoting bacteria (Suárez-Moreno et al. 2012, Kielak et al. 2016). Burkholderia (recently classified as Paraburkholderia) is often detected in poor soils, residing in the rhizosphere rather than in the roots, and has been suggested to play an important role as nutrient solubilizer and biofilm producer (Da Costa et al. 2014, Kuramae et al. 2020). In our study, the positive relationship between orobanchol and BCP was also found only in the rhizosphere.

Whether orobanchol is really the signaling molecule that is recruiting BCP and Granuciella under nutrient-poor environmental conditions cannot be unequivocally concluded from the present study and there are several possible explanations for the statistical correlations that we find. A first explanation could be that orobanchol is indeed a recruitment signal for these bacteria. A recent study showed that the Burkholderiaceae and unclassified Acidobacteria GP1 (subdivision 1) were more abundant in sorghum cultivar SRN-39, which exudes 300–1100-fold more orobanchol than six other sorghum genotypes (Schlemper et al. 2017). Further study will be needed to confirm this, for example by testing the effect of SLs on growth or chemotaxis of candidate isolates and comparing colonization of the roots of WT and SL biosynthesis mutants, preferably also mutants with differences in SL composition (Wang and Bouwmeester 2018). A second explanation could be that these microbes were recruited by other metabolites that are closely linked to SLs. Transcriptome studies revealed a correlation between the SL pathway and phenolic compound biosynthesis, such as phenylalanine metabolism and phenylpropanoid biosynthesis (Nasir et al. 2019, Wang et al. 2021). Indeed, we found a positive correlation between orobanchol and predicted microbial enzyme functions in aromatic compound degradation, especially catechol-related degradation pathways (Table S12, Supporting Information). Catechol is known to decrease microbial respiration due to its toxicity (Zwetsloot et al. 2018), suggesting that microbes positively correlated with orobanchol are exploiting a niche in the rhizosphere through their capacity to detoxify phenolics. To further unravel the link between SLs and other metabolites in the root exudate, such as the phenolics, under Pi starvation and how they affect the microbial community, further study of the root exudate is required, using metabolomics. A third explanation could be that colonization by AMF, for which we have strong indications that it is driven by orobanchol, also causes an increase in other microorganisms, such as BCP, in the rhizosphere. In fact, hyphae of AMF also release compounds into the soil, just as the plant roots do, and recruit their own distinct bacterial community near the hyphae that has been suggested to benefit the plant host as well (Artursson et al. 2006, Toljander et al. 2007). Burkholderia is one of the common AMF-associated bacteria found in the hyphosphere and on the surface of spores (Levy et al. 2009, Scheublin et al. 2010, Agnolucci et al. 2015).

Conclusion

In summary, using data integration on the level of SLs produced by 16 rice genotypes and their root-associated microbial communities, we provide evidence that SLs influence the diversity and composition of the root and rhizosphere microbial communities. We could pinpoint associations between the level of SLs and specific microbial taxa and putative bacterial functions. Furthermore, we observed that structurally different SLs have different effects, with orobanchol displaying the strongest correlations with the relative abundance of soil microbes, including AMF. This insight into the structural specificity of SLs in the relationship with root-associated microbiome gives not only a new perspective on the ecological role of SLs as rhizosphere signals but also a new evolutionary perspective on SL diversity in nature. Further research should underpin the mode of action of SLs, on the side of both the transmitter, the host, and the receiver, the microorganisms. Testing isolates of microbes, positively correlating with SLs in the present study, in bioassays with SLs would be one approach to try to confirm our observations. Another approach to disentangle the chemical/microbial network linked to SLs in the rhizosphere is to create subtler mutants with altered SL composition (instead of full knockouts). However, this requires identification of the genes involved in SL diversification. A better understanding of the SL–microbiome network will provide novel tools for the breeding of crop genotypes producing a tailored cocktail of SLs optimal for the recruitment of a beneficial, stress-protective microbiome.

Acknowledgements

We thank R. de Ruiter at Staatsbosbeheer for providing the forest soil, Plant Physiology Group at Wageningen University and Research for providing rice SL mutants seeds, J. van den Hoek for help with sample preparation and A. Chojnacka for a preliminary SL measurement. Publication number 7369 of the Netherlands Institute of Ecology (NIOO-KNAW).

Funding

This work was supported by the European Research Council through ERC Advanced Grant CHEMCOMRHIZO (Chemical communication in the rhizosphere of plants) to HJB (grant number 670211).

Conflict of interest

None declared.

Author notes

Harro J. Bouwmeester and Anouk Zancarini are authors contributed equally to this work.