1.1 The effect of the growing year on the content of bioactive compounds in the roots of medicinal licorices
The results of two-way ANOVA showed that the contents of the bioactive compounds (glycyrrhizic acid (GIA), liquiritin (LI) and total flavonoid (GTF)) were not significantly affected by the interaction effect between growing year (1–3) and plant species (Glycyrrhiza uralensis, Glycyrrhiza inflata, and Glycyrrhiza glabra) (P > 0.05) (Table 1). However, the contents of the bioactive compounds were significantly affected by the main effect growing year (1–3) (P < 0.05) (Table 1 and Figure 1), with a trend of stable increase in the contents observed with each growing year. As shown in Figure 1, the contents of GIA, LI and GTF in the year 3 were significantly higher than those in year 1, the contents of GIA and GTF in year 2 were significantly higher than those in year 1, and the contents of GIA and LI in year 3 were significantly higher than those in year 2.
Details of the temperature, rainfall, leaves and soil factors during the experimental period are presented in Supplementary Table S1. In addition, Spearman correlation analysis showed that the content of bioactive compounds was significantly correlated with soil physicochemical properties and nutritional components of leaves (Table 2). In terms of soil physicochemical factors, GIA and LI had a very significant positive correlation with soil ammonium nitrogen (SAN) (R2 > 0; P < 0.01), but had a very significant negative correlation with total salt (TS) (R2 < 0; P < 0.01); GTF had a very significant positive correlation with SAN (R2 > 0; P < 0.01), but had a very significant negative correlation with nitrate nitrogen (SNN) (R2 < 0; P < 0.01), and had a significant negative correlation with TS (R2 < 0.01; P<0.05). In terms of leaf nutrition (Table 2), GIA, LI and GTF were very negatively correlated with water content (PWC), total potassium (PTK) (Figure 2 a, b, c) and crude fiber (CF). In addition, GTF also showed a very significant negative correlation with total phosphorus (PTP) (Figure 2d). There was a significant relationship between leaf nutrients and a small number of variables. Specifically, PWC was positively correlated with SNN, but negatively correlated with SAN; organic carbon (POC) was positively correlated with soil total potassium (STK); total nitrogen (PTN) was positively correlated with POC; total phosphorus (PTP) was very positively correlated with SNN, TS and PTN; PTK was very positively correlated with SNN, TS, PWC, PTN and PTP, but negatively correlated with SAN; CF was very positively correlated with PWC and PTK.
1.2 Composition of bacterial community in the root of medicinal liquorices
In 27 samples, a total of 1,979,531 effective sequences were obtained after quality control. The sequencing results for each sample are listed in Supplementary Table S2. The effective sequences were clustered into operational taxonomic units (OTUs) with 97% identity, and a total of 2,432 OTUs were obtained. The rarefaction curves were normalized to the lowest number of sequences, showed that the number of OTU in each sample increased gradually with quantity of sequence, thus confirming that the amount of sequencing data is adequate (Figure 3). According to the OTUs, 28 phyla, 45 classes, 106 orders, 216 families, 508 genera and 313 species were annotated. Figure 4a shows the 10 bacterial phyla with the greatest abundance in the bacterial community in the roots of medicinal licorices. Proteobacteria dominated the observed sequences at the phylum level, representing 79.917%, 64.420%, 80.265%, 82.848%, 65.733%, 87.886%, 64.544%, 75.024% and 85.847% of the total number of species in E.W, R.W, S.W, E.D, R.D, S.D, E.G, R.G and S.G, respectively. In addition, Actinobacteria occupied a large part of the relative abundance in E.G (18.960%) and R.D (22.171%), respectively. In addition, the abundance of Bacteroidetes was high in the R.W, accounting for 19.742%; the abundance of Firmicutes was also high in the R.G, E.G and R.W samples, accounting for 8.735%, 7.466% and 6.287%, respectively (Figure 4a).
In addition, t-test analysis of the two groups showed significant differences in the relative abundance of Tenericutes in samples E, G and W (p < 0.05) (Figure 4b). As shown in Figure 4b, the relative abundance of Tenericutes in E.G was significantly higher than that in E.D, E.W and S.G (p < 0.05). The relative abundance of Tenericutes in E.D was significantly higher than that in E.W (p < 0.05). Moreover, the relative abundance of Tenericutes in S.W was significantly higher than that in E.W (p < 0.05). These findings indicated significant differences in the relative abundance of Tenericutes among plant species and growing year (p < 0.05).
In terms of genus, unidentified-Rhizobiaceae was more abundant than other genera (Figure 5), with the abundance in individual samples ranging from 3.398% (E.D) to 33.985% (S.W). The abundance of Pseudomonas was high in the E.W, E.D, R.W, S.D and R.G samples, accounting for 15.195%, 7.165%, 6.326%, 16.584% and 7.772%, respectively. Pantoea (15.386%) and Halomonas (7.630%) were found to be the most dominant in E.D sample.
We used linear discriminant analysis (LDA) effect size (LEfSe) to identify discriminative taxa among different species of medicinal licorices and different growing years. As shown in Figure 6, the results of LEfSe analysis of all years and species based on the rank sum test revealed a total of 16 biomarkers with significant differences that were contained in the E.D (5 taxa), E.W (7 taxa), R.G (2 taxa) and R.W (2 taxa) groups. No discriminative taxa were observed in year 3. These biomarkers included Variovorax, Nocardiopsis, Methylophaga, Pelagibacterium, Halomonas and Sinomicrobium at the genus taxonomic level (Figure 6b).
Details of the composition of the top 10 dominant bacteria at other classification levels are listed in Supplementary Table S3. Specifically, Alphaproteobacteria, Gammaproteobacteria, unidentified-Actinobacteria dominate at the class taxonomic level; the dominant species at the order taxonomic level are Rhizobiales, Sphingomonadales and Gammaproteobacteria; the dominant species at the family taxonomic level are Rhizobiaceae, Pseudomonadaceae and Sphingomonadaceae; the dominant species at the species taxonomic level are Pantoea-brenneri, Neorhizobium-huautlense, and Pseudomonas-psychrotolerans.
1.3 Effects of year and species on alpha diversity and beta diversity in the bacterial community of roots.
The alpha diversity index of each group is shown in Supplementary Table S4. The results of two-way ANOVA showed that the alpha diversity indexes (Shannon, Simpson, Chao1 and ACE) were not significantly affected by the interaction between plant species and growing year and were not significantly affected by the main effect species and main effect growing year (Table 3).
However, beta diversity analysis showed significant differences in the endophytic bacterial community among the different groups. As shown in Figure 7, Wilcox rank sum test based on the UniFrac distances showed significant differences in beta diversity between E.W and S.W (P < 0.05), E.W and R.W (P < 0.05), E.W and E.G (P < 0.05), R.W and R.D (P < 0.05), and E.D and E.G (P < 0.01), which indicated the existence of significant differences in endophytic bacterial community of roots of medicinal licorices between different species and different growing years.
1.4 The relationship between the dominant phylum and genus of endophytic bacteria and the bioactive compounds, soil physicochemical properties and leaf nutrition.
Spearman correlation analysis showed that there was a significant relationship between dominant bacteria phylum and bioactive compounds, soil physicochemical properties and leaf nutrition (Table 4). Specifically, Proteobacteria showed a very significant negative correlation with PTN; Actinobacteria showed a significant negative correlation with SAN, GIA and LI (R2 < 0; P < 0.05); Firmicutes showed a significant positive correlation with PTN and CF; and Acidobacteria showed a significant positive correlation with POC (R2 > 0; P < 0.05).
As shown in Figure 8, there was a significant relationship between the dominant bacteria genus and the bioactive compounds of roots, soil physicochemical factors and leaf nutrition. Specifically, unidentified-Rhizobiaceae had a significant positive correlation with SAN, GIA, GTF and LI, but a very significant negative correlation with CF; Pantoea had a significant positive correlation with PTK; Pseudomonas had a significant positive correlation with PTP; Halomonas had a significant positive correlation with SNN, TS, PWC, PTK and CF, but a very significant negative correlation with SAN, GIA, GTF and LI; Nocardiopsis had a significant positive correlation with TS, PWC, PTP, PTK and CF, but a very significant negative correlation with SAN, GIA, GTF and LI; Novosphingobium had a significant negative correlation with PWC, but a very significant positive correlation with GTF; and Sinomicrobium had significant positive correlation with SNN, PWC, PTK and CF, but a significant negative correlation with GIA, GTF and LI. It is worth noting that we found that the biomarkers with significant differences at the genus taxonomic level of the dominant bacterial community (including Nocardiopsis, Halomonas and Sinomicrobium) had a significant negative correlation with GIA, GTF, and LI (R2 < 0; P < 0.05). Details of the interrelationships between the biomarkers with significant differences in the class, orders, families and genera of dominant bacteria and bioactive compounds, soil physicochemistry and leaf nutrients are shown in Supplementary Table S5. These results showed that the biomarkers with significant differences (except Methylophaga) at each classification level were significantly negatively correlated with GIA, GTF, and LI, while all showed a significant positive correlation with PTK.
1.5 Total flavonoids explained the difference in composition and distribution of endophytic bacteria community in the roots of planting licorices to the greatest extent.
Distance-based redundancy analysis (db-RDA) based on the Bray–Curtis distance showed that the bioactive compounds, soil physicochemical and leaf nutrient components had significant effects on the endophytic bacterial community in the roots of medicinal licorices (Figure 9, Table 5). Specifically, the contents of GIA, GTF, and LI all had a significant effect on the endophytic bacterial community in the roots (P < 0.05). Among them, the content of GTF explained the difference in the composition and distribution of endophytic bacterial communities in the roots of cultivated medicinal liquorices to the greatest extent (r2 = 0.638, P < 0.01). Among the soil environment factors, TS of soil was identified as the factor that most significantly affects the endophytic bacterial community, followed by SAN and SNN. Among the leaf nutrients factors, PWC was identified as the factor that most significantly affects the endophytic bacterial composition, followed by PTK and PTP. In addition, Mantel tests revealed that the combination of medicinal components is the environmental factor with the greatest correlation with the endophytic bacterial community, indicating that the combination has the greatest impact on the microbial community (r = 0.260; P = 0.008) (Supplementary table S6).