Genome-Wide Identification of Genes Related to Biosynthesis of Phenolic Acid Derivatives in Bletilla striata at Different Suspension Culture Stages

To screen the genes regulating the biosynthesis of phenolic acid derivatives from the genome of Bletilla striata, we designed a suspension culture system to sample the cells for the following experiments. The contents of four phenolic acid derivatives were determined by high-performance liquid chromatography, and several full-length transcriptome sequencings of RNA samples at 10 time points were performed for bioinformatics analysis. The correlation analysis was used to identify and verify the key DEGs involved in the biosynthesis of the four phenolic acid derivatives. The results showed that the contents of p-hydroxybenzylalcohol (HBA), Dactylorhin A, Militarine, and Coelonin peaked at 33 days postinoculation (Dpi), 18 Dpi, 39 Dpi, and 39 Dpi of the culture system, respectively. Based on transcriptome data, 80 DEGs involved in the biosynthesis of phenolic acid derivatives were obtained. The KEGG pathway enrichment analysis classified them mostly into five metabolic pathways: phenylpropane biosynthesis, starch and sucrose metabolic, cyanoamino acid metabolism, gluconeogenesis and glycolysis, and phenylalanine metabolism. qPCR analysis revealed that the relative gene expression levels were consistent with the overall trend of transcriptome sequencing results. Among them, 14, 18, 23, and 41 unigenes were found to be involved in the synthesis of HBA, Dactylorhin A, Coelonin, and Militarine, respectively. These unigenes laid a solid foundation for elucidating the biosynthesis mechanism of phenolic acid derivatives in suspension cells of B. striata.


INTRODUCTION
Bletilla striata is a perennial herb of Orchidaceae andalso a traditional and precious Chinese herbal medicine with a high concentration of medicinal metabolites. For thousands of years, it has been prescribed in traditional Chinese medicine to treat hematemesis and promote wound healing (Zhao et al., 2021). A recent report detailed that B. striata has therapeutic effects on gastrointestinal bleeding, hemoptysis, traumatic bleeding, and postpartum hemorrhage (Zhang et al., 2022). Phenolic acid derivatives are a class of compounds with a polyhydroxyphenol structure, which are important secondary metabolites for supporting most pharmacological activities, such as antioxidant, anti-tumor, and anti-inflammatory, of B. striata (Kassim et al., 2010;Anantharaju et al., 2016;Nigro et al., 2017;Milovanović et al., 2021;Yuan et al., 2021). Among the metabolites, p-hydroxybenzylalcohol (HBA), Dactylorhin A, Millitarine, and Coelonin were the most abundant active components in B. striata. In our previous study, we found that the content of these secondary metabolites could reach 0.793 mg/g, 7.792 mg/g, 9.447 mg/g, and 0.345 mg/g of dried suspension cells, respectively (Pan, 2019). A number of studies further demonstrated that HBA has anti-inflammatory, antioxidation, anti-convulsion, anti-depression, and ameliorating effects on memory impairment, as well as protective benefits against brain damage caused by cerebral ischemia and sedative and hypnotic effects (Luo et al., 2017;Ding et al., 2019;Zhang D. et al., 2020). In recent years, many kinds of research have further confirmed the anti-cancer and antioxidant properties of dihydrophenanthrene compounds (Boudjada et al., 2019;Jiang et al., 2019). In conclusion, HBA, Dactylorhin A, Millitarine, and Coelonin are natural chemicals with the potential for drug development.
However, a variety of phenolic acid derivatives have been identified in species such as B. striata, Flickingeria fimbriata, (Wu et al., 2017) and Dendrobium scabrilingue (Sarakulwattana et al., 2020), but their biosynthetic mechanisms and the synergistic regulation mechanisms of multiple genes are still unclear. Therefore, it is important to find the genes that regulate the biosynthesis of secondary metabolites through genetic improvement way or engineering bacteria.
Former studies found that the biosynthesis pathways of phenolic acids mainly include the shikimic acid metabolic pathway, the phenylpropanoid metabolic pathway, and the flavonoid metabolic pathway (Marchiosi et al., 2020). However, the biosynthetic pathways of phenolic acid derivatives have not been clearly explained. By analyzing the metabolic pathways of key time nodes related to HBA, Dactylorhin A, Millitarine, and Coelonin, we see that they may provide a theoretical basis for exploring the synthetic methods for phenolic acid derivatives.
In this study, aided by the full-length transcriptome sequencing of PacBio Sequel and Illumina short-read sequencing technology (Zhou et al., 2021), we first profiled gene expression in different growth stages of the lag phase, exponential phase, deceleration phase, and decline phase of the suspension culture system of B. striata. Then, we screened out the genes involved in the biosynthesis of phenolic acid derivatives. Finally, we verified the expression patterns of selected genes in the suspension system. The obtained conclusions may provide a theoretical foundation for the genetic breeding of B. striata as well as an experimental basis for exploring the molecular biosynthesis mechanism of phenolic acid derivatives.

