Transcriptome sequencing of Salvia miltiorrhiza after infection by its endophytic fungi and identification of genes related to tanshinone biosynthesis

Abstract Context: Salvia miltiorrhiza Bunge (Labiatae) is a traditional Chinese herb. Endophytic fungi, which are biotic elicitors, can induce accumulation of secondary metabolites in their host plants. Objective: To analyze the interaction mechanism between S. miltiorrhiza and endophytic fungi. Materials and methods: Endophytic fungi U104 producing tanshinone IIA were isolated from the healthy disease-free tissue of root of S. miltiorrhiza by conventional methods. The endophytic fungus U104 of S. miltiorrhiza was co-cultured with the sterile seedlings of S. miltiorrhiza for 20 d (temp:day/night = 26 °C/18 °C, photoperiod:12/12 h, illuminance:2000 Lx). Transcriptome sequencing of S. miltiorrhiza seedlings after 20 d of co-cultivation was performed using the Illumina platform. Results: A total of 3713 differentially expressed genes (DEGs) were obtained. These different expression genes, such as STPII, LTP2, MYB transcription factors, CNGC, CDPK, Rboh, CaM, MAP2K1/MEK1, WRKY33, SGT1/SGT and Hsp90/htpG, showed that host S. miltiorrhiza had biological defence response in the initial stage of interaction. Under the induction of endophytic fungi, 14 key enzyme genes were up-regulated in the tanshinone biosynthesis pathway: DXS, DXS2, DXR, HMGR3, AACT, MK, PMK, GGPPS2, GPPS, KSL, IDI, IPII, FDPS and CPS. Discussion and conclusions: A total of 14 key genes were obtained from the tanshinone component synthesis and metabolic pathways, providing a reasonable explanation for the accumulation of tanshinone components, an accumulation induced by endophytic fungi, in the host plants. The large amounts of data generated in this study provide a strong and powerful platform for future functional and molecular studies of interactions between host plants and their endophytic fungi.


