The endophytes Bacillus cereus LgD2 and Fusarium oxysporum TPB were isolated from P. fargesii, preserved in 50% glycerol, and stored at −80°C (Yan et al., 2022). Both LgD2 and TPB strains were removed and rapidly shaken at 37°C in a water bath for 10 s for activation.

The activated LgD2 bacterial solution was aspirated and inoculated in LB agar medium using the four-zone delineation method and incubated overnight at 37°C. A single colony was inoculated into 2 mL LB liquid medium (Solarbio Life Sciences, China) and incubated at 37°C for 5 h with shaking. Subsequently, 1 mL of this culture was inoculated into a triangular flask containing 500 mL liquid medium and shaken at 200 rpm at 37°C until the OD₆₀₀ was 0.5–0.6. The culture was then centrifuged at 8.93 g for 5 min, and 500 mL of sterile water was added to resuspend the pellet and create a bacterial suspension (Costa Júnior et al., 2020; Liu et al., 2022).

An activated TPB mycelium was selected and inoculated into Potato Dextrose Agar medium (Solarbio Life Sciences, China) and cultured at 28°C for 5–7 days. The surface of the mycelium was rinsed with sterile water, and then filtered through four layers of sterile gauze to obtain a spore suspension (Andrade et al., 2018; Wang et al., 2021). The spore suspension was diluted in sterile water to 1 × 10⁷ cells/mL and kept for subsequent experiments.

On April 1, 2022, two-year-old P. polyphylla seedlings were purchased from Yunnan Suixitang Biotechnology Co. (China). They were planted in sterile plastic pots containing nutrient soil sterilized at 121°C for 1 h. Five plants were placed in each plastic pot and provided with 100 mL sterile water every 2–3 days. The pots were placed in the backyard greenhouse of Yunnan Agricultural University, 50% of the area of which was covered by a shade net, and its top was closed by a transparent plastic film. The strain inoculation treatments were carried out on May 1, 2022, with roots inoculated with bacterial or spore suspension every 15 dpi; five pots were inoculated for each treatment, with 50 mL of bacterial or spore suspension being poured into each pot separately. Sterile water (50 mL) was used as the control treatment. The rhizomes of P. polyphylla inoculated with endophyte strains three and six times (45 dpi and 90 dpi treatments) were collected on June 14, 2022 and July 29, 2022, respectively, cleaned, and immediately stored in liquid nitrogen.

Three replicates of each of 45 dpi and 90 dpi rhizomes of P. polyphylla treated with LgD2, TPB, or sterile water (18 samples in total) were removed from liquid nitrogen and ground (0.1 g for each rhizome sample). RNA was extracted from the rhizomes of P. polyphylla using the HiPure HP Plant RNA Mini Kit (Guangzhou Magen Biotechnology, China) RNA integrity was assessed using 1% agarose gel electrophoresis and an Agilent 5,300 Bioanalyzer (Agilent, USA).

Briefly, 10 μg of total RNA of acceptable integrity, concentration, and purity was selected from each rhizome sample; 18 samples were used to construct libraries using the TruSeq Stranded mRNA LT Sample Prep Kit. Transcriptome sequencing was performed at Wuhan SeqHealth Technology Co., Ltd. (China), using an Illumina HiSeq 4,000 platform. The raw data obtained from MGISEQ-T7 sequencing were converted into sequence data using base calling (FASTQ format) to obtain the most original sequencing data file.

De novo transcriptome assembly and splice variant calling were performed using the Trinity software to obtain transcript sequences, with the longest transcripts as Unigenes. The SortMerna software was used to filter rRNA reads, which were error-corrected and deduplicated using the Rcorrector and fastUniq software. Data quality control was performed using fastp (version 0.23.0). Transcript splicing quality was evaluated using BUSCO (version 5.5.0).

Protein sequences obtained from the Coding sequence (CDS) prediction by Unigene were annotated using the primary mission of Universal Protein Resource (UniPprot), Non-Redundant (NR), Pfam, evolutionary genealogy of genes: Non-supervised Orthologous Groups (eggNOG), Gene Ontology (GO),¹ Kyoto Encyclopedia of Genes and Genomes (KEGG),² and CAZYmes databases. Unigene sequences predicted without CDS were created in the NR and Rfam databases for annotation. Fragments per kilobase per million (FPKM) and read count values were calculated for each Unigene using Bowtie2 and eXpress. The absolute value of log2(FC) > 1 and value of p <0.05 were used as criteria to indicate that the gene was a differentially expressed Unigene. Hierarchical clustering analysis of differentially expressed genes (DEGs) was performed to determine the expression patterns of genes among different groups and samples. GO and KEGG pathway enrichment analyses were performed for DEGs using hypergeometric distribution-based R.