Material
The capsules were harvested from the B. striata Germplasm Garden of Zunyi Medical University, Xinpu District, Zunyi City, Guizhou Province, China (27 • 42' N, 107 • 01' E). After disinfection , matured seeds obtained from capsules were induced for the following suspension culture. The information on the compound standards used for measuring the contents of phenolic acid derivatives is shown in Supplementary Table S3. In addition, the compound structures are shown in Supplementary Figure S1.

Culture of Suspension Cells
Suspension culture was carried out for the efficient biosynthesis of phenolic acid by selecting B. striata with relatively uniform growth and good induction. Based on the contents of various compounds, we used MS (Murashige & Skoog Medium) as the base medium, with the addition of 1.0 mg/L of 6-BA, 3.0 mg/L of 2,4-D, and 30.0 g/L of sucrose. The pH was adjusted to 5.8-6.0. After 45 days of induction and 30 days of subculture, the content of secondary metabolites of suspension-cultured cells was measured every 3 days for 45 days postinoculation (Dpi). In addition, the content of phenolic acid derivatives was calculated according to the standard curve (Pan, 2019). The suspensioncultured cells were stored at 3 Dpi,9 Dpi,18 Dpi,27 Dpi,30 Dpi,33 Dpi,36 Dpi,39 Dpi,and 42 Dpi for measuring the chemical content and RNA-Seq sequencing as described below.

Content Determination of Phenolic Acid Derivatives
Starting from the day of the directional suspension culture, the contents of phenolic acid derivatives were measured every 3 days, and the sampling was repeated three times. After being drained and dried, the cell clusters were crushed with a mortar, and 0.2 g was accurately weighed and placed in a round-bottom flask. Then, the cells were heated and refilled with 70% methanolwater for 2 h, transferred to a 5-ml volumetric flask filled with 70% methanol-water to maintain a constant volume, shaken well, and filtered using a 0.45 µm filter membrane for subsequent experiments. The content of phenolic acid derivatives in cells was determined by HPLC (Waters e2695).

Transcriptome Data Assembly and Sequencing
After grinding the selected callus materials of 10 time points with liquid nitrogen, total RNA was extracted by Trizol and stored at −80 • C for subsequent sequencing on PacBio Sequel and Illumina Hiseq 2500 platforms, which were entrusted to us by Novogene Co., Ltd. To obtain clean reads, high-throughput sequencing from a large amount of raw data was obtained, from which the adaptor sequence, primer sequence, and low-quality reads were removed. The LoRDEC (Salmela and Rivals, 2014) software was used to correct the PacBio Sequel data with the Illumina Hiseq 2500 data. The CD-HIT (Fu et al., 2012) software was used to de-redundant the corrected consensus sequence, and the obtained transcript set was applied as the reference sequence (ref) for subsequent analysis. Then, the clean reads of each sample obtained by Illumina sequencing were compared to the reference to determine the transcripts detected in each sample.

Screening of DEGs
Using the RSEM (Li and Dewey, 2011) software to carry out the statistical analysis on the comparison results of Bowtie2 and perform FPKM conversion to determine the expression levels of the genes detected in the samples. The resulting q-values were adjusted by Benjamini and Hochberg's method. DEGs among different samples were screened based on |log 2 (FoldChange)|≥ 2 and a q-value < 0.05.

