Introduction

Licorice (Glycyrrhiza) is a genus of herbaceous perennial plants in the family Leguminosae, with approximately 20 species known as popular medicinal plants in arid and semi-arid regions worldwide (Lewis et al., 2005; Xie et al. 2018). Glycyrrhiza uralensis Fisch. is the most common species of licorice and its roots contain two major secondary metabolites, glycyrrhizin and liquiritin, which belong to triterpene saponins and flavonoids, respectively (Hayashi and Sudo 2009). Licorice roots are widely used in the food, beverage, and cosmetics industries (Zhang and Ye 2009). The roots are also used in the pharmaceutical industry for antiviral (human immunodeficiency virus [HIV] and severe acute respiratory syndrome [SARS]) and antiulcer medication, and for regulating blood coagulation and thrombosis (Cinatl et al. 2003; Hosseinzadeh and Nassiri-Asl 2016). Wild licorice have been exhausted and are gradually being replaced by cultivated licorice, which decreases the quality of licorice (Chen et al. 2017).

The growth of plants and production of secondary metabolites are intimately related to microbial communities. Therefore, the dynamics of the microbiome associated with medicinal plants have been explored extensively (Huang et al. 2018; Koberl et al. 2013). There are several studies that have reported the beneficial effects of arbuscular mycorrhizal symbiosis on plant growth and secondary metabolite accumulations in the roots of licorice (Chen et al. 2017; Xie et al. 2018). Other studies have shown that the growth and efficacy of the medicinal herb Danshen (Salvia miltiorrhiza) can be promoted by beneficial plant-associated microbes, such as those of the Pantoea, Pseudomonas, and Sphingomonas genera (Huang et al. 2018). However, to our knowledge, there is no information available on licorice root-associated bacterial communities and their links with secondary metabolites.

The temporal dynamics of soil microbial communities are usually investigated in chronosequences of land restoration in natural conditions (Lauber et al. 2013; Liu et al. 2019) and different developmental stages of crop plants root microbiota (Zhang et al. 2018). A vital indicator for the dynamics of microbial succession is the temporal turnover of the time–decay relationships, which have been used to describe decreasing community similarities over time (Nekola and White 1999). So far, the time–decay turnover patterns have not yet been investigated in root-associated bacterial communities. Additionally, the structural and functional variability in root-associated microbiomes have been studied in some perennial plants, such as ryegrass (Lolium perenne), Drummond’s rockcress (Boechera stricta), alpine rockcress (Arabis alpina), and fruit (Citrus) (Chen et al. 2016; Dombrowski et al. 2017; Wagner et al. 2016; Xu et al. 2018). However, there is currently a lack of research on the temporal dynamics of root-associated bacterial communities in perennial licorice.

Soil represents the most diverse ecosystem on earth with an exceptionally high microbial diversity (Fierer and Jackson 2006). During the growth process, plants select and enrich microorganisms from bulk soil to root-associated compartments (such as the rhizosphere and endosphere), each of which harbors a distinct microbiome (Edwards et al. 2015). The rhizosphere soil is profoundly influenced by the secretion of plant root exudates, which accordingly results in a ‘hotspot’ environment (Peiffer et al. 2013; Schlaeppi et al. 2014). In contrast, the root endosphere features a highly specific microbiome that is mainly affected by plant genotype (Vandenkoornhuyse et al. 2015). Several studies have reported on the root-associated dominant or core bacteria using Arabidopsis thaliana as the model plant (Bulgarelli et al. 2012; Lundberg et al. 2012). Information on the core microbiota is critical for understanding the assembly and stability of root-associated microbial communities. Licorice plants of different ages may be accompanied by consecutively and specifically enriched bacterial communities among distinct soil–root compartments. However, the selective core enrichment process of root-associated bacterial communities in licorice has rarely been investigated.

Although root-associated microbiota have been investigated for a long time, there is still little agreement on how these communities are shaped and what factors regulate their composition (Bulgarelli et al. 2012; Fan et al. 2018). The temporal dynamics of microbial communities in natural environments are influenced by multiple factors simultaneously. Therefore, it is difficult to determine the relative contributions made by individual factors to the overall microbial succession process (Jiao et al. 2017). The assembly and composition of the root-associated microbiome could be influenced by soil conditions such as the pH, nitrogen and phosphorus contents (Edwards et al. 2015; Xiao et al. 2017). Therefore, it is desirable to obtain a predictive understanding of how root-associated bacterial communities in licorice plants of different ages respond to soil characteristics. In this context, variation in soil characteristics including primary nutrients, soil texture, and micronutrients, is essential for normal growth and secondary metabolite production in plants (Fageria et al. 2002). Therefore, understanding how soil variables regulate the individual taxa of the microbial communities associated with the accumulation of root secondary metabolites in licorice is crucial for the optimization of plant–soil interactions for practical licorice cultivation.

Many studies have analyzed soil microbial communities using Illumina paired-end sequencing of 16S rRNA and internal transcribed spacer (ITS) amplicons (Jiao et al. 2017; Xiao et al. 2017). In this study, single-end sequencing was adopted using the IonS5™XL platform. This provided longer sequences to read without splicing (Mehrotra et al. 2017), and thus could help to elucidate the root-associated bacterial communities more comprehensively and accurately. The aims of this study were to (1) elucidate the temporal dynamics of root-associated bacterial communities together with variation in soil characteristics and secondary metabolite concentrations in roots; (2) investigate the core-enriched taxa and their time–decay relationships; and (3) provide a comprehensive understanding of the key factor(s) regulating the temporal dynamics of individual taxa related to root secondary metabolites in licorice.

Materials and methods

Study site