Introduction
As a traditional Chinese herb, Salvia miltiorrhiza Bunge (Labiatae), has more than 1000 years of application history. Its active medicinal ingredients mainly include two groups: hydrophilic phenolic acids and lipid-soluble tanshinones. The phenolic acids possess various bioactivities including antioxidant, antiinflammatory, anticancer, antibacterial, antivirus, antifibrotic activities . Phenolic acids such as salvianolic acid B Jing et al. 2016) and salvianolic acid A (Pan et al. 2014;Ding et al. 2016) have antitumor and antioxidant activities. Tanshinones have antitumor (Dong et al. 2011) and anti-inflammation  activities; cryptotanshinone can prevent and treat atherosclerosis (Suh et al. 2006), tanshinone IIA can be used for treating osteoporosis and reducing blood lipids (Kwak et al. 2006) and cholesterol content . Tanshinone IIA is a major lipid-soluble compound having promising health benefits. In vivo and in vitro studies showed that the tanshinone IIA and salvianolate have a wide range of cardiovascular and other pharmacological effects, including antioxidative, anti-inflammatory, endothelial protective, myocardial protective, anticoagulation, vasodilation and antiatherosclerosis, as well as significantly reducing proliferation and migration of vascular smooth muscle cells (Ren et al. 2019).
S. miltiorrhiza has many biological activities and excellent prospects. The majority of studies focus on using genetic engineering methods to regulate the expression of key genes in secondary metabolic pathways to directly influence the accumulation of end products. Studies have shown that overexpression of SmGGPPS and SmDXSII in hairy roots produces higher levels of tanshinone than control and single-gene transformed lines . Two elicitor treatments suggested that tanshinone accumulation positively correlated to the expression of key genes such as SmGGPPS, SmCPS and SmKSL (Hao et al. 2015); In the hairy roots of S. miltiorrhiza, overexpression of SmERF1L1 significantly increased tanshinones production by comprehensively upregulating tanshinone biosynthetic pathway genes , silencing of SmERF115 reduced the phenolic acid level, but increased tanshinone content ). In addition, elicitors such as methyl jasmonate, salicylic acid, heavy metal ions (Co 2þ , Ag þ and Cd 2þ ), sorbitol and ultraviolet can be used to increase the accumulation of tanshinones and phenolic acids in S. miltiorrhiza hairy roots (Zhao et al. 2010;Xing et al. 2014;Wang et al. 2016), whereas, biotic elicitors such as yeast extracts can also stimulate its hairy roots to produce more tanshinones (Shi et al. 2007;Wu et al. 2008).
Plant endophytic fungi refer to fungi that live inside the various tissues and organs of healthy plants during certain stages or all stages of their life cycle without causing apparent symptoms of infection in host plants. During the long-term evolutionary process, endophytic fungi are important components in medicinal plants. They form a stable and mutually beneficial symbiotic relationship with medicinal plants and can produce the same, or similar, secondary metabolites in the host plants (Venugopalan and Srivastava 2015). Endophytic fungi can also act as elicitors to rapidly activate specific genes in the secondary metabolic pathway in medicinal plants to accumulate a large amount of active ingredients (Zhai et al. 2017).
RNA sequencing (RNA-Seq) is the high-throughput sequencing of mRNA in a species. Its resolution has the accuracy of a single nucleotide, it can dynamically reflect gene transcription levels and it provides specific sequence-structure information of transcripts in samples (Hansen et al. 2011). Currently, RNA-Seq is being extensively applied in all fields, including basic biological research, medical research and drug development (Kawahara et al. 2012;Foth et al. 2014;Zhang et al. 2014). This study performed RNA-Seq on sterile plantlets of S. miltiorrhiza and endophytic fungi to examine the differential gene expression after infection of tissue-cultured plantlets with endophytic fungi, to understand the underlying molecular mechanism of interaction, then analyzed to provide new ideas and methods for studying the regulation of secondary metabolism in medicinal plants.
In the laboratory, two endophytic fungi producing tanshinones were isolated from the roots of S. miltiorrhiza: TR21 and U104. TR21 is a wild strain isolated from S. miltiorrhiza and U104 was induced by TR21. The strains of TR21 and U104 were identified as Ascomycota, Eurotiomycetes, Eurotiales, Trichocomaceae, Eurotium, Emericella foeniculicola Udag (Trichocomaceae). In the early stage of the experiment, the endophytic fungi TR21 and U104 of S. miltiorrhiza were co-cultured with the sterile S. miltiorrhiza seedlings for 10 d and 20 d, respectively, and the content of tanshinone in the plants was determined. The accumulation amount of tanshinone was used as the index to screen the best inducible strain and time. The results showed that the content of tanshinone in aseptic tissue culture seedlings of S. miltiorrhiza was 103.89, 151.08 and 155.56 mg/mL, respectively, after 10 d of co-culture with endophytic fungi TR21 and U104. After a total of 20 d of culture, the contents of tanshinone in the control, TR21 and   The content of tanshinone increased to different degrees. The effect of 20 d interaction was significantly better than that of 10 d interaction, and in the process of 20 d interaction, the induction effect of U104 strain was significantly higher than that of TR21 strain. Taking the content of tanshinone as an indicator, it was found that when the sterile seedling of S. miltiorrhiza was treated with U104 strain for 20 days, the induction effect was the best.

Materials and methods
Endophytic fungus and host plant materials S. miltiorrhiza sterile plantlets were obtained by tissue culture and identified by Professor Zhezhi Wang of Shaanxi Normal University, Xi'an, Shaanxi, China. S. miltiorrhiza sterile plantlets that had been preserved in the laboratory and cultured for one week were used as the materials. The test strain U104 is an endophytic fungus, which was isolated from S. miltiorrhiza and preserved in our laboratory. The stem segments of sterile plantlets with opposite leaves were inserted in a triangular configuration directly into culture flasks containing 100 mL of MS culture medium, with three plants in each flask. The day/night culture temperature was 26 C/18 C, the light cycle was 12 h, and the luminance was 2000 lÂ. The U104 strain was inoculated into potato dextrose agar (PDA) culture medium and cultured at 28 C for 7 d. The sterile plantlets (which had been cultured for 7 d, were 2-3 cm in height and showed consistent development status) were divided into two groups. In one group, one piece of fungal disc with a diameter of 3 ± 0.1 mm was inoculated into culture medium at 1-2 cm from the plantlets. In the other group, sterile plantlets without inoculation with fungi were used as the controls. The whole plant was taken as a sample in this study. After 10 d and 20 d of co-culture, the samples were immediately frozen in liquid nitrogen and stored in a À80 C freezer for future use. All interaction experiments had two biological replicates. The results are shown in Figure 1.