GO Enrichment and KEGG Pathway Analysis of DEGs
We used GO (http://www.geneontology.org/) and KEGG databases (http://www.genome.jp/kegg/) to perform functional enrichment analysis on the selected DEGs. It was considered that the GO function or KEGG pathway was significantly enriched (q-value < 0.05), and the detection was annotated to DEGs in different GO functional groups or metabolic pathways.

Screen the key Genes Involved in Regulating the Synthesis of Phenolic Acid Derivatives
Using the K-means cluster analysis method, the DEGs with similar expressions were grouped into one category, and the differential gene expression trend was analyzed. Using the R package to analyze the correlation between DEGs and phenolic acid derivatives of HBA, Dactylorhin A, Militarine, and Coelonin at 10 time points (Spearman correlation coefficient), the screening criteria were q-value < 0.05 and R > 0.05. Subsequently, according to the accumulation and change curves of the four compounds on different days, DEGs at different time points were analyzed by the KEGG pathway and GO function analysis. The DEGs with significant upregulation or downregulation were screened. The main metabolic pathways involved in DEGs were determined. Finally, combined with correlation analysis, the key genes involved in the synthesis of phenolic acid derivatives were screened.

QPCR Verification of DEGs
The detected DEGs closely related to each derivative were randomly selected for verification via qPCR. Approximately 1 µg total RNA was used as the template for synthesizing the first cDNA strand by reverse transcription. The SYBR PrimeScriptTM RT-PCR Kit (TaKaRa) fluorescence quantitative kit was used to carry out qPCR reactions on an AGS 4800 real-time qPCR instrument (Hangzhou Anyang Technology Co., Ltd., Hangzhou, China). Primer Premier software was used to design specific primers for qPCR reactions of different randomly selected genes ( Table 1). The approaches for cDNA reverse transcription and the qPCR reaction system are detailed in Supplementary Tables S5, S6. Besides, the cDNA reverse transcription reaction conditions were 42 • C for 60 min and 95 • C for 3 min. In addition, qPCR reaction conditions were enzyme activation at 95 • C for 5 min; denaturation at 95 • C for 5 s, and extension at 60 • C for 30 s, 45 cycles.

Data Statistical Analysis
Each experiment was conducted with three biological replicates and three technical replicates. The 2 − Ct method was used to analyze gene expression levels in qPCR arrays. The statistical analysis of the data was FIGURE 1 | The content changes of four phenolic acid derivatives. * Represents the significant difference between the two groups compared with the previous group at a level of P < 0.05. ** represents the significant difference at a level of P < 0.01. run in SPSS, and the statistical methods of the one-way ANOVA and multiple comparisons were applied to the harvested data.

Determination of the Content of Phenolic Acid Derivatives in Suspension Cells of B. striata
The content of the four metabolites varied significantly in different growth stages of suspension culture cells (Figure 1).
With the growth of culture, the accumulation of HBA in the suspension culture cells increased first and then decreased before finally reaching the maximum peak at 33 Dpi.

Transcriptome Sequencing Analysis of Suspension Cells of B. striata
By performing full-length transcriptome sequencing on suspension cells of B. striata at different time points, a total of 50.02 Gb of raw data were obtained, which were deposited in the SRA database under the accession number SRR18045794. A total of 902,688 CCS (Circular Consensus Sequence) readings were obtained, including 719,981 reads with 5'-primers, 785,110 reads with 3'-primers, and 781,562 reads with Poly A tail. There were 620,018 full-length (FL) reads. The number of full-length non-chimeric reads (FLNC) was 472,211, and the average length of FLNC was 2,423 bp. The FLNC sequence of the same transcript was converted using the ICE algorithm Clustering, which was performed to eliminate redundancy, and a total of 246,933 identical transcripts were obtained. Using CD-HIT to perform de-redundancy analysis on high-quality transcripts, a total of 100,276 transcripts were obtained, which were used as the final transcript sequence for further analysis .

Distribution of DEGs
In order to explore the expression trend of DEGs during the development of suspension cells of B. striata, we conducted a Kmeans analysis on DEGs, which resulted in four gene expression trend graphs (Figure 2A). Although gene cluster 3 contained the largest number of genes, the changes in expression levels in each period were not obvious. The gene expression levels of gene cluster 2 were upregulated throughout the development of suspension cells of B. striata, whereas gene clusters 1 and 4 were reduced.
The FPKM value was used to quantify gene expression, and it can be seen from the results that there is little difference in gene expression levels among different samples ( Figure 2B). This showed that, after FPKM quantification, the analysis results of DEGs between each sample are accurate and reliable. Subsequently, DESeq2 software was used to analyze the significant difference in the DEGs between different samples, and pairwise comparisons of suspension cells at eight time points were made. It was found that the number of DEGs gradually increased during the growth and breeding of suspension cells. A total of 2,052 DEGs were selected from the 3 Dpi vs. 30 Dpi comparison, which was significantly higher than the other comparisons. In the following comparison groups (3 Dpi (Figure 3), indicating that these DEGs may play an important role in the growth and development of B. striata.

Correlation Analysis Between DEGs and Four Main Secondary Metabolites at Each Time Point
We detected correlations between the four phenolic acid derivatives and 17,500 DEGs. Among them, HBA, Dactylorhin A, Militarine, and Coelonin were significantly correlated with 5,134, 3,567, 1,687, and 599 DEGs, respectively (P < 0.05) (Supplementary Table S1).

DEGs Analysis of HBA
Based on the change in HBA content, we found that it increased gradually from 0 Dpi to 27 Dpi, and then the HBA content began to increase exponentially, where it stabilized at 33 Dpi. Three time points (18 Dpi, 27 Dpi, and 33 Dpi) with significant differences in HBA content were selected for differential expression analysis. At the 18 Dpi vs. 27 Dpi time point, the number of upregulated genes was 78, and the number of downregulated genes was 127. At the 27 Dpi vs. 33 Dpi time point, the number of upregulated genes was 712, while the number of downregulated genes was 209.
The functional clustering arrays of GO and KEGG on the DEGs imply a mechanism for explaining how HBA was biosynthesized. GO function analysis found that, at the 18 Dpi vs. 27 Dpi time point, DEGs were mainly enriched in the binding process, cellular process, metabolic process, and catalytic activity. At the 27 Dpi vs. 33 Dpi time point, DEGs were mainly enriched in the hydrolase activity and nucleoside phosphatase activity. A KEGG pathway analysis of DEGs at three time points found that a total of 39 DEGs were significantly enriched in four pathways (q-value < 0.05), namely  Table S2). We selected the pathways of phenylpropanoid biosynthesis and cyanoamino acid metabolism, which are closely involved in the synthesis of HBA. There were 11 upregulated genes with significant differences in the phenylpropane metabolic pathways, encoding β-glucosidase (BGLU) and peroxidase (PRX); there were seven downregulated genes encoding phenylalanine ammonia-lyase (PAL) and BGLU,. In the cyanoamino acid metabolism pathway, the up and downregulated genes with bigger differences were involved in encoding BGLU. Compared with the results of correlation analysis, 18 unigenes (8 enzymes) were finally screened out as being potentially involved in HBA synthesis ( Table 2).

DEGs Analysis of Dactylorhin A
Three time points−18 Dpi, 36 Dpi, and 42 Dpi-were detected to have significant differences with Dactylorhin A content, and these time points were selected for further differential expression analysis.  metabolic process and oxidoreductase activity. While the KEGG pathway analysis found that, at the 18 Dpi vs. 36 Dpi time point, a total of 51 DEGs were significantly enriched in four pathways (qvalue < 0.05), namely cyanoamino acid metabolism, starch and sucrose metabolism, phenylpropanoid biosynthesis, and alanine, aspartate, and glutamate metabolism pathways. In addition, the DEGs at the 36 Dpi vs. 42 Dpi were mainly annotated for the TCA cycle and starch and sucrose metabolism, among others (Supplementary Table S2). Selecting the metabolic pathways of starch and sucrose, which were closely related to the synthesis of Dactylorhin A, we found that there were 13 upregulated genes with significant differences, encoding UGP, BGLU, SPS, TPS, scrK, and 13 downregulated genes, encoding INV, SUS, TPS, BAM, and BGLU. Combined with the above analysis, 23 unigenes (related to 11 enzymes) were screened out as being involved in the synthesis of Dactylorhin A ( Table 2).

DEGs Analysis of Coelonin
Based on the change in the content of Coelonin, we found that its content increased slowly from 3 Dpi to 30 Dpi and to 33 Dpi, remained relatively stable, and then increased exponentially before reaching its peak at 39 Dpi.  The functional clustering arrays of GO and KEGG on the DEGs imply a mechanism for explaining how Coelonin was biosynthesized. Through the comparative analysis of the GO database, the unigenes annotated by the GO can be divided into three categories: biological process, molecular function, and cellular component. At the 3 Dpi vs. 30 Dpi time point, DEGs were mainly involved in the metabolic process, cellular process, and organic substance metabolic process in biological processes; hydrolase activity, catalytic activity, and binding in molecular function. In the suspension cells of the 30 Dpi vs. 33 Dpi and 33 Dpi vs. 39 Dpi, DEGs were significantly enriched in the plant-type cell wall organization, plant-type cell wall organization, cell wall organization, and external encapsulating structure organization in biological processes; significantly enriched molecular functions were the structural constituents of the cell wall, such as structural molecule activity and others. The cell component that was significantly enriched was the histone deacetylase complex. Among them, the unigenes that are more abundant in biological processes may be related to the pathway of Coelonin synthesis concerned in this study.
KEGG enrichment results showed that, in comparison with 3 Dpi, 53 DEGs were involved in the metabolic pathways of starch and sucrose metabolism (ko00500) Table S2).
We selected the phenylpropanoid biosynthesis, phenylalanine metabolism, and cyanoamino acid metabolism pathways, which were closely related to the synthesis of Coelonin. And according to the correlation analysis between gene expression and Coelonin content, the differential genes related to Coelonin metabolism were found. The results showed that the key regulatory genes mainly encoded PAL, 4CL, POD, SUS, GAUT, ADH, PDC, PFK, and ALDO. Finally, 14 unigenes (9 enzymes) were identified as potential Coelonin synthesis candidates ( Table 2).