The study area is located in the Yuzhong North Mountains (104°18′–104°19′ E, 36°8′–36°9’ N), in Lanzhou, Gansu province (~150 km from downtown Lanzhou) in northeast China. This area has an altitude of 2400–2500 m in the mountains and it is characterized by a temperate semi-arid continental climate. It is known as the cradle and cultivation base of licorice (Glycyrrhiza uralensis Fisch.) with wide distributions of licorice plants. Its average annual and accumulated temperatures are 6.7 °C and 2754–3285 °C, respectively. The average annual sunshine duration is 2500–2600 h and average annual precipitation is 400 mm. The soil is classified as loessial (Calcaric Cambisol according to FAO classification). The licorice species in the study area were all identified as G. uralensis Fisch. The general fertilization regime was conducted with nitrogen (urea: 85 kg ha−1), phosphorus (P2O5: 40 kg ha−1), and potassium (KCl: 35 kg ha−1). The cultivated sites were subjected to similar fertilization and management practices, with close distance (< 500 m) from the wild site and each other.

Experimental design and sampling strategy

One-year-old uniform licorice seedlings were transplanted from a nursery into cultivated sites for one to four years until sampling in July 2018 (the transplanting was staggered so that all samples were collected in the same year). Licorice seedlings were properly cultivated after transplantation, but the plants that had been cultivated for three years had been excavated ahead of our sampling. Therefore, the samples selected from the cultivated sites were referred to as 1-y, 2-y, and 4-y (Supplementary Table S1). According to the descriptions of local residents, the wild licorice had been growing for over 10 years.

The sampling comprised five plots as biological replicates following a Z-shaped pattern in each site, where licorice grew evenly. For each plot, five random soil cores from 0 to 20 cm depth and with a 10 cm diameter were collected and mixed to form a single sample as bulk soil. A group of three to five licorice plants were excavated by digging around the group to keep the roots as intact as possible in each plot. All samples were placed into sterile plastic bags, kept on dry ice, and taken to the laboratory immediately. Bulk soils were homogenized by passing the samples through a 2 mm sieve. The soils were then divided into two subsamples. Licorice roots were processed to obtain the rhizosphere soil according to a previously described method (Xiao et al. 2017). The scrubbed roots (the root endosphere) were also divided into two subsamples. Subsample of the bulk soils and the licorice roots were used for the analysis of soil characteristics and root secondary metabolites, respectively. The remaining (sub)samples were stored at −80 °C in preparation for microbial analysis. In total, 60 samples (four sites × three soil–root compartments × five plots) were obtained for further processing.

Analysis of soil characteristics and root secondary metabolites

Bulk soil characteristics were divided into three groups including primary nutrients, micronutrients, and texture. Soil pH was determined using a routine method with a fresh soil to water ratio of 1:5. The soil water content (SWC) was measured gravimetrically by drying 5 g fresh soil until the soil reached a constant weight. For the analysis of soil organic matter (SOM), total carbon (TC), total nitrogen (TN), total phosphorus (TP), and total potassium (TK), the samples were air-dried and sieved (1 mm mesh), and the contents were determined by combustion (Fan et al. 2018). The soil texture was tested using a laser particle sizer (LS13320, Beckman Coulter, USA) with air-dried soil (Fan et al. 2018). For the soil available nitrogen (AN), phosphorus (AP), and potassium (AK) analyses, the samples were dispersed or extracted, and the contents were determined using an atomic absorption spectrometer (Fan et al. 2018). Specifically, the TN was quantified following the Kjeldahl method and the AN was quantified by extraction with 1 mol L−1 KCl, steam distillation, and titration following the Alkaline diffusion method. The TP was extracted using 1 N HCl after ignition at 550 °C, the AP were extracted using 0.5 mol L−1 NaHCO3 (pH = 8.5), and then they were measured following the Mo-Sb colorimetric method. The TK was digested in a nickel crucible with sodium hydroxide at 750 °C, the AK was extracted using 1 mol L−1 ammonium acetate, and then they were measured by using flame photometry.

Soil micronutrients including copper (Cu), iron (Fe), and zinc (Zn) were analyzed using atomic absorption spectroscopy (Fageria et al. 2002). Moreover, three representative root secondary metabolites in licorice were analyzed, namely glycyrrhizin, liquiritin, and isoliquiritin (liquiritin isomer). Chromatography–mass spectrometry was performed according to standard procedure (Zhang and Ye 2009) with slight modifications. The detailed method is available in Supplementary S1.

DNA extraction, high-throughput sequencing, and bioinformatic analysis

The total DNA was extracted from bulk and rhizosphere soil samples (0.5 g each) using the Fast DNA®SPIN Kit (MP Biochemicals, Solon, USA). The total DNA was extracted from the root samples (0.5 g each) using a DNA secure Plant Kit (Tiangen Biotech, Beijing, China) following the manufacturer’s procedures. The DNA concentration and purity were estimated using a Nanodrop 1000 spectrophotometer (Thermofisher Inc., Wilmington, Germany) and electrophoresis in 1% (w/v) agarose gel (Xiao et al. 2017). The hypervariable V4-V5 region of the 16S rRNA gene was amplified using the primers pair 515F (5′-GTG CCA GCM GCC GCG GTA A-3′) and 907R (5′-CCG TCA ATT CCT TTG AGT TT-3′) (Edwards et al. 2015). The polymerase chain reaction (PCR) amplifications were performed in triplicate for each sample, and PCR products were mixed in equidensity ratios. Then, the mixed PCR products were purified using a GeneJET™ Gel Extraction Kit (Thermo Scientific). Sequencing libraries were generated using Ion Plus Fragment Library Kit (48 rxns; Thermo Scientific). The library quality was assessed using the Qubit@ 2.0 Fluorometer (Thermo Scientific). Finally, the library was sequenced on an IonS5™XL platform (Thermofisher Inc., Massachusetts, USA) and 400 bp single-end reads were generated by Novogene (Beijing, China).