RNA extraction and mRNA-Seq library construction and sequencing
Total RNA in samples from two groups was extracted using a MiniBEST Plant RNA Extraction kit (TaKaRa, Dalian). After RNA samples were qualified using electrophoresis, NanoDrop, Qubit 2.0 and Agilent 2100 analyses, the cDNA library was constructed by Biomarker Technology Co., Ltd (Beijing, China). Paired-end sequencing was performed using an Illumina HiSeq 2000 sequencer and the read length of sequencing was determined with PE150.

Sequence processing and unigene library
After the sequencing was complete, the linker sequences and low-quality reads of raw data were removed to obtain high- Figure 2. The result comparison between qRT-PCR and RNA-Seq. The expression of Smc34207, Smc36046 and Smc43522 genes was down-regulated, and the remaining genes were up-regulated, which was consistent with the RNA-Seq results. Table 3. Primer sequences used for qRT-PCR amplification.

Gene
Primer sequences ATCGGCATTCCACAGACT CTTACATCCTCCACACCAAT quality clean data. Trinity software (Grabherr et al. 2011) and paired-end method was used for sequence assembly. BLAST software (Altschul et al. 1997) was used to compare unigene sequences in the NR, Swiss-Prot (Apweiler et al. 2004), GO (Ashburner et al. 2000), COG (Tatusov et al. 2000), KOG (Koonin et al. 2004), KEGG (Kanehisa et al. 2004) and Pfam databases (Finn et al. 2014) to obtain the annotation information of the unigenes. BLAST parameter E-value is not greater than 10 À5 . The differential expression among samples was analyzed using DESeq (Anders and Huber 2010). A false discovery rate (FDR) less than 0.01 and a fold change (FC: the ratio between base mean value of the treatment group and base mean value of the blank group) no less than 2 were used as the screening standards. In addition, r 2 (Schulze et al. 2012) was used as the evaluation indicator for the correlation between the biological replicates.

Quantitative real-time RT-PCR
Fluorescence quantitative polymerase chain reaction (PCR) was performed using SYBRV R Premix Ex Taq TM II (Tli RNaseH Plus).
The amplification system included 1 mL of cDNA (10Â), 0.5 mL of each of the 10 mmol/L upstream and downstream gene-specific primers and 10 mL of 2 Â SYBR Premix Ex Taq TM , with the total  volume brought to 20 mL using ddH 2 O. The amplification condition list is shown in Table 1. Each sample had three technical replicates and three biological replicates. Data analysis was performed using the 2 ÀDDCt method.

RNA-Seq analysis of transcriptome samples
After data filtering and quality-control analysis of the raw data, a total of 22.73 Gb clean data were obtained. At least 5.42 Gb clean data were attained for every sample. The percentages of Q30 bases in all samples were no less than 91.46% (Table 2). Sequence assembly was performed and a total of 200,043 transcripts and 96,802 unigenes were obtained, and 18,364 unigenes had lengths longer than 1 kb. Through the comparison of seven major databases, a total of 48,621 unigenes with annotation information were ultimately obtained. Based on these results, 3713 differentially expressed genes (DEGs) were obtained using the differential expression analyses. A total of 3451 genes showed up-regulated expression and the FC was more obvious at 5-to 10-fold. A total of 262 genes showed down-regulated expression, and the FC was generally between 1-and 5-fold. The scatter plot of the gene correlation of the samples showed that the expression trends of most genes in the samples of the two biological replicates were similar, and the correlation between the replicates was high (both r 2 > 0.90). Actin was used as the internal control gene and the relative expression levels of all genes obtained using qRT-PCR were compared and validated with the FC values in the RNA-Seq results ( Figure 2). The upstream and downstream gene-specific primers used in qRT-PCR are shown in Table 3. The results were basically consistent, indicating the reliability of the RNA-Seq results.