DEGs Analysis of Militarine
According to the change in the content of Militarine, it was found that from 3 Dpi to 30 Dpi, its content increased slowly and then increased exponentially before reaching its peak at 39 Dpi and then dropping sharply. Pairwise comparison of the DEGs of 3 Dpi, 30 Dpi, 39 Dpi, and 42 Dpi (3 Dpi vs. 30 Dpi, 30 Dpi vs. 39 Dpi, 39 Dpi vs. 42 Dpi) showed that there were 2052, 652, and 589 DEGs, respectively, and 1825, 496, and 410 DEGs appeared in only one DEG set, respectively; 33, 48, and 182 DEGs appeared in two DEG sets, respectively, and a total of 12 DEGs appeared in 3 DEG sets. Among them, at the 30 Dpi vs. 39 Dpi time point, the number of upregulated genes was 296, and the number of downregulated genes was 356; at the 39 Dpi vs. 42 Dpi time point, the number of upregulated genes was 387, and the number of downregulated genes was 202.
According to GO analysis, 1,344 DEGs were annotated into the GO database at 3 Dpi vs. 30 Dpi time point, with the majority of them enriched in hydrolase activity, asparagine synthase (glutamine-hydrolyzing) activity, catalytic activity, single-organism metabolic process, and carbohydrate metabolic process, among others. At the 30 Dpi vs. 39 Dpi time point, 476 DEGs were significantly enriched in cysteine-type peptidase activity, covalent chromatin modification, and histone modification. The KEGG pathway analysis of DEGs at three time points found that a total of 136 DEGs were significantly enriched in five pathways at 3 Dpi vs. 30 Dpi time point, namely starch and sucrose metabolism, phenylpropanoid biosynthesis, glycolysis/gluconeogenesis, photosynthesis, and phenylalanine metabolism pathways; at the 30 Dpi vs. 39 Dpi time point, a total of 30 DEGs were significantly enriched in Spliceosome; at the 39 Dpi vs. 42 Dpi time point, DEGs were mainly enriched in plant hormone signal transduction and phenylpropanoid biosynthesis (Supplementary Table S2). Based on the above analysis, 41 unigenes (18 key enzymes) of the main regulatory genes were screened out, which may be involved in the synthesis of Militarine, as shown in Table 2.