Low-quality sequences were sheared using Cutadapt (Martin 2011) and quality-filtered using the QIIME pipeline (v1.7.0) (Caporaso et al. 2010). After the removal of chimeric sequences using the USEARCH tool in the UCHIME algorithm (Edgar et al. 2011), the remaining sequences were assigned to operational taxonomic units (OTUs) at similarities of 97% using the UPARSE pipeline (Edgar et al. 2011). OTUs lacking more than two sequences were removed. Taxonomic information was annotated for a representative sequence of each OTU using the RDP classifier at a confidence level of 80% (Wang et al. 2007) using the SILVA database release 128.

Statistical analysis

All statistical analyses were conducted in the R environment (v3.6.1; http://www.r-project.org/). Most of the results were visualized using the ‘ggplot2’ package (Gómez Rubio 2017), unless otherwise indicated. The soil characteristics and root secondary metabolite concentrations were log-transformed before the significance of factors (plant sites) was tested by analysis of variance (ANOVA), followed by comparisons of means using Tukey HSD parametric tests (‘multcomp’ package) (Hothorn et al. 2008). Principal component analysis (PCA) was used to investigate the distribution of these variables on a scaled parameter matrix and visualized using the ‘FactoMineR’ package (Lê et al. 2008). The Spearman’s correlation of these variables was analyzed and visualized using the ‘corrplot’ package (Wei and Simko 2013).

Alpha-diversity indices were calculated using QIIME (v1.7.0) and significant differences were found using Kruskal-Wallis non-parametric tests along with ANOVA (‘agricolae’ package) (De Mendiburu 2014). Nonmetric multidimensional scaling (NMDS) was implemented, with significant differences in bacterial community composition tested by three different but complementary non-parametric multivariate statistical analysis methods (‘vegan’ package) (Oksanen et al. 2007): permutational multivariate analysis of variance using distance matrices (Adonis), analysis of similarities (Anosim), and multiple response permutation procedure (MRPP). Heatmaps were illustrated based on Z-score-normalized relative abundance of taxa using the ‘pheatmap’ package (Kolde 2015). A circos plot was created based on the relative abundance of taxa using the ‘circlize’ package (Gu et al. 2014). Boxplots were used to illustrate the variation in alpha- and beta-diversity values.

The most common and ubiquitous bacterial taxa across all sites were identified as dominant taxa using the criteria of a recent study (Delgado-Baquerizo et al. 2018) with slight modifications. Analysis of the differential taxa was performed using a negative binomial generalized linear model in the ‘edgeR’ package (Nikolayeva and Robinson 2014). The OTU counts from bulk soil were taken as a control, and corresponding P-values were corrected for multiple tests using a false discovery rate (FDR) set at 0.05 (Benjamini and Hochberg 1995). Core-enriched taxa were selected from overlapping OTUs in Venn diagrams using the ‘Venndiagram’ package (Chen and Boutros 2011). Ternary plots were visualized to show the distribution of these taxa using the ‘ggtern’ package (Hamilton 2015). The functions of these taxa were annotated using FAPROTAX (Louca et al. 2016). Time–decay relationships were estimated by fitting a linear model between the changes in community similarity (based on 1 – [Bray-Curtis dissimilarity]) and the age of licorice plants; the slopes were used to measure the rate of community turnover (Nekola and White 1999).

The correlations among the matrices of the three compartment communities, soil characteristics, and secondary metabolites in the roots were examined using the Mantel test (‘vegan’ package). Multiple regression on distance matrices (MRM) analysis and the random forest (RF) regression algorithm were used to investigate the effect level of significantly correlated variables on bacterial communities (assessed by Bray-Curtis distance and Beta-NMDS1) using the ‘ecodist’ (Goslee and Urban 2007) and ‘rfPermute’ (Jiao et al. 2019) packages, respectively. Linear regression was used to investigate the correlations between log-transformed TK content and the average relative abundance of core-enriched taxa; significant correlations were calculated for each taxon (‘corrplot’ package). Individual taxa were then selected from these significantly correlated taxa again, using the ‘maSigPro’ package (Conesa et al. 2006). These taxa were visualized using curve plots to illustrate their temporal variability; their correlations with secondary metabolites were calculated using the ‘corrplot’ package.

Results

Temporal variation in soil characteristics and root secondary metabolites

The ANOVA results showed that the age of licorice plants had significant effects on most variables (P < 0.05). Most soil characteristics such as SOM, TC, TN, and Cu, significantly increased with increasing age of cultivated plants; the C/N ratio significantly decreased correspondingly (Table 1). Other soil characteristics such as SWC, TP, AP, TK, AK, Zn, silt, and clay also increased, whereas AN, Fe, and sand decreased. In addition, the isoliquiritin concentrations in licorice roots steadily decreased with increasing age of cultivated plants. However, the liquiritin and glycyrrhizin concentrations first decreased and then recovered slightly (Table 2). This trend was similar to the temporal variation in soil AN, TP, and TK levels. Many soil characteristics (SWC, TP, AP, TK, AK, Fe, and silt) occurred at their lowest levels in the wild site that lacked fertilization and agricultural management. However, all three secondary metabolites had higher concentrations in wild than cultivated licorice.

Table 1 Soil characteristics in sites with different licorice plant ages (n = 5); values are mean ± standard deviation
Full size table
Table 2 Concentrations of secondary metabolites in the roots of licorice plants of different ages (n = 5); values are mean ± standard deviation
Full size table

The PCA based on scaled soil characteristics and secondary metabolite concentrations showed that the samples were clearly separated according to the age of licorice plants (Fig. 1A). The contributions of variables to principle components were different and their correlations were determined according to the angles among the variables (Fig. 1B), which were further confirmed by the correlation analysis. The concentrations of the three secondary metabolites were strongly positively correlated with each other, but were significantly negatively correlated with most of the soil primary nutrients (TK, TP, AK, SWC, AP, SOM, TC, and TN). Among these, only TK was related to all three secondary metabolites. Significant correlations were also illustrated between soil characteristics (Supplementary Fig. S1). In general, soil characteristics and root secondary metabolites had varying degrees of temporal variation and exhibited intimate connections.

Fig. 1
figure 1

Principal component analysis (PCA) of soil characteristics and secondary metabolites in the roots of licorice plants of different ages. (a) Ordering distribution and interpretation of the first ten principal components. (b) The contribution of the variables to the principal components. The depth of the color of the arrow line represents the degree of contribution. The angle of the arrow line represents correlations. pH, soil pH; SWC, soil water content; SOM, soil organic matter; TC, total carbon; TN, total nitrogen; AN, available nitrogen; C/N, total carbon/nitrogen ratio; TP, total phosphorus; AP, available phosphorus; TK, total potassium; AK, available potassium; Cu, copper; Fe, iron; Zn, zinc; Liq, liquiritin; Isoliq, isoliquiritin; and Acid, glycyrrhizin

Full size image

Temporal succession of root-associated bacterial communities

Based on the sequencing of 16S rRNA amplicons, 4,464,796 quality reads were obtained (after filtering) from the 60 samples (20 samples per soil–root compartment). The number of reads per sample ranged from 51,305 to 87,560 with an average of 74,413 ± 9175 reads. In total 3,223,764 reads could be classified as bacteria, and most (92.7%) could be classified at the phylum level. A total of 11,430 OTUs were clustered and the number of OTUs per sample ranged from 72 to 3233, with an average of 1626 ± 1087 OTUs. After homogenization, subsequent analysis was carried out based on a minimum of 26,262 reads per sample.

The most of abundant phyla were Oxyphotobacteria (BS: 0.1%; RS: 3.0%; R: 82.6%), Proteobacteria (BS: 21.6%; RS: 60.3%; R: 1.2%), Actinobacteria (BS: 21.7%; RS: 18.9%; R: 0.1%), Acidobacteria (BS: 19.3%; RS: 3.1%; R: 0.01%), Bacteroidetes (BS: 6.9%; RS: 8.5%; R: 0.1%), Planctomycetes (BS: 8.7%; RS: 1.7%; R: 0.004%), Chloroflexi (BS: 6.2%; RS: 1.8%; R: 0.003%), and Gemmatimonadetes (BS: 6.9%; RS: 0.7%; R: 0.002%). Different soil–root compartments had distinct relative abundance of taxa, which also varied across the licorice sites (Supplementary Fig. S2). Overall, the bulk soil communities were dominated by Actinobacteria, Acidobacteria, Planctomycetes, Chloroflexi, Gemmatimonadetes, and Deltaproteobacteria. The rhizosphere communities were dominated by Alphaproteobacteria, Gammaproteobacteria, Bacteroidetes, and Firmicutes. In particular, Oxyphotobacteria accounted for the majority (82.6%) of the root endosphere communities (Fig. 2A).

Fig. 2
figure 2

(a) Circos plot showing the phylum distribution of licorice root-associated bacterial communities in the three soil–root compartments. Because the relative abundance of Proteobacteria was greater than 25%, bacteria in this phylum are shown at the class level. (b) Boxplots of the alpha-diversity (Shannon index and operational taxonomic unit [OTU] richness) of root-associated bacteria across different sites in the three soil–root compartments. (c) Nonmetric multidimensional scaling (NMDS) of root-associated bacterial community composition across different sites in three soil–root compartments based on Bray-Curtis distances. The three insets represent the three soil–root compartments. (d) Boxplots of root-associated bacterial community similarity among the three soil–root compartments. Different letters indicate significant differences (P < 0.05, Kruskal-wallis test). BS, bulk soil; RS, rhizosphere soil; R, root endosphere

Full size image

The alpha diversity of the root-associated bacterial communities was significantly influenced by soil–root compartment and licorice age (Supplementary Table S2). The Shannon index and OTU richness decreased with root proximity, while they increased in the rhizosphere soil communities with increasing age of licorice plants. Similarly, both indices of the bulk soil communities increased from the 1-y to 4-y cultivated licorice (this was not observed for wild licorice). For the root endosphere communities, the indices were not markedly different, apart from the OTU richness in the wild site compared with that of the 1-y and 2-y sites, which significantly differed (Fig. 2B).

In the NMDS plot, different soil–root compartments were clearly separated. The distinct sampling sites were well clustered in the bulk and rhizosphere soils. Specifically, rhizosphere soil communities were completely divided into four clusters corresponding to the three ages of the cultivated licorice and the wild licorice (Fig. 2C). The community similarity was ranked as rhizosphere soil < bulk soil < root endosphere (Fig. 2D). The stress values of NMDS analysis ranged from 0.028 to 0.111 in the whole bacterial communities and the three compartments, which reflected the accuracy of the method. The NMDS results were also confirmed by three multiple statistical approaches (Table 3). For example, compartment and age both had significant effects on the composition of whole communities (P < 0.05). When the whole communities were split into three compartments, the age of licorice plants had greater effects (P < 0.05) on the rhizosphere soil communities (Adonis: R2 = 0.601, Anosim: R = 0.842, MRPP: δ = 0.355) than on the bulk soil communities (R2 = 0.485, R = 0.660, δ = 0.379). However, age had no significant effects on the composition of the root endosphere communities (P > 0.05) with respect to community similarity. Taken together, these results revealed distinct temporal successions of the diversity and composition of the root-associated bacterial communities in the three soil–root compartments with increasing age of licorice plants.

Table 3 Three statistical analyses of the whole root-associated bacterial community composition among different factors contributing to the variation in sequencing data based on Bray-Curtis distances
Full size table

Time–decay relationships of core-enriched taxa

First, the dominant taxa that were highly abundant and ubiquitous in the licorice plants of different ages were screened. The taxa in the top 20.0% of relative abundance (highly abundant taxa) and occurring in more than half of all soil samples (ubiquitous taxa) were collected. This yielded 360 OTUs, accounting for 3.2% of all observed taxa. However, on average, these dominant taxa accounted for 65.3% of the sequences. Additionally, the dominant communities had close relationships (r = 0.98; P < 0.001; Supplementary Fig. S3A) with the remaining communities and exhibited similar distribution patterns with the whole communities (Supplementary Fig. S3B and Table S3).

Second, the number of enriched taxa was analyzed in the rhizosphere soil and root endosphere compared with the bulk soil based on the dominant taxa. In total, 104 (1-y), 110 (2-y), 118 (4-y), and 130 (wild) OTUs were enriched in rhizosphere soil across the different licorice sites. The growing number of enriched taxa revealed an increasing trend with temporal succession. On the contrary, the number of enriched OTUs displayed a decreasing trend in the root endosphere, with 15, 14, 7, and 7 OTUs in the 1-y, 2-y, 4-y, and wild sites, separately. Interestingly, the number of root endosphere-enriched OTUs was reduced by half in the 4-y site and then became stable. Despite these differences, 62 OTUs were found to be overlapped and consecutively enriched in the rhizosphere soil across the four ages of licorice. Similarly, there were 5 OTUs enriched consecutively in the root endosphere (Fig. 3A).

Fig. 3
figure 3

Venn diagrams of the enriched taxa across different licorice sites in the rhizosphere soil (RS), root endosphere (R), and bulk soil (BS). (a) The rhizosphere soil and root endosphere individually and consecutively enriched and depleted OTUs compared with the bulk soil in the corresponding sites. The bulk soil enriched OTUs are depleted OTUs in the rhizosphere soil and root endosphere. (b) Overlapping between the rhizosphere soil and root endosphere consecutively enriched and depleted OTUs. Red font, enriched OTUs; green font, depleted OTUs. The number of individually, consecutively, and core-enriched discriminating OTUs is indicated in brackets. (c) Time–decay relationships of the whole, dominant, and core-enriched taxa. The rate of community turnover across time (slopes) and the fitting degree (R2) of the line are provided. Significance of the linear fitting model is marked with *** (P < 0.001)

Full size image

Third, core-enriched OTUs were selected in the rhizosphere soil (58 OTUs) and root endosphere (1 OTUs) from the overlapping of consecutively enriched OTUs between these two compartments. It was found that more OTUs were consecutively enriched in the bulk soil when compared with the root endosphere (216 OTUs) than when compared with the rhizosphere soil (57 OTUs) across the four licorice ages. Thus, the bulk soil had 57 core-enriched OTUs (Fig. 3B). A substantial number of these OTUs in the bulk and rhizosphere soils belonged to the phyla Acidobacteria and Proteobacteria, respectively (Fig. 4A). The core-enriched OTUs also had distinct distributions and relative abundances in the three compartments across the four ages of licorice (Fig. 4B). The results from FAPROTAX showed that these core-enriched OTUs were related to 17 functional pathways (Fig. 4C). In general, most of the OTUs had chemoheterotrophic features. Some core-enriched OTUs of the bulk soil were related to nitrification and fermentation. In contrast, some core-enriched OTUs of the rhizosphere soil were related to denitrification, plant pathogen, and aromatic compound degradation.

Fig. 4
figure 4

Distribution and functions of root-associated core-enriched taxa in licorice plants of different ages. (a) Taxonomic composition of the core-enriched OTUs in the bulk and rhizosphere soils at the phylum level. (b) The ternary diagrams of core-enriched OTUs across different licorice sites (1-y, 2-y, 4-y and wild) in the three soil–root compartments. (c) Metabolic and ecological functions of core-enriched OTUs in the bulk and rhizosphere soils based on FAPROTAX. Each column represents an OTU. The presence of functions is shown in red color. BS, bulk soil; RS, rhizosphere soil; and R, root endosphere

Full size image

Finally, the time–decay relationships of the whole, dominant, and core-enriched communities were elucidated among the three soil–root compartments (Fig. 3C). The whole and dominant bacterial communities had similar time–decay relationships (P < 0.001) in the bulk and rhizosphere soils. The rhizosphere community displayed a stronger rate of decay (slopes = −0.037) than the bulk soil community (slopes = −0.028). However, the core-enriched community of both the bulk and rhizosphere soils displayed decreased turnover slopes at almost the same rate of decay (slopes = −0.026). In the root endosphere, the whole, dominant, and core-enriched communities had no significant time–decay relationships due to the relatively stable community composition.

Effect of soil potassium on individual taxa

Through visualizing of the Mantel test, it was found that soil TK, TC, TN, C/N, Cu, and clay as well as root isoliquiritin and glycyrrhizin levels, had significant correlations with the core-enriched bacterial communities of the bulk soil. These variables were deemed latent factors affecting these communities. Distinct patterns occurred in the rhizosphere soil and root endosphere (Fig. 5A). Among the latent factors, soil TK had higher correlations with the core-enriched bacterial communities of the bulk and rhizosphere soils than the other factors. Additionally, soil TK consistently displayed a significantly higher explained variance for core-enriched bacterial communities in the two compartments according to the results of MRM (Table 4) and RF (Supplementary Table S4) analyses, which confirmed the results of the Mantel test. Therefore, soil TK could be a key factor regulating the temporal dynamics of core-enriched bacterial communities in the bulk and rhizosphere soils.

Fig. 5
figure 5

(a) Correlation analysis among soil characteristics, root secondary metabolites, and root-associated core-enriched bacterial communities (Bray-Curtis distances) in the three soil–root compartments based on the Mantel test. The color of the line represents the significance of the differences (P-values). The size of the line represents correlation coefficients (Mantel’s r). (b) Correlation heatmaps between TK and the core-enriched OTUs in the bulk and rhizosphere soils. The number of significant correlations is indicated in brackets. (c) Temporal variation curves of individual taxa with significant correlations and temporal variability in the bulk and rhizosphere soils. pH, soil pH; SWC, soil water content; SOM, soil organic matter; TC, total carbon; TN, total nitrogen; AN, available nitrogen; C/N, total carbon/nitrogen ratio; TP, total phosphorus; AP, available phosphorus; TK, total potassium; AK, available potassium; Cu, copper; Fe, iron; Zn, zinc; Liq, liquiritin; Isoliq, isoliquiritin; and Acid, glycyrrhizin. BS, bulk soil; RS, rhizosphere soil; and R, root endosphere

Full size image
Table 4 Results of the multiple regression on distance matrices (MRM) analysis of the core-enriched bacterial communities (Bray–Curtis distances) in the bulk soil (BS) and rhizosphere soil (RS) for each latent factor
Full size table

Two divergent clusters of individuals from the core-enriched taxa were then characterized: negatively (TK-Neg) and positively (TK-Pos) related to TK in the bulk and rhizosphere soils (Fig. 5B). Their characterization was confirmed by the linear regression model (Supplementary Fig. S4). Based on the ‘maSigPro’ analysis, the taxa that were significantly related to soil TK were further selected along with time series. In particular, 11 individual taxa showed significant temporal variation in the bulk and rhizosphere soils with various annotations (Supplementary Table S5). Next, the overall temporal variation in these taxa and their significant correlations with secondary metabolites were illustrated (Fig. 5C; Supplementary Fig. S5). In detail, all TK-Neg taxa showed consistent variation in bulk soil and displayed positive correlation with some of the secondary metabolites. These taxa all belonged to the phyla Acidobacteria and Actinobacteria. However, some taxa (including OTU 20, 43, and 10,491) showed the same trends with temporal variation in soil TK in the rhizosphere soil regardless of their negative correlations. These taxa were affiliated to the phyla Proteobacteria and Actinobacteria, which were further annotated as Inquilinus, Nocardioides, and Promicromonospora at the genus level, respectively. These taxa were all positively correlated with root glycyrrhizin and isoliquiritin, and their highest relative abundances were found in the wild site, which exhibited the lowest soil TK content. Therefore, all TK-Pos taxa displayed their lowest relative abundances in the wild site.

In addition, the temporal variation of some TK-Pos taxa in the bulk soil (including OTU 189, 325, 5144, and 271) opposed the TK variation. Some of these taxa showed negative correlations with the three secondary metabolites; they belonged to the phyla Acidobacteria and Bacteroidetes (in which they were annotated to the class of Ignavibacteria). This result was also found in the rhizosphere soil. OTU 16 and 13 were negatively correlated with the secondary metabolites. They were annotated to the Xanthomonadaceae family and the Variovorax genus, respectively, both of which belong to the phylum Proteobacteria. Collectively, these individual taxa displayed significant temporal variation and habitat adaptation related to soil TK. Additionally, they may play different roles in adjusting the accumulation of root secondary metabolites in licorice, which provides clues for later isolation of these individual taxa.

Discussion

Specific root exudates have been shown to exhibit high selectivity for the soil microbiome in rhizosphere soil (Edwards et al. 2015; Huang et al. 2018). In the current study, the rhizosphere soil communities were clearly discriminated from the bulk soil communities. This was due to the greater effects of licorice age and habitat filtering on the rhizosphere soil. The research focus was narrowed down to a few hundred dominant taxa in the root-associated bacterial communities, and the core-enriched taxa in different soil–root compartments were investigated. Soil TK was defined as a key factor affecting the core-enriched bacterial communities according to multiple but complementary methods. On this basis, we further explored the effects of the key factor on individual taxa associated with the roots of licorice. This approach differed from other studies that have always focused on the community level when investigating the relationships between soil characteristics and microbial communities (Chen et al. 2016; Fan et al. 2018). This analysis facilitated rapid and accurate forecasting of soil characteristics that mediated the links between root secondary metabolites and the individual root-associated taxa (Fig. 6). Consequently, we proposed a focused and novel framework for revealing the regulatory effects of soil characteristics on core-enriched and individual bacterial taxa associated with the synthesis of root secondary metabolites in licorice.

Fig. 6
figure 6

Conceptual diagram of total soil potassium (TK) mediating the links between root secondary metabolites in licorice and individual taxa (OTUs) in the bulk soil (BS) and rhizosphere soil (RS). The color of the line represents significantly positive (blue) and negative (red) correlations. The size of line represents correlation coefficients (r). Liq, liquiritin; Isoliq, isoliquiritin; and Acid, glycyrrhizin. Etabolites in licorice

Full size image

Effects of soil characteristics on root secondary metabolites

Normally, the content of nutrients in the soil decreases during perennial growth due to plant nutrient uptake (Chen et al. 2016; Shakya et al. 2013) and removal in the harvested portion of the crop. In contrast, most of the soil characteristics tested in this study significantly increased as the age of licorice plants increased. For example, the increase in SOM, TN, and TC may be mainly due to the effects of normal irrigation and fertilization every year, but also on account of the accumulation of litters from plant leaf and stem. The SWC displayed an increasing trend in the cultivated sites and was negatively correlated with the concentrations of secondary metabolites across all sites. This may be explained by the fact that licorice is resistant to water deficient environment and widely used for ecological restoration due to its drought-tolerant features (Li et al. 2011; Xie et al. 2018). In addition, there is more accumulation of licorice root secondary metabolites under environmental stress conditions (Xie et al. 2018). The increasing trend in silt and clay contents in the cultivated sites compared with the wild site not only confirmed licorice’s role in soil water conservation and improvement, but identified the impacts of the better field management in the cultivated sites.

The three secondary metabolites had higher values in the wild site, which was not well managed or fertilized. Therefore, licorice plants may be inherently excellent at maintaining higher secondary metabolite concentrations in roots, despite lower levels of SWC and other soil primary nutrients. In the present study, soil TK was negatively related to the glycyrrhizin, liquiritin, and isoliquiritin concentrations in the licorice roots, implying its potential roles in the accumulation of secondary metabolites. Potassium is also involved in many enzyme activation systems in plants, which promotes the stem strength and improves the stress resistance of plants (Wang and Wu 2017). Soil potassium is directly related to the production of phytonutrients in oregano (Origanum vulgare L.) (Jan et al. 2018).

The distributions of soil characteristics and secondary metabolites were divided into four groups corresponding to the composition of bulk soil bacterial communities. This division confirmed that the age of licorice plants had significant effects on some of the soil characteristics and secondary metabolites (as was found from the results of the ANOVA and significance tests). Root secondary metabolites can be massively accumulated in licorice under drought or nutritional stress conditions (Chen et al. 2017; Xie et al. 2018). This is consistent with the results that the root secondary metabolites were all significantly negatively correlated with some soil primary nutrients in our study. Thus, excessive soil nutrients are not necessarily beneficial to the accumulation of root secondary metabolites in licorice. The three secondary metabolites (glycyrrhizin, liquiritin, and isoliquiritin) tested in this study had mutually reinforcing relationships because of their similar metabolic pathways (Xie et al. 2018).

Distribution of different root-associated bacterial communities

The phylum-level composition of the bacterial communities associated with the licorice roots was similar to that described in previous studies, regardless of plant species (Chen et al. 2016; Peiffer et al. 2013). Interestingly, the phylum Oxyphotobacteria accounted for 82.6% of the root endosphere communities in the present study. Oxyphotobacteria is the product of the evolution of oxygenic photosynthesis in the ancestors of Cyanobacteria and it is the most abundant phylum in oceanic picoplankton (Bibby et al. 2001). This result has expanded the distribution range of this phylum and needs to be further clarified. Additionally, the fact that some Cyanobacteria have photosynthetic nitrogen-fixing function is interesting and inner diazotrophic bacteria have been well documented for rice (Zhang et al. 2019), maize (Roesch et al. 2008), sugarcane (Thaweenut et al. 2011) and some gramineous energy plants (Kirchhof et al. 1997), but the function of Oxyphotobacteria in licorice root should be further identified. However, the root nodules on the licorice in different age groups were not enough to investigate the rhizobia communities, resulting from the effects of licorice growing environment and artificial management. Thus, the root nodules and rhizobia were not included in our study. The symbionts were not recovered in the root endophytic population, which might be as consequence of our data processing and rare root nodules during our sampling. However, we would pay further attention to the symbionts in root endophytic communities according to the recent study (Zhang et al. 2019) in the future.

The alpha-diversity of the root-associated bacterial communities significantly decreased with increasing root proximity. This result was consistent with other study that has explained ecological niche selection by plants (Edwards et al. 2015). The Shannon index and OTU richness both increased in the bulk and rhizosphere soils with the increasing age of licorice plants. This trend may be related to the increased soil primary nutrients. In the root endosphere, the Shannon index was stable and the OTU richness displayed little change among the four licorice sites. This confirmed the stability of root endosphere communities.

The composition of the root-associated bacterial communities was well separated according to the three soil–root compartments. Additionally, the rhizosphere soil communities were well divided into four clusters compared with those of the bulk soil due to a greater effect of age. This is because licorice root exudates have greater effects on rhizosphere bacteria than bulk soil bacteria; the higher input of rhizodeposits into the rhizosphere may result in substantial changes in the bacterial community (Shakya et al. 2013; Wagner et al. 2016). Similarly, the composition of the root endosphere communities was stable and no significant differences were observed in their alpha-diversity between the licorice sites. A plausible reason is that these communities, which inhabited the roots of the same licorice species, were less affected by bulk soil characteristics and mainly determined by plant genotype. Similar results have also been obtained in the model plant A. thaliana (Bulgarelli et al. 2012; Hardoim et al. 2015). The genotypic variation amongst licorice plants needs to be further studied to confirm our inference.

Enrichment process and time–decay relationships of core-enriched taxa

Dominant microbial communities play important roles in nutrient cycling in agricultural soils and in microbial colonization of plants (Jiao et al. 2019; Lundberg et al. 2012). Accordingly, the dominant bacterial communities associated with the roots of licorice were identified representatively in this study. Licorice root compartments selectively enriched bacterial communities from the bulk soil and harbored different numbers of enriched taxa across various ages. This enrichment process could be due to the development stages of licorice plants, which was supported by a previous study (Zhang et al. 2018). With increasing licorice age, the rhizosphere soil enriched more taxa, whereas the root endosphere enriched fewer taxa. The amount and chemical composition of root exudates have been found to change during distinct rice growth stages (Zhang et al. 2018). Thus, it is inferred that the licorice plants selected specific communities over time, enriched a variety of taxa in the rhizosphere soil to prevent environmental interference or pathogenic infection. On the contrary, the root endosphere communities were relatively stable. It is possible that more taxa were enriched at first to adapt to the new environment, whereas some unstable taxa were eliminated during voluntary selection by the plants to recover stable original conditions similar to those in wild licorice.

The core-enriched taxa were then defined. This method has been used to select commonly enriched taxa of the same rice species growing in different farmlands (Zhang et al. 2019). However, only one taxon without annotations was obtained from the root endosphere and the root endosphere communities were always stable as mentioned above. Thus, the conserved core-enriched taxa in the bulk and rhizosphere soils were worthy of attention because of precise selections for them. The functional pathways of these core-enriched taxa also had habitat preferences. For example, rhizosphere soil is a relatively anaerobic microenvironment suitable for the metabolism of denitrifying bacteria. Additionally, rhizosphere communities are more affected by plant root exudates that can recruit pathogenic microorganisms (Dombrowski et al. 2017; Huang et al. 2018). These habitat characteristics corresponded to the functions of the taxa present. These functions include chemoheterotrophy, which was found to be the basic lifestyle of most of the bacterial communities in the bulk and rhizosphere soils of the licorice in this study.

Time–decay relationships in microbial communities are consistent among similar environments (Shade et al. 2013). Herein, it was found that the significant turnover slopes of root-associated bacterial communities were much smaller (−0.028 to −0.037) compared with previous studies, which reported the turnover slopes to be between 0 and − 0.3 in air, soil, and freshwater habitats (Jiao et al. 2017). A possible reason for this wide difference is that licorice selected relatively conserved root-associated bacterial communities. There were no decay patterns in root endosphere communities, which was consistent with their stability. The whole and dominant bacterial communities had steeper turnover slopes in the rhizosphere soil, which indicates a faster rate of community temporal turnover (Nekola and White 1999). In addition, this result indicated that the rhizosphere soil bacterial communities were more vulnerable and underwent rapid elimination and replacement of species due to interactions with different root exudates and soil characteristics (Huang et al. 2018; Koberl et al. 2013). However, the turnover slopes of the core-enriched communities both decreased, because these communities were relatively conserved during licorice development.

The links between secondary metabolites and individual taxa mediated by soil potassium

The core-enriched bacterial communities in the three soil–root compartments were affected by different latent factors resulting from niche differentiation. Many studies have reported that the dynamics of the rhizomicrobiome are affected by different soil characteristics, such as soil pH, nitrogen, and water conditions, to varying degrees (Edwards et al. 2015; Xiao et al. 2017; Xie et al. 2018). Similar relationships are elucidated in another study, indicating that soil characteristics can influence the metabolite content in herbs through different assembly of bacterial communities (Huang et al. 2018).

In this study, soil TK was defined as a key factor driving the temporal dynamics of core-enriched bacterial communities in the bulk and rhizosphere soils. Some studies have reported that soil potassium is significantly correlated with rhizosphere bacterial communities in wild oat (Nuccio et al. 2016) and with the relative abundance of Proteobacteria in Japanese barberry (Coats et al. 2014). Some individual taxa were then selected from the core-enriched taxa based on standards of significant correlation with TK and temporal variation in these taxa for better targeting. These individual taxa were sensitive to temporal variation in soil TK and were significantly correlated with the concentrations of different secondary metabolites. This finding demonstrated that soil TK mediates the links between root secondary metabolites and the individual taxa in this study.

Many studies have reported that Acidobacteria is a phylum of oligotrophic bacteria (Chen et al. 2016; Xu et al. 2018). Therefore, it was likely that these taxa that were negatively correlated with TK had a higher relative abundance due to the lower TK content in the wild site. The positively correlated taxa displayed lower relative abundances in the wild site, which was similar to the result of the rhizosphere soil. These taxa may regulate the accumulation of root secondary metabolites in licorice at lower TK levels as indicated by the distinct correlations between their relative abundances and secondary metabolites. Moreover, the class Ignavibacteria, which belongs to the phylum Bacteroidetes, showed opposing temporal variation to TK. Bacteria in this class were negatively related to secondary metabolites in the bulk soil. This is because soil TK is mainly consisted of inorganic potassium, which cannot be used by Ignavibacteria. Ignavibacteria prefers organic matter, including secondary metabolites, available in the bulk soil. This is supported by a study reporting that Ignavibacteria are chemoorganotrophic bacteria (Podosokorskaya et al. 2013).

Proteobacteria and Actinobacteria species are always copiotrophic in bulk soil (Edwards et al. 2015). However, the genus Inquilinus of the phylum Proteobacteria preferred lower TK in the rhizosphere soil in this study. Inquilinus has been isolated from ginseng fields and is able to promote the growth of ginseng (Chun Hwi et al. 2012; Hae-Min et al. 2011). This may also be true for licorice as Inquilinus was positively correlated with the secondary metabolites in the licorice roots. Some studies have found that Nocardioides and Promicromonospora of the Actinobacteria phylum have a growth-promoting effect on ginseng and pine tree roots, respectively (Chun Hwi et al. 2012; Kaewkla and Franco 2017). In addition, the inoculation with the plant growth-promoting rhizobacteria improve the tolerance of licorice to salt stress (Egamberdieva et al., 2016). These taxa may therefore have the potential ability to stimulate the licorice resistance to environmental disturbance and the accumulation of root secondary metabolites in licorice, given their positive correlation with glycyrrhizin and isoliquiritin.

In the rhizosphere soil, the family Xanthomonadaceae and the genus Variovorax displayed opposing temporal variation to soil TK and were negatively correlated with secondary metabolites in the licorice roots. It has been reported that Xanthomonadaceae is enriched during the consecutive monoculture of Rehmannia glutinosa, which is detrimental to plant growth (Wu et al. 2018). Additionally, Variovorax is one of the bacterial genera that displays increased abundance in response to apple replant disease (Lucas et al. 2018). These two taxa may also be harmful to licorice growth and inhibit the production of secondary metabolites despite of the Variovorax paradoxus is commonly described as plant growth promoting rhizobacteria (Kudoyarova et al., 2019). That may be due to the different environmental conditions and plant species. Collectively, these individual taxa exhibited habitat-specific features and were associated with the adjustment of secondary metabolite accumulation in the licorice roots under regulation by TK. The interactions of these taxa are very complex and may not be isolated from each other. Further research is needed to investigate the vital role of TK in the regulation of root secondary metabolites in licorice under controlled conditions.

Conclusions

This study investigated the temporal succession of root-associated bacterial communities and simultaneous variation in soil characteristics and root secondary metabolites in licorice plants of different ages. The core-enriched communities differed among the root-associated bacterial communities, and their time–decay relationships also varied. Soil TK was defined as the key factor regulating individual core-enriched taxa. These individual taxa with distinct habitat adaptation to soil TK played different roles in adjusting the accumulation of root secondary metabolites in licorice. Overall, the results provide valuable information for better understanding the soil potassium-regulated temporal dynamics of root-associated core bacterial communities in licorice. Not only do the results shed light on the core community level, but also on the level of individual taxa. Future research should focus on the isolation of these individual taxa and on confirming their interactions with soil potassium and root secondary metabolites in licorice in controlled experiments. This research should be applied to manage the fertilization and soil bacterial communities efficiently in practical licorice cultivation.