Functional annotation and enrichment analysis of the DEGs
The annotation of the unigenes in the three GO categories showed that nodes with more obvious differences might be associated with the observed differential expression. GO enrichment analysis on unigenes with differential expression showed that the DEGs involved in the oxidation-reduction process accounted for most. Additionally, 99 genes were involved in oxidoreductase activity. Because the late stage of biosynthesis of tanshinone components involves a large amount of oxidation-reduction reactions, these genes could be used as potential candidate genes of key enzymes in the biosynthesis pathway of tanshinone components. In the cellular components, in addition to those undergoing significant enrichment in the cell nuclei, cytoplasm and ribosomes, genes enriched in the integral components of the membrane and the extracellular region also accounted for a large quantity, this might be associated with recognition proteins (receptor proteins) on the cell surface during interactions between endophytic fungi and the host plant, S. miltiorrhiza. In addition, a fair number of DEGs showed significant enrichment in some transcription and regulatory factors. The significance of this discovery guided our studies on signal transduction during the interaction. The KEGG pathways showed enrichment of DEGs mainly included primary metabolic processes, such as amino acid metabolism, glucose metabolism, lipid metabolism and carbon fixation, as well as secondary metabolic processes, such as phenylalanine metabolism (ko00360), terpenoid backbone biosynthesis (ko00900) and phenylpropanoid biosynthesis (ko00940). In addition, plant-pathogen interactions (ko04626) and plant hormone signal transduction (ko04075) also showed significant enrichment.

Responsive expression of the host plant at the initial stage of induction
In the response to the induction process (Figure 3), which was induced by the endophytic fungi, in the host S. miltiorrhiza, eight differential genes, including CNGC (cyclic nucleotide-gated channel), CDPK (calcium-dependent protein kinase [EC:2.7.11.1]), Rboh (respiratory burst oxidase), CaM (calmodulin), WRKY33 (WRKY transcription factor), MAP2K1/MEK1 (mitogen-activated protein kinase kinase 1 [EC:2.7.12.2]), SUGT1/SGT1 (suppressor of G2 allele of SKP1) and HSP90A/htpG (molecular chaperone HtpG), all showed upregulated expression. CNGC is a nonselective cation channel and is a component of the signal transduction pathway in plant systems (Jha et al. 2016). When plant is induced by its endophytic fungi, the CNGC channels will open and Ca 2þ influx occurs (Verret et al. 2010;Ma 2011). Therefore, on one hand, CaM activation causes feedback inhibition on CNGC activities and prevents a rapid increase of intracellular Ca 2þ concentrations. On the other hand, CDPK activation (Yoon et al. 1999) leads to the phosphorylation of downstream target proteins such as Rboh (Kobayashi et al. 2006;Suzuki et al. 2011).
The WRKY protein is one of the substrates of the mitogenactivated protein (MAP) kinase signalling cascade reaction (Mao et al. 2011;Zhou et al. 2015). Therefore, the WRKY transcription factor can be inferred to activate WRKY regulatory genes, especially defence-related genes. Comprehensive analyses showed that certain defence-response reactions occur in host plants at the initial stage of induction by the endophytic fungi.

Key tanshinone synthesis-related genes and pathways
Among the 3713 obtained DEGs, 75 genes were annotated to the host plant S. miltiorrhiza (Table 4) and could be classified into seven groups: biological stimulus-response, transcription factor MYB, terpenoid synthesis, phenolic acid synthesis, cytochrome P450 and oxidation-reduction reaction and others. The plant transcription factor MYB is one of the largest transcription factor families in plants . They generally serve as positive regulatory factors and exert their functions in stress response and the phenylpropanoid metabolism pathway in plants.
Cytochrome P450 extensively participates in plant development and metabolism regulation (Seki et al. 2011;Li et al. 2013;Wu et al. 2013). The CYP76AH1 annotated in this study is the first P450 gene in the tanshinone biosynthesis pathway (Guo et al. , 2016Ma et al. 2015). The upregulation of its expression might be associated with accumulation of tanshinone components. We speculated that under the induction by U104 endophytic fungi, the plant first developed defence response reactions to upregulate the expression levels of genes encoding lipid transfer protein-2 (LTP2) (Gom es et al. 2003;Wu et al. 2004), allergen, SMLII (Peumans and Van Damme 1995) and other stress response genes. Next, changes in genes encoding the transcription factor MYB and other related proteins were activated. Finally, the expression levels of the genes of key enzymes in the terpenoid biosynthesis pathway changed to promote the accumulation of S. miltiorrhiza active ingredients. The annotated key enzymes such as DXS (Kai et al. 2011;Ma et al. 2012), DXR (Wu et al. 2009), HMGR (Dai et al. 2011;Kai et al. 2011;Ma et al. 2012;Shi et al. 2014), MK, GGPPS (Kai et al. 2011;Ma et al. 2012), GPPS, KSL , IDI, IPII, FDPS and CPS Xu et al. 2015) in S. miltiorrhiza were also present in the terpenoid compound biosynthesis pathway ( Figure  4, Table 5). In addition, AACT and PMK in the mevalonic acid (MVA) upstream pathway also had upregulated expression.
These key enzymes provided significant guidance in our studies on the terpenoid metabolic pathway. In future studies, the molecular functions of these enzymes can be studied using genetic methods such as gene silencing and gene overexpression.

Discussion
By implementing RNA-Seq in host plants, all involved defence response DEGs were exhibited, as well as vital DEGs that promoted the accumulation of active ingredients after induction of endophyte fungi. The large amount of data generated in this study provides a powerful platform for functional and molecular studies of future interactions between host plants and their endophytic fungi. In these DEGs, CNGC, CDPK, Rboh, CaM, MAP2K1/MEK1, WRKY33, SUGT1/SGT1 and HSP90A/htpG are the DEGs that involved in biological response stimulation, in which WRKY33 belongs to the WRKY gene family, one of the largest family of plant transcription factors currently studied (Suttipanta et al. 2011;Phukan et al. 2016;Chen et al. 2017 respond to the induction of MeJA and improve tanshinone production after the induction of S. miltiorrhiza using MeJA, indicating that SmWRKY2 may be involved in stress-regulated processes . SmWRKY1 can respond to the induction of salicylic acid (SA), methyl jasmonate (MeJA) and nitric oxide (NO), and improve the yield of tanshinone by positively regulating SmDXR expression (Cao et al. 2018).
Overexpression of NtWRKY50 upregulated the expression level of related defence genes and increased tobacco resistance to Ralstonia solanacearum . Therefore, we hypothesized that the up-regulated expression of WRKY33 may be a key gene for regulating tanshinone production in response to fungal induction in plants.
In this study, seven key enzymes were up-regulated to varying degrees in the upstream of the tanshinone biosynthesis pathway: MEP pathway (DXS, DXS2, DXR), MVA pathway (AACT, MK, PMK, HMGR3). Previous studies have shown that the biosynthesis of tanshinone is mainly through the MEP pathway whereas our study found that key genes in the upstream pathway of MVA showed significant differences (FC > 6) after U104 endophytic fungi induction. Therefore, we speculate that the upstream pathway of MVA may play a key role in the synthesis of tanshinone precursors in the S. miltiorrhiza seedlings induced by endophytic fungi and the specific regulatory mechanisms remain to be further explored.
Induced by endophytic fungi, the host plant can respond to the fungus and increase the yield of tanshinone by regulating the genes involved in the tanshinone synthesis pathway. Our work reported RNA-seq of S. miltiorrhiza by endophytic fungi. In addition, some studies have compared RNA-seq of S. miltiorrhiza by other elicitors such as methyl jasmonate (MeJA) and yeast extract (YE). MeJA and YE responsive genes related to tanshinones and phenolic acids biosynthesis. Compared to MeJA, YE had a more significant effect on genes involved in biosynthesis of tanshinone. The expression patterns of genes involved in phenolic acid biosynthesis pathways were diverse. PAL, C4H, 4CL, TAT, HPPR, RAS and CYP98A14 were all induced by MeJA. Nevertheless, YE didn't show any clear effect on these genes. It was also consistent with change of active ingredients contents after treatment by these two elicitors, in which the content of tanshinones in hairy root could be induced by these two elicitors, but to phenolic acid, the contents could only be induced by MeJA, not by YE (Zhou et al. 2017). However, the mechanism of interaction between endophytic fungi and host still needs a lot of research and repeated verification owing to the lack of studies on the interaction between fungi and S. miltiorrhiza and studies on the metabolic pathway of tanshinone in endophytic fungi with the same active ingredients in the host.

Disclosure statement
The authors declare that they have no competing interests.

Authors' contributions
The professor Xiying Wei designed the research. Yan Jiang, Lei Wang, Shaorong Lu, Yizhe Xue, Juan Lu and Yanyan Zhang performed the research, contributed to the development of material and PCR and qPCR analysis, then contributed to the writing of the article. All authors read and approved the final manuscript.