Expression Analysis of Related Genes Regulating the Synthesis of Phenolic Acid Derivatives
In order to verify the accuracy of gene expression in the transcriptome sequencing results, 10 DEGs were randomly selected from samples of four periods and verified by qPCR, in which the β-actin was applied as the internal reference gene. The results showed that the relative gene expression levels were consistent with the transcriptome sequencing results, indicating that the data obtained by transcriptome sequencing were accurate and reliable (Figure 5).

DISCUSSION
The secondary metabolites were produced by plants during their growth and development, as well as to defend against environmental stresses. In addition, the synthesis and accumulation of these metabolites can be significantly improved through suspension cell culture (Yang et al., 2018;Isah, 2019;Jiang et al., 2020). The secondary metabolites in suspension cells provide the material basis for their growth and reproduction. We analyzed the DEGs in different growth stages and found that many genes involved in the biosynthesis of secondary metabolites expressed significant differences in different growth stages. Moreover, the content of secondary metabolites also changed significantly, indicating that, with the change in suspension cell morphology and continuous growth and development, the synthesis of secondary metabolites of suspension cells is regulated by a variety of genes.
As time goes on, cultured suspension cells synthesize and accumulate secondary metabolites through different metabolic processes, thereby improving their ability to synthesize and accumulate HBA, Dactylorhin A, Militarine, and Coelonin in suspension cells of B. striata (Pan, 2019). This study found that the contents of HBA at 18 Dpi,27 Dpi,and 33 Dpi,Dactylorhin A at 18 Dpi,36 Dpi,and 42 Dpi,Coelonin at 3 Dpi,30 Dpi,33 Dpi,and 39 Dpi,and Militarine at 3 Dpi,30 Dpi 39 Dpi, FIGURE 5 | Comparison of expression of 10 genes obtained by qPCR analysis and by RNA-seq (FPKM values). * Represents the significant difference between the two groups compared with the previous group at a level of P < 0.05. ** represents the significant difference at a level of P < 0.01. and 42 Dpi were significantly different by HPLC, respectively. It was found that, compared with the high-efficiency culture lag zone, the content of the four metabolites in the exponential phase increased significantly; the content of Dactylorhin A, Militarine, and Coelonin decreased significantly in the declining phase, while HBA slowly increased. According to the chemical structures of these four compounds (Supplementary Figure 1), the main structures of Dactylorhin A and Militarine were the same, which implies that they may have connections in the biosynthesis process. The increase in Militarine content resulted in a reduction of Dactylorhin A synthesis. HBA was simple in structure, which could be an integral part of forming Militarine and Coelonin. Therefore, HBA contributes to the increase of Militarine and Coelonin.
Genes involved in various metabolic processes of cultured suspension cells are critical in the process of the synthesis of secondary metabolites. The metabolism of phenolic acid compounds was a relatively complicated process. The shikimic acid pathway, the phenylpropanoid metabolic pathway, and the flavonoid metabolic pathway were the main pathways for the synthesis of phenolic acid compounds (Saxe et al., 2021). Among them, phenylalanine and tyrosine are direct precursors of phenolic acid biosynthesis, which are converted from glucose into a multistep metabolism via the shikimic acid pathway and then into the phenylpropanoid metabolic pathway (Feduraev et al., 2020;El-Shora et al., 2021). Additionally, phosphoenolpyruvate and erythrose-4-phosphoric acid were the main precursor substance for the synthesis of phenolic acid compounds (Deng et al., 2020). Phenylalanine was deaminated and transformed into trans-cinnamic acid by phenylalanine ammonia-lyase (PAL) (Trivedi et al., 2022). Subsequently, phenolic acids such as trans-cinnamic acid, pcoumarate, and ferulic acid are formed under the action of cinnamate 4-hydroxylase (C4H), 4-coumaratecoa-ligase (4CL), O-methyltransferase (OMT), and other enzymes catalyzed further into caffeic acid and chlorogenic acid (Wohl and Petersen, 2020). HBA was one of the main components of phenolic compounds and the main precursor of aromatic compounds (Zhang Z. L. et al., 2020). Former studies have confirmed that Dactylorhin A and Militarine were glycoside compounds isolated from plants, which were produced by the loss of water from the hydroxyl amino thiol group of monosaccharides or oligosaccharides and the hemiacetal hydroxyl group of another molecule. With the growth and development of suspensioncultured cells, the expression levels of genes involved in the metabolism of starch and sucrose and the biosynthetic pathway of phenylpropanoid differed significantly. The differences in the expression of these genes may be the main reason for the significant differences in the content of glycosides in different periods. BGLU is a group of heterogeneous hydrolases that have the ability to hydrolyze and synthesize glycosidic bonds and have great application value in various fields (Santos et al., 2018;Sun et al., 2020). The enhanced expression of genes involved in BGLU may promote the synthesis of carbohydrates and provide energy for the growth and development of suspensioncultured cells. At the same time, BGLU was regulated by a large number of genes at different time points and was involved in the synthesis of phenolic acids, which indicated that BGLU may be involved in the synthesis of phenolic acids and other substances and maintaining cell morphology during the development of B. striata. Coelonin was a dihydrophenanthrene compound. The methoxy and hydroxyl groups were the main substituents on the benzene ring of the dihydrophenanthrene compound. From the point of view of molecular structure, phenolic acid and its derivatives can be obtained after the oxidation of dihydrophenanthrene compounds.
In this study, we detected several interesting co-changes between gene expression and chemical content, which may infer the genes' function. The expression of i1_LQ_BSbp1_C76772/f1p0 increased steadily from 0 to 8 Dpi but decreased dramatically later, which was positively correlated with the content changes of Dactylorhin A. The expression of i3_LQ_BSbp1_C942 /f1p19 was continuously downregulated from 0 to 36 Dpi but upregulated from 36 Dpi to 39 Dpi and then downregulated again, which was positively correlated with the content change of Dactylorhin A. The expression of i3_LQ_BSbp1_C5868/f1p2 was continuously downregulated at the stage of 0-18 Dpi but then upregulated, which showed a positively significant correlation with the content of Militarine. Based on these results, we predict that the TPS plays a very important role in plant metabolic regulation. In addition, we also speculated that these genes are closely related to the biosynthesis of Dactylorhin A and Militarine. Similar to this analysis and speculation, we built the relationship between genes i2_HQ_BSbp1_C5228/f3p1 and Dactylorhin A and gene i2_LQ_BSbp1_C117598/f1p12 and Militarine. In addition, the INV may be key to the synthesis of Dactylorhin A and Militarine. Therefore, during the growth of B. striata suspension culture cells, with the passage of time, the expression of key enzyme genes regulating the synthetic biological pathway is different in different periods, and through interaction, the content of phenolic acid derivatives is different in different periods.
In different growth periods, the metabolic pathways involved in the synthesis of these four secondary metabolites of DEGs mainly include phenylpropanoid biosynthesis, phenylalanine metabolism, cyanoamino acid metabolism, and starch and sucrose metabolism. The results of the enrichment of metabolic pathways indicated that carbohydrate metabolism was needed as a basis during development, and hormones and some secondary metabolites were required to co-regulate. According to the KEGG pathway, GO enrichment analysis, and correlation analysis of DEGs, the final screening of enzymes involved in the synthesis of HBA, Dactylorhin A, Militarine, and Coelonin were 8, 11, 18, and 9, respectively, and the related unigenes were 18, 23, 41, and 14, which led to the speculation that these genes which were screened out may be the key regulators of phenolic acid derivatives. Therefore, the differential expression of key enzyme genes in the biosynthetic pathway at different periods may ultimately lead to differences in the content of phenolic acid derivatives at different periods.
At present, four phenolic acid derivatives of secondary metabolites have been isolated and identified, but the metabolic mechanism of their synthesis and accumulation has not been studied. It was found that the content of Dactylorhin A, Millitarine, and other compounds in B. striata was higher than the cultivated ones (Wang et al., 2014). In this study, gene expression profiles of suspension cells at different growth stages were obtained through the transcriptome method, and a total of 17,500 DEGs were screened. The preliminary screening showed that the genes of PAL, 4CL, BGLU, INV, TPS, BAM, and PDC may regulate the key enzyme of phenolic acid derivatives. Through the gene identification and regulatory mechanism exploring of four main secondary metabolites, the key regulated genes in this study could be overexpressed or silenced in the future to further verify their functions in the biosynthesis of four phenolic acid derivatives in B. striata. Thus, the metabolic mechanism of synthesis and accumulation of the four secondary metabolites will be thoroughly explained. This provided data for an in-depth study of the synthesis and regulation mechanisms of phenolic acid derivatives in suspension cells of B. striata.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: National Center for Biotechnology Information (NCBI) BioProject database under accession numbers PRJNA807395 and PRJNA807396.

AUTHOR CONTRIBUTIONS
DX conceived, supervised, wrote, and reviewed the manuscript. HL, CH, WW, and LL originally wrote and reviewed the draft. HL, CH, QL, MW, SX, JS, and YH conducted the experiments and carried out the analysis. DX founded and administrated the project. All authors read and approved the final version.