Sequencing and De novo assembly
Transcriptome sequencing of D. kotschyi was performed by the Illumina HiSeq technology (HiSeq2500 platform). A total of 23,599,495 raw pair-end reads of 150 bp length were generated. After removal of ambiguous nucleotides, low-quality sequences, and contaminated sequences, all high-quality reads were normalized. Finally, a total of 6,529,853 paired-end reads were used for de novo assembly (Supplementary Table 1). The de novo assembly resulting from trinity contains 67,859 genes and 165,597 transcripts with a GC percentage of 42.64 and an average length of 1066 bp. The size distribution of the transcripts and unigenes was shown in Supplementary Fig. 1a (More details in Supplementary Table 2).
Functional annotation
In this investigation, 165,597 de novo assembled transcripts of the D. kotschyi were compared with TrEMBL, Swiss-Prot, and Pfam (Protein family) databases to assign putative functions and Gene Ontology browser, KEGG, as well as EggNOG databases for further gene evolutionary histories and functional annotations. Comparison of the D. kotschyi transcripts with TrEMBL database as an extensive protein database using the BLASTx search with the stringency of E-value 1e−10 revealed that 104,606 sequence transcripts had high levels of similarity to the ones from related plant sequence (Supplementary Fig. 1b). Among the 104,606 annotated transcripts, 46.20% (46,993 transcripts) showed high similarity to Handroanthus impetiginosus, followed by 25.79% (26,236 transcripts) to Erythranthe guttata, 2.16% (2195 transcripts) to Salvia miltiorrhiza. Abundance estimation of each gene and isoform were calculated based on FPKM (fragments per kilobase of transcript per million mapped reads) and TPM (transcripts per million) values (Supplementary excel file 1). The distribution of TPM values for genes and transcripts is shown in Supplementary Fig. 2.
Moreover, all assembled unigenes were aligned with other public sequence databases, including Swiss-Prot, Pfam, KEGG, Gene Ontology browser, and EggNOG. The numbers of annotated assembled transcripts mapping to each database are summarized in Supplementary Table 3. A total of 47,776 transcripts were annotated in the six databases. A Venn diagram displayed shared and unique unigenes of D. kotschyi according to Swissprot, Pfam, EggNOG, GO, and KEGG public databases (Supplementary Fig. 3a). GO assignment was used to reduce complexity and classify the functions associated with the D. kotschyi unigenes. Of 165,597 unigenes, 78,188 biological process (47.2% of all transcripts), 80,356 molecular function (48.5%) and 79,769 cellular components (48.1%) terms were identified (Supplementary Fig. 3b). Based on gene ontology analysis, 91,017 transcripts could be assigned to one or more GO terms. Within the cellular component category, the three most enriched terms were “Membrane (GO:0016020)” with 28,763, “Nucleus (GO:0005634)” with 21,922, and “Cytoplasm (GO:0005737)” with 14,320 transcripts. In the biological process domain, the three most common categories were “Cellular process (GO:0009987)” with 62,195, “Metabolic process (GO:0008152)” with 55,070, and “Nucleobase-containing compound metabolic process (GO:0006139)” with 22,967 transcripts. In the molecular function domain, the three most abundant groups were “Binding (GO:0005488)” with 59,971, “Catalytic activity (GO:0003824)” with 48,681, and “Nucleic acid-binding (GO:0003676)” with 20,957 transcripts (Supplementary Fig. 4 and Supplementary excel file 2).
To further investigate the D. kotschyi transcriptome data, we searched the annotated sequences for COG/NOG classifications according to the popular EggNOG database. This analysis showed that 46,951 (28.35%) transcripts were assigned to 22 COG categories. The cluster for “Posttranslational modification, protein turnover, chaperones” (6886) represented the largest group, followed by “Signal transduction mechanisms” (6142) and “Transcription” (5994). The category of “Nuclear structure” (5) was the smallest group (Supplementary Fig. 5 and Supplementary Table 4). Also, 32,263 (19.48%) transcripts were poorly characterized in the COGs categories and grouped in the “Function unknown” classification. Complete the summary of this annotation is in Supplementary excel file 3. To find the active biological pathways of the unigenes, we also conducted a search of all transcripts against the KEGG database with the BlastKOALA tool (KEGG Orthology And Links Annotation) (https://www.kegg.jp/blastkoala/). Of 106,731 sequences, 42,715 had significant KO term matches with known enzymes in the KEGG database and were assigned to 324 KEGG pathways (Supplementary Fig. 6 and Supplementary excel file 4). The most strongly represented biological pathways were “metabolic pathways” (958), “Signal transduction” (754), and “Biosynthesis of secondary metabolites” (519).
Moreover, our analysis of the D. kotschyi transcriptome data to search the whole set of transcription factors (TFs) in PlnTFDB revealed that 6,626 assembled transcripts (6.20%) encode putative TFs. The predicted TFs can be classified into 50 TF families based on their DNA-binding domains (DBDs) according to the criteria of PlantTFdb (Supplementary excel file 5). Among predicted TF families, the C2H2-type zinc finger family showed the most abundant predicted TF (1385, 20.90%) followed by WD40-like (985, 14.87%), MYB-HB-like (626, 9.45%), bHLH (349, 5.27%) (Supplementary Fig. 7). Identifying this large set of TFs provides a rich resource for future characterization of specific TFs in various biochemical pathways in D. kotschyi.
Simple sequence repeats (SSRs) analysis
SSR markers provide many advantages over the other marker systems, and their moderate density still serves as the best co-dominant marker system for constructing framework linkage maps16. SSRs have been extensively used for various genotyping applications because they are abundant, easy to develop and detect, and highly polymorphic17. Transcriptome SSRs also exhibit high inter-specific transferability18,19,20. We counted the frequency of SSRs with different numbers of tandem repeats, and among the 165,597 examined sequences (171,215,531 bp total length), 68,044 SSRs in 46,695 sequences (28.1% of total sequences) were identified. Additionally, 14,862 sequences contained more than 1 SSR locus. Herein, the SSR motifs such as mono-, di-, tri-, tetra-, penta-, hexanucleotides were identified. The mono- 37.71% (25,661) and dinucleotide 45.96% (31,273) repeat motifs showed the highest frequency, respectively (Supplementary Fig. 8a). The dinucleotide repeat motif AG/CT was the most abundant, followed by A/T, AC/GT, AT/AT, AAG/CTT, and C/G (Supplementary Fig. 8b). The most abundant type was SSRs with ten tandem repeats, followed by six tandem repeats and five tandem repeats (Supplementary Fig. 8c and Supplementary excel file 6).
Analysis of candidate genes involved in the biosynthetic pathways of (poly) methoxylated flavones
Until now, no research aiming at the isolation of the genes involved in the methoxylated flavones biosynthesis pathway has been carried out. Therefore, for the first time, we identified the full-length genes and then combined the literature review13,14,15,21,22, and our experimental analysis, subsequently re-constructed the whole methoxylated flavones biosynthesis pathway, from phenylalanine to all methoxylated flavones which were detected in this plant5,6,23. Based on available information about flavones biosynthesis22,24, after the apigenin formation, the various type of flavonoid hydroxylases (FH) and flavonoid O-methyltransferases (FOMTs), leading the biosynthesis pathway to the production of various compounds which may be species-specific (or family). Identification of candidate genes, including FOMTs and FHs, was carried out based on structural similarity and strong protein sequence similarities of putative transcripts against related sequences from other species. In 2019, Mohammadi and co-workers considering transcriptome and metabolic data as well as literature sources could identify the genes involved in the diosgenin biosynthesis pathway, therefore, proposed the most possible pathway of diosgenin biosynthesis25. In this study, the same method and procedure were used to identify the genes involved in the biosynthesis pathway of methoxylated flavones. Finally, we investigated and identified different pathways in D. kotschyi and proposed an entire pathway for the biosynthesis of valuable anticancer and antioxidant compounds in this plant (Fig. 1).
Proposed pathway for biosynthesis of valuable compounds in D. kotschyi. Include Xanthomicrol, Calycopterin, and Penduletin (yellow box), Rosmarinic acid (purple box), and MEP pathway for terpenoids (green box). Enzyme abbreviations: PAL, phenylalanine ammonia-lyase; C4H, cinnamate 4-hydroxylase; 4CL, 4-coumarate: CoA ligase; CHS, chalcone synthase; CHI, chalcone isomerase; FNS, flavone synthase; F3H, flavanone 3-hydroxylase; F3ʹH, flavonoid 3ʹ-hydroxylase; F6H, flavonoid 6-hydroxylase; F3OMT, flavonoid 3-O-methyltransferases; F4ʹOMT, flavonoid 4ʹ-O-methyltransferases; F3ʹOMT, flavonoid 3ʹ-O-methyltransferases; F6OMT, flavonoid 6-O-methyltransferases or CRS, cirsimaritin synthase; F7OMT, flavonoid 7-O-methyltransferases; F8OMT, flavonoid 8-O-methyltransferases; HPPR, hydroxy phenylpyruvate reductase; TAT, tyrosine aminotransferase; RAS, rosmarinic acid synthase; TAL, tyrosine ammonia-lyase; DXS, 1-deoxy-D-xylulose-5-phosphate synthase; DXR, 1-deoxy-D-xylulose-5-phosphate reductoisomerase; MCT, 2-C-methyl-D-erythritol-4-(cytidyl-5-diphosphate) transferase ; CMK, 4-cytidine 5′-diphospho-2-C-methyl-D-erythritol kinase ; MCS, 2-C-methyl-D-erythritol-2,4-cyclodiphosphate synthase ; HDS, Hydroxy-2-methyl-2-(E)-butenyl 4-diphosphate synthase ; HDR, Hydroxy-2-methyl-2-(E)-butenyl 4-diphosphate reductase ; IDI, Isopentenyl diphosphate isomerase ; GPPS, Geranylgeranyl diphosphate synthase; OCS, Ocimene synthase. The proposed pathways were constructed by combining the literature review13,14,15, our experimental data, and bioinformatics data analysis. The final figure was created by Inkscape (1.0.2, https://inkscape.org).
In this study, using transcriptome sequencing data of D. kotschyi, sequences of all genes involved in the target biosynthetic pathways, including methoxylated flavones, terpenoids (MEP; methylerythritol phosphate pathway), and rosmarinic acid, were identified based on the high similarity of the putative genes with the related sequences from closely related species (Supplementary Table 5). Herein, fourteen genes including early biosynthetic genes PAL, C4H, 4CL, CHS, CHI and FNS (blue box), five genes encoding FOMTs, i.e., F3OMT, F4ʹOMT, F6OMT, F7OMT, and F8OMT (yellow box), as well as three FHs, i.e., F3H, F3ʹH and F6H (yellow box), in the biosynthesis of methoxylated flavones were identified (Fig. 1). On the other hand, TAT, RAS, and HPPR genes involved in the rosmarinic acid biosynthesis pathway were identified (purple boxes) (Fig. 1 and Supplementary Table 5). We also identified ten genes functioning in the biosynthesis of plastidic MEP pathway and main monoterpene compound in the essential oil of D. kotschyi, that is, ocimene, including DXS, DXR, MCT, CMK, MCS, HDS, HDR, IDI, GPPS, and OCS, which shown in the green box (Fig. 1). In the sequel to this study, all FOMTs and FHs, which involved in the production of valuable flavones in this plant, as well as the RAS gene in the rosmarinic acid metabolic pathway, were selected for further analysis (Supplementary Fig. 9). The in silico analysis revealed that the studied genes possess functional motifs, which are required to activity and identity of these enzymes (Fig. 2). All FOMTs have five conserved motifs in their structure that are highly similar to each homologous gene from related species26. These motifs were shown in Fig. 2, and their multiple sequence alignment is presented in Supplementary Fig. 10. Based on PredictNLS and DeepLoc results, the FOMT enzymes were predicted to be located in the cytoplasm21, and their methyltransferase activity was approved by GO:0,008,168 term. Also, the phylogenetic tree showed that the FOMTs genes of Zarrin-Giah were closely linked with those in the Lamiaceae family (Supplementary Fig. 9a). F3H contains five conserved motifs with two functional domains, including ferrous iron-binding site with HxDxnH motif (in His216, Asp218, and His274) and 2-oxoglutarate binding site with RxS motif (in Arg284 and Ser286) (Fig. 2 and Supplementary Fig. 11). This enzyme is predicted to be localized to cytoplasm27 with known oxidoreductase activity (GO:0,016,491) (Supplementary Fig. 9b). F3ʹH belongs to the CYP450 Family and possesses six motifs28 which contains the Proline reach domain (Hing) in P35PGPRPWP42, Heme binding site in F443GAGRRICAG452, and Oxygen binding site in A307GTDTT312 (Fig. 2 and Supplementary Fig. 12). This enzyme was predicted to be located in the endoplasmic reticulum29, and its monooxygenase activity was confirmed with GO:0,004,497 term (Supplementary Fig. 9c).
Represent of domains and motifs of FOMTs, F3H, F3ʹH, F6H, and RAS genes identified from D. kotschyi using Interpro database and literature data by webLogo tool.
F6H is a member of A-type CYP450s30,31 and has four common conserved family motifs (Fig. 2 and Supplementary Fig. 13), GO:0004497 represents the monooxygenase activity of this enzyme, and subcellular localization prediction indicated that this enzyme is associated with the endoplasmic reticulum membrane (Supplementary Fig. 9d). The RAS enzyme belongs to the BAHD acyltransferases family and contains two conserved regions that possess the HXXXD motif, which contains the catalytically active histidine and the DFGWG motif that is thought to be responsible for the steric position of the active site32 (Fig. 2 and Supplementary Fig. 14). This enzyme is predicted to be located in the cytoplasm, and its acyltransferase activity was approved by GO:0,016,746 term32,33. Phylogenetic analysis of the RAS enzyme also revealed that it has close relation with those in the Lamiaceae family (Supplementary Fig. 9e).
Expression analysis by qRT-PCR
In this study, we used isoform-specific primers of nine candidate unigenes (F3OMT, F4ʹOMT, F6OMT, F7OMT, F8OMT, F3H, F3ʹH, F6H, and RAS genes) associated with methoxylated flavones and rosmarinic acid biosynthetic pathway in three different tissues of D. kotschyi at flowering stage (leaves, flowers, and buds). The flower tissue showed low expression for all candidate genes, so relative expression for each gene in other tissues was compared with this tissue. Figure 3a and b show the expression of DkF3OMT, DkRAS, and DkCRS genes in the leaf (10.5, 6.8, and 6.6 respectively) and bud tissues (12.8, 9.8, and 10.4 respectively) was higher than those in the flower tissue. In addition, the expression of DkF3H and DkF3OMT genes in the leaf (1.9 for each gene) and bud tissues (4.7 and 5.1 respectively) was lower than those in the flower tissue. Relative expression analysis of candidate genes in the bud tissue to the leaf tissue showed that DkF7OMT had the higher expression and DkF3OMT had the lower expression in bud tissue (Fig. 3C).
Expression pattern of selected genes from methoxylated flavones and rosmarinic acid biosynthetic pathway using qRT-PCR in D. kotschyi. Nine genes involved in methoxylated flavones and rosmarinic acid biosynthesis, namely; F3OMT, flavonoid 3-O-methyltransferases; F4ʹOMT, flavonoid 4ʹ-O-methyltransferases; F6OMT, flavonoid 6-O-methyltransferases or CRS, cirsimaritin synthase; F7OMT, flavonoid 7-O-methyltransferases; F8OMT, flavonoid 8-O-methyltransferases; F3H, flavanone 3-hydroxylase; F3ʹH, flavonoid 3ʹ-hydroxylase; F6H, flavonoid 6-hydroxylase; RAS, rosmarinic acid synthase, were selected for qRT-PCR analysis. The results represent the means ± standard error of experiments performed in triplicate.
Flavonoids and methoxylated flavones contents
Initially, the correlation coefficients (R2) were calculated for all standards of compounds, which were about 0.99 (Supplementary Fig. 15). The chromatograms of HPLC peaks relevant to the standard of each metabolite in different tissues are shown in Fig. S16 and Fig. 4A. These results demonstrate the accumulation pattern and RT for each compound at 280 nm. As shown in Fig. 4B, the flavonoids contents varied among the different organs and tissues. In general, all of the measured metabolites had a higher amount in the mixture of flowers and buds and flower tissue. Calycopterin had a significant amount in leaf tissue (15.98 mg/g DW). After that, a high quantity of penduletin (9.30 mg/g DW) was achieved in the mixed tissue. Then apigenin (8.29 mg/g DW) and rosmarinic acid (3.34 mg/g DW) in the mixed tissue and isokaempferid (1.95 mg/g DW) in bud tissue had the highest values, respectively. The result is consistent with previous results5,8,23.
Chromatogram of HPLC peaks and accumulation of flavonoids in different organs of D. kotschyi. (A) Chromatograms of HPLC peaks corresponding to: Ros, Rosmarinic acid; Api, Apigenin; Cir, Cirsimaritin; Iso, Isokaempferid; Pen, Penduletin; Cal, Calycopterin in leaf tissue. (B) Contents of flavonoids in different organs, including flower, bud, leaf, and a mixture of tissues. The values and error bars represent the mean and standard error of three biological replicates, respectively.