The BLAST software and Pfam databases were used to identify transcripts belonging to CYP450 and UGTs. These transcripts were then translated into protein sequences using TBtools (version 1.120). The CYP450 protein sequences of Arabidopsis thaliana downloaded from TAIR,³ whereas the UGT sequences of A. thaliana, Isatis tinctoria, Malus domestica, Carthamus tinctorius, Avena strigosa, and Arachis hypogaea were downloaded from the NCBI database. A phylogenetic tree of CYP450 and UGT gene members was constructed using the maximum likelihood method (Chen et al., 2020). Phylogenetic topology was assessed using 1,000 bootstrap replicates.

Gene co-expression networks were constructed using the weighted gene co-expression network (WGCNA) method in R software (version 3.2.2) (Langfelder and Horvath, 2008; Li et al., 2022). Genes with FPKM <1 in 90% of samples were filtered out, and the expression matrices of the top 8,000 genes were selected as input files in WGCNA for the identification of gene modules with strong co-expression. By calculating the correlation coefficients and clusters of the levels of expression of each sample, samples with low correlation or those that could not be clustered on the tree diagram were removed. The WGCNA network was constructed and modules were detected based on the phenotypic traits of polyphyllin. Candidate CYP450 and UGT genes were selected from the important templates most relevant to the phenotypic data, based on edge and KME values assigned to the assumed hub genes of the polyphyllin biosynthesis network. Final visualization was performed using Cytoscape (version 3.10.0) (Doncheva et al., 2019).

Total RNA of acceptable purity and integrity was used for the synthesis of first-strand cDNA using the All-in-One RT 5× MasterMix (ABM, Canada) alongside the AccuRT Genomic DNA Removal Kit. The EvaGreen 2X qPCR MasterMix kit (ABM, Canada) was used to validate the levels of expression of five genes related to the steroid saponin biosynthesis pathway. Primers were designed using Primer3Plus⁴ (Supplementary Table S1). The reaction system (20 μL) consisted of 10 μL of 2× SYBR Premix Ex Taq™, 2 μL cDNA template, 0.5 μL of upstream and downstream primers, and 7 μL of ddH₂O. PCR conditions were as follows: initial denaturation at 95°C for 2 min, followed by 40 cycles of denaturation at 95°C for 10 s, annealing at 60°C for 15 s, and extension at 72°C for 15 s. Four biological replicates were used for each treatment. The relative gene expression was analyzed using the 2^-ΔΔCq method and β-actin was used as the internal reference gene (Bustin et al., 2005; Nolan et al., 2006).

According to the standard for the determination of the medicinal content of P. polyphylla in “Chinese Pharmacopoeia 2020 Edition,” high-performance liquid chromatography (General 0512) was used to determine the content of PPI, PPII, PPVII, PPH, and PPD in P. polyphylla. Briefly, rhizome samples of P. polyphylla treated with TPB and LgD2 for 90 dpi were dried at 45°C to a constant weight, and then ground into powder using a Chinese herbal medicine pulverizer. Next, 0.5 g of powder (passed through sieve No. 3) was obtained, precisely weighed, and placed in a 50 mL volumetric flask. Then, 20 mL methanol was added to the flask, immersed for 30 min, ultrasonicated for 30 min for extraction, and cooled at 25°C. Additional methanol was used to make up for the loss of weight, and the mixture was shaken well, centrifuged at 8.93 g for 5 min, and then filtered through a 0.45 μm micropore filter membrane. The content of PPI, PPII, PPVII, PPH, and PPD in the solution was determined by high performance-liquid chromatography (HPLC) using a C18 column (Zorbax ODS 4.6 × 250 mm, 5 μm; Agilent Technologies, USA) at a detection wavelength of 203 nm, column temperature of 30°C, and flow rate of 1.0 mL/min. PPI, PPII, PPVII, PPH, and PPD controls were purchased from Sichuan Vicki Biotechnology Co. Ltd. (China). Appropriate amounts of PPI, PPII, PPVII, PPH, and PPD controls were obtained, precisely weighed, and mixed with methanol to obtain solutions containing 0.4 mg of each per 1 mL of methanol. For analysis, 20 μL of each of the control and test solutions were precisely aspirated, using octadecylsilane-bonded silica gel as filler and acetonitrile (solvent A) and water (solvent B) as the moving phase in a linear gradient: 0 min: 30% A; 40 min: 60% A; 41 min: 100% A; 44 min: 80% A (Yang et al., 2013; Yan et al., 2022).

IBM SPSS Statistics (version 25) was used to calculate the mean and standard deviation. An independent samples t-test was used to statistically evaluate the differences between the two treatments. *Asterisks indicate significant differences from the control group (t-test; ∗ p < 0.05; ∗∗ p < 0.01). The “heatmap” package in R (version 4.2.2) was used for clustering heatmap plots, whereas the “ggplot2” package was used for visualization.

To investigate the potential effect of endophytes on polyphyllin biosynthesis in P. polyphylla, the content of polyphyllin in P. polyphylla was determined before and after treatment with TPB and LgD2. Inoculation with TPB and LgD2 led to the increased accumulation of polyphyllin in P. polyphylla (Supplementary Table S2). The contents of PPVII, PPD, PPH, and PPI were significantly increased by 42.77, 26.95, 21.96, and 15.27%, respectively, in the TPB treatment group. Likewise, the contents of PPVII, PPD, and PPH were significantly increased by 24, 69.38, and 19.06%, respectively, whereas the PPII content was significantly decreased by 45.97% in the LgD2 treatment group. These results indicated that inoculation with TPB or LgD2 significantly induced the biosynthesis and accumulation of polyphyllin in P. polyphylla.

To identify the intrinsic mechanism by which the TPB and LgD2 strains promote polyphyllin biosynthesis, total RNA isolated from the rhizomes of P. polyphylla was sequenced under different treatments using the Illumina HiSeq 4,000 platform. Eighteen sequencing libraries constructed from three biological replicates each of 45 and 90 dpi TPB-, LgD2-treated, and control rhizome samples yielded a total of 861,041,222 clean data reads. The Q30 of each sequencing data set was >95%, and 788,438,886 dedup reads were retained after UID deduplication. Accordingly, 1,329,726 contigs were obtained by assembling the dedup reads using Trinity, with an N50 of 638 bp, average length of 537.49 bp, longest read length of 16,964 bp, and average GC content of 45.11%. BUSCO analysis suggested 81% coverage and completeness for transcriptome assembly (Table 1). Then, 1,262,695 (94.96%) of these Unigenes were annotated using the KEGG, UniProt, NR, Pfam, Rfam, eggnog, and GO databases (Supplementary Table S3). The transcriptome sequence data were stored in the NCBI Short Read Archive (SRA) under the BioProject registration number PRJNA935848.

DEGs were screened from plants inoculated with TPB and LgD2 strains and functional annotation was performed. Gene expression was significantly different between plants under treatment with different strains and inoculation times (Supplementary Figure S1). GO annotation of DEGs revealed that TPB treatment significantly upregulated phosphorus metabolism, protein modification, oxidation–reduction, and phosphotransferase, methyltransferase, acyltransferase, and metal ion transmembrane transporter activities, whereas nitrogen metabolism, RNA metabolism, metal ion binding, and cation binding were significantly downregulated. Additionally, LgD2 treatment upregulated oxidation–reduction processes, responses to abiotic stimuli, cellular redox homeostasis, carbohydrate and ATP biosynthesis, phosphotransferase, and UDP-glycosyltransferase activities, whereas carbohydrate metabolism, purine nucleotide metabolism, and metal and oxygen binding were downregulated (Supplementary Figure S2). KEGG pathway analysis showed that TPB treatment upregulated DEGs mainly enriched in the pathways of carbon fixation, glycolysis/gluconeogenesis, and plant-pathogen interaction in photosynthetic organisms. Conversely, downregulated DEGs were mainly enriched in the pathways of bacterial invasion of epithelial cells, pathogenic E. coli and Salmonella infections, and shigellosis. LgD2 treatment also upregulated carbon metabolism, the tricarboxylic acid (TCA) cycle, plant interactions with pathogens, Vibrio cholerae infection, nitrogen metabolism, oxidative phosphorylation, and other pathways (Figure 1). These analyses revealed that TPB and LgD2 treatment promoted carbohydrate accumulation and redox homeostasis, increased the levels of phosphorylation, glycosylation, and acylation, and reduced the binding capacity to metal ions and cations. In particular, plant-pathogen interactions were significantly upregulated, whereas the ability of some pathogenic bacteria for infection and bacterial invasion of epithelial cells were decreased, suggesting that endophytic inoculation may promote the defense response of P. polyphylla against some bacteria and pathogens.

The terpene backbone and steroid biosynthesis pathways are located upstream of polyphyllin biosynthesis. According to the annotation using the KEGG database, 722 Unigenes were identified as involved in “steroid biosynthesis” (ko00100), encoding 31 enzymes, and while 823 Unigenes were identified as related to “terpenoid skeleton biosynthesis” (ko00900), encoding 36 enzymes. Twenty Unigenes were highly homologous to the functional signature and DEGs under different strain treatments. More specifically, eight DEGs were present in the terpene skeleton biosynthesis pathway, whereas 12 DEGs were present in the steroid biosynthesis pathway (Figure 2).

Unexpectedly, the remaining genes, except farnesyltransferase type-1 subunit alpha (FNTA), 1-deoxy-D-xylulose-5-phosphate synthase (dxs), and TAG lipase acyltransferase (TGL4), were downregulated under strain treatment compared with those in the control. An increased number of genes and differential genes were annotated to the MVA pathway, corroborating that the MVA pathway is the main pathway for polyphyllin biosynthesis (Liu et al., 2016; Upadhyay et al., 2018). These genes were downregulated at 90 dpi, except for hydroxymethylglutaryl-CoA reductase (HMGCR), geranylgeranyl diphosphate synthase (GGPS), plant 4,4-dimethylsterol C-4alpha-methyl-monooxygenase (SMO1), and plant 4alpha-monomethylsterol monooxygenase (SMO2). Previous studies have shown that PPI, PPII, and PPD are diosgenins, whereas PPVII and PPPH are pennogenins, all of which are synthesized from cholesterol (Deng et al., 1999; Li et al., 2001). Therefore, cholesterol accumulation is key to the biosynthesis and accumulation of polyphyllin. Strain treatment did not significantly promote cholesterol biosynthesis; it had the opposite effect, showing a decreasing tendency in the levels of cholesterol over time.

To further validate the transcriptomic results, the expression of five key genes involved in polyphyllin biosynthesis was investigated. The qRT-PCR analysis revealed that the expression of four genes involved in cholesterol biosynthesis, 3-hydroxy-3-methylglutaryl coenzyme A reductase (HMGR), squalene synthase (SQS), cycloartenol synthase (CAS), and oxidosqualene cyclase (SE), was downregulated in P. polyphylla under TPB and LgD2 treatment. The only exceptions were the CAS gene at 45 dpi rhizomes under TPB treatment and the SE gene at 90 dpi rhizomes treated with LgD2 (Supplementary Figure S3). Overall, the genes upstream of polyphyllin biosynthesis mainly showed a trend of downregulated expression under different strain inoculations and time treatments, which was basically consistent with our RNA sequencing analysis, indicating reliable the RNA-seq results. CYP90B27 is responsible for C-22 hydroxylation in the steroidal alkaloid biosynthesis pathway, a downstream pathway in polyphyllin biosynthesis (Yin et al., 2018). The expression of the CYP90B27 gene was upregulated in P. polyphylla treated with either strain at 90 dpi, which may be important for the accumulation of polyphyllin.

In this study, Pfam annotation and the BLAST algorithm were used to identify 257 CYP450s and 71 UGTs. Forty-five CYP450s were downloaded from the Arabidopsis thaliana database and used to construct a phylogenetic tree; in total, 302 CYP450s were characterized into 13 subfamilies (Figure 3A). Through phylogenetic analysis, the 257 CYP450s of P. polyphylla were assigned to the CYP71, CYP94, CYP81, CYP704, CYP90, CYP85, CYP86, and CYP55 families. Among them, CYP71 and CYP94 had the largest number of CYPs, followed by CYP86, CYP90, and CYP55, whereas CYP704 had the lowest number of CYP450s. Additionally, 132 UGTs of A. thaliana, Isatis tinctoria, Malus domestica, Carthamus tinctorius, Avena strigosa, and Arachis hypogaea were downloaded and used to construct a phylogenetic tree. In total, 203 UGTs were classified into 11 subfamilies (Figure 3B). The 71 UGTs of P. polyphylla belonged to the UGT74, UGT75, UGT85, UGT87, UGT91, UGT80, UGT71, and UGT73 families. Among them, UGT87 and UGT80 had the highest number of UGT genes, followed by UGT91, UGT74, and UGT75, whereas UGT85 had the lowest number of genes.

In total, 257 CYP450s and 71 UGTs genes were screened for differential expression and 32 CYP450s and 11 UGTs were differentially expressed between the control and TPB and LgD2 treatments. These differentially expressed CYP450s were grouped into three categories, with the most common being CYP71, followed by CYP86 and CYP94 (Figure 4A). Interestingly, genes in the first category were highly expressed in the control at 45 dpi, those in the second category were highly expressed in the control at 90 dpi, whereas those in the third category were upregulated at both 45 and 90 dpi under TPB and LgD2 treatments. Similarly, UTGs were grouped into two categories, with ten genes belonging to UGT87, whereas one was assigned to UGT80 (Figure 4B). The first class of UGTs was upregulated at 90 dpi, whereas the second class was highly expressed at 45 dpi in plants under TPB treatment. Differentially expressed CYP450s were both up- and down-regulated, whereas differentially expressed UGTs were mostly upregulated at different times and under different treatments.

Annotation of these CYP450 and UGT genes to the Uniprot, NR, eggNog, GO, KEGG, and CAZYmes databases (Supplementary Table S4) revealed that the functions of these differentially expressed UGTs were mainly related to the formation of glycosidic bonds (GTs), with GT1 accounting for the major part. GO annotation showed that among them, TRINITY_DN2024_c0_g1 and TRINITY_DN2707_c1_g1 belonged to the metabolic pathway of secondary metabolite biosynthesis, functioning in the transfer of hexosyl groups. Moreover, eggNOG database annotations revealed that TRINITY_DN2707_c1_g1 has a flavonol 3-O-glucosyltransferase activity. TRINITY_DN2024_c0_g1 was annotated as scopoletin glucosyltransferase-like protein in the NR database. Annotation of TRINITY_DN541052_c0_g1 and TRINITY_DN117284_c0_g1, which are differentially expressed CYP450 genes, suggested that they function as GT31 and GT1, respectively, in the CAZYmes database. NR annotations identified TRINITY_DN541052_c0_g1 as a lanosterol 14-alpha-demethylase and TRINITY_DN117284_c0_g1 as an isoflavone 2′-hydroxylase-like enzyme.

The coexpression network was constructed by inputting transcript data from P. polyphylla under different treatments and choosing a soft threshold, β = 14, at which R² = 0.9, indicating the successful generation of the model (Figure 5A). Unigenes with similar expression patterns were clustered in the same branch, with different colors indicating different coexpression modules in the branch (Figure 5B). A threshold of 30 was chosen for each module, and the 8,000 Unigenes were divided with the highest correlation into 26 modules (Figure 5C). Then, the genes in the steroid and terpene backbone biosynthesis pathway were identified by selecting the expression modules that were significantly associated with polyphyllin content, and the correlation coefficients between the biosynthesis of polyphyllin and gene modules were calculated (Figure 6A). The closer the correlation coefficient was to one, the more similar the expression or distribution patterns of the metabolites and genes in the modules of P. polyphylla inoculated with different strains. Among them, PPVII and PPH were most correlated with the pink module, with correlation values of 0.75 and 0.76, respectively; PPD was most correlated with the yellow module (0.68); PPII was most correlated with the brown module (0.44); and PPI was most correlated with the light green module (0.51), indicating that the genes in the pink, yellow, brown, and light green modules were highly correlated and might be involved in polyphyllin biosynthesis. In total, 1880 genes were co-clustered with the expression modules in the four modules.

Of the 1880 genes, were selected and screened those involved in the terpene skeleton and steroid biosynthesis pathways to identify genes shared in the CYP450 and UGT phylogenetic trees. Using the four prediction modules, six genes in the terpene skeleton biosynthesis pathway, six genes in the steroid biosynthesis pathway, 17 CYP450s, and 2 UGTs were identified. The clustering heatmap shows the expression patterns of these genes under different treatments (Supplementary Figure S4). The expression of these genes was upregulated under TPB and LgD2 treatment, especially that of CYP90B1, which may play a key role in promoting polyphyllin biosynthesis.

To explore the function of these CYP450 and UGT genes, these candidate genes were annotated into the Uniprot, NR, eggNog, GO, KEGG, and CAZYmes databases (Supplementary Table S5). Annotation in the CAZYmes database revealed that TRINITY_DN2002_c0_g1 and TRINITY_DN19349_c0_g1 were both GT1, whereas NR database annotation resulted in their identification as UDP-glycosyltransferase 92A1-like and UDP-glycosyltransferase 89B1, respectively. TRINITY_DN7787_c1_g1 was annotated to the diterpenoid biosynthesis pathway in the KEGG database and to GT31 in the CAZYmes database, and Uniprot and NR annotations revealed that this gene encodes an Ent-kaurenoic acid oxidase.

To further investigate the roles of CYP450s and UGTs in polyphyllin-related modules, the top genes with the highest weights in the pink and yellow modules were selected for network correlation analysis. CYP450, KAO, CYP734A1, and ispE play key roles in the pink module, whereas CYP450, acetyl-CoA C-acetyltransferase (atoB), 4-diphosphocytidyl-2-C-methyl-D-erythritol kinase (ispE), CYP90B1, and mevalonate kinase (MVK) play key roles in the yellow module (Figure 6B).