Integrated transcriptomic and metabolomic data reveal the flavonoid biosynthesis metabolic pathway in Perilla frutescens (L.) leaves

Perilla frutescens (L.) is an important medicinal and edible plant in China with nutritional and medical uses. The extract from leaves of Perilla frutescens contains flavonoids and volatile oils, which are mainly used in traditional Chinese medicine. In this study, we analyzed the transcriptomic and metabolomic data of the leaves of two Perilla frutescens varieties: JIZI 1 and JIZI 2. A total of 9277 differentially expressed genes and 223 flavonoid metabolites were identified in these varieties. Chrysoeriol, apigenin, malvidin, cyanidin, kaempferol, and their derivatives were abundant in the leaves of Perilla frutescens, which were more than 70% of total flavonoid contents. A total of 77 unigenes encoding 15 enzymes were identified as candidate genes involved in flavonoid biosynthesis in the leaves of Perilla frutescens. High expression of the CHS gene enhances the accumulation of flavonoids in the leaves of Perilla frutescens. Our results provide valuable information on the flavonoid metabolites and candidate genes involved in the flavonoid biosynthesis pathways in the leaves of Perilla frutescens.

Perilla frutescens (L.), which is a self-compatible annual herb, belongs to the family Lamiaceae. This species has been widely cultivated in China, Japan, and Korea for centuries. Perilla frutescens is an important medicinal and edible plant in China with medical and nutritional uses 1 . Its leaves can be utilized as a transitional medicinal herb, as a vegetable, and as a spice, and its seeds can be processed into foods and nutritional edible oils 2 . The extract from leaves of Perilla frutescens is composed of many chemical components containing flavonoids, polysaccharides, amino acids, and trace elements, which are mainly used in traditional Chinese medicine 3 . Flavonoids serve specific functions to the plant under specific developmental or biotic/abiotic conditions and are the first line of defense against ultraviolet rays and pathogens 4 . Previous studies reported that flavonoids have many biological functions such as anti-inflammatory, anti-oxidative, anti-diabetic, and anti-hypertensive activities [5][6][7][8] .
The biosynthesis of flavonoids is regulated by a series of enzymes related to the phenylpropanoid and flavonoid pathways 9 . Among these enzymes, phenylalanine ammonia-lyase (PAL), cinnamic acid 4-hydroxylase (C4H), 4-coumarate coenzyme A ligase (4CL), and acetyl-CoA carboxylase (ACC) catalyze the biosynthesis of phenylpropanoids. Chalcone synthase (CHS), chalcone isomerase (CHI), favonol synthase (FLS), flavanone 3-hydroxylase (F3H), flavonoid 3′-hydroxylase (F3′H), flavonoid 3′5′-hydroxylase (F3′5′H), dihydroflavonol 4-reductase (DFR), anthocyanidin synthase/leucocyanidin dioxygenase (ANS/LDOX), and flavonoid 3-O-glucosyltransferase (UFGT) are key enzymes that control the biosynthesis of flavonoids. The above-mentioned genes involved in the biosynthesis of flavonoids have been reported in many plant species such as Arabidopsis and Ginkgo biloba 10,11 . In addition, the expression levels of genes and transcription factors also play key roles in regulating flavonoid biosynthesis 12 . For example, a high expression of CHS and F3H genes increases the accumulation of flavonoids Flavonoid metabolites in the leaves of Perilla frutescens. To compare the differential flavonoid metabolites between JIZI 1 and JIZI 2, the data obtained from UPLC/ESI-Q TRAP-MS/MS were analyzed. In this work, a total of 223 different flavonoid metabolites were identified in the leaves of Perilla frutescens (Table S1). The heatmap of metabolites was drawn by R software after unit variance scaling (UV), and hierarchical cluster analysis (HCA) was performed on the accumulation pattern of metabolites among different samples (Fig. 2a). The 223 flavonoid metabolites were classified into 6 categories, including 109 flavones, 33 flavonols, 27 flavonoids, 22 anthocyanins, 20 flavanones, and 12 isoflavones (Fig. 2b). Among the 223 flavonoid metabolites, chrysoeriol, apigenin, malvidin, cyanidin, kaempferol, and their derivatives were abundant in the leaves of Perilla frutescens, which were more than 70% of total flavonoid contents. Using the identification criterion of the absolute Log 2 FC ≥ 1 and VIP value ≥ 1, 57 (25.6%) flavonoid contents were significantly different among the

Differentially expressed genes in the leaves of Perilla frutescens.
To understand the molecular basis of flavonoid biosynthesis metabolic pathway, transcriptomes were analyzed to identify differentially expressed genes in the leaves. A total of 54.75 million clean reads were produced from the leaves of Perilla frutescens. These clean reads were further assembled into 60,458 unigenes with a mean length of 1194.19 bp using Trinity software. All unigenes were searched in the Nr, GO, COG, Swiss-Prot, egg-NOG, Pfam, KOG, and KEGG databases with the BLASTX (E-value < 1.0E-5) program for functional annotations. Among of 60,458 assembled unigenes, a total of 42,580 (70.43%) unigenes were annotated functionally based on the above-mentioned databases ( Fig. 3a and Table S2). With the filter criteria of |Log 2 FC| ≥ 1 and FDR < 0.05, 9277 differentially expressed genes (DEGs) were identified in the leaves of Perilla frutescens, including 1642 upregulated genes and 7635 downregulated genes (Fig. 3b). In GO enrichment analysis of DEGs, 3711 DEGs out of 9277 were involved in three major GO categories, i.e., biological process, cellular component and molecular function (Fig. 3c), and 12 phenylpropanoid biosynthetic process, 29 flavonoid biosynthetic processes, 15 flavonoid glucuronidation processes, and 9 anthocyanin-containing compound biosynthetic processes were identified in the biological process category.
To verify the credibility of transcriptome information, we further selected 15 DEGs to validate the sequencing results. The qRT-PCR results showed that 11 genes showed higher expression levels, 4 genes were lower expression in JIZI 2 than in JIZI 1, our qRT-PCR results were consistent with those obtained with the RNA-Seq method (Fig. 5).     Table S4). Among these transcription factors, the most abundant were MYBs (28) and AP2/ERFs (17), followed by WRKYs (14), bHLH (9), MADs (9), WD40s (6), and NACs (3). Interestingly, WD40s and NACs were all upregulated. These transcription factors might contribute to flavonoid metabolites in the leaves of Perilla frutescens.

Discussion
UPLC/ESI-Q TRAP-MS/MS is popular in the field of identification and analysis of plant metabolites, which has the advantages of high sensitivity and throughput, fast separation, and wide coverage. So far, this technology is widely applied to analyze the metabolites in tomato, strawberry and asparaguses [16][17][18] . In recent years, metabolomics integrated with transcriptomics is widely used to investigate the biosynthesis of metabolites to reveal the biosynthesis pathways of metabolites in plants 19 . Perilla frutescens (L.) has been used for centuries as a traditional medicinal herb in China. Increasing attention has been given to the anti-allergic, antioxidant, anti-inflammatory, and anti-tumor activities of perilla plants [20][21][22][23][24] . It was reported that the flavonoids in Perilla frutescens leaves have hypolipidemic and antioxidant effects 25 . In particular, luteolin and apigenin can improve several hyper-monoaminergic neuropsychological disorders as monoamine transporter activators 26 . In our study, the total flavonoid content of JIZI 2 was significantly higher than that of JIZI 1 (Fig. 1). In order to elucidate flavonoid biosynthesis in the leaves of Perilla frutescens, metabolomic and transcriptomic data were collected for JIZI 2 and JIZI 1 leaves. A total of 223 flavonoids were identified in the leaves of Perilla frutescens by UPLC/ESI-Q TRAP-MS/MS. Of these flavonoids, 57 were significantly different in JIZI 2 vs. JIZI 1, including 46 upregulated flavonoids and 11 downregulated flavonoids (Fig. 2). Chrysoeriol is a flavonoid metabolite with anti-inflammatory, anti-tumor, and cardioprotective agent effects. Choi et al. reported that chrysoeriol isolated from the leaves of Digitalis purpurea could inhibit the induction of nitric oxide synthase by blocking AP-1 activation 27 . Yang et al. reported that chrysoeriol had anti-tumor activity against human multiple myeloma cells 28 . Chrysoeriol could potentially serve as a novel cardioprotective agent against doxorubicin-induced cardiotoxicity 29 . In our study, we found the contents of chrysoeriol and its glycosides were highest in the leaves of Perilla frutescens, which were almost 50% of total flavonoid contents. In addition, the apigenin and its derivatives were also abundant in the leaves of Perilla frutescens, which are recognized as bioactive flavonoids possessing anti-oxidant, anti-inflammatory, and anti-cancer activities 30 (Table S1).
Combinatorial analysis of transcriptomic and metabolomic are useful for understanding of molecular regulation of secondary metabolites. Fukushima et al. reported that high expression levels of genes encoding F3′H, DFR, and ANS leaded to accumulate more anthocyanins in red perilla than in green perilla 31 . In this study, transcriptome analysis of the leaves of Perilla frutescens identified unigenes involved in the flavonoid biosynthesis process and showed differentially expressed genes. In the flavonoid biosynthesis pathway, the expression levels of genes encoding PAL, C4H, CHS, F3H, F3′H, DFR, ANS, UFGT, and UGT75C1 were higher in JIZI 2 than in JIZI 1 (Fig. 4) which resulted in a higher flavonoid content in JIZI 2. CHS is a key enzyme of the flavonoid pathway that plays an important role in the phenylpropanoid pathway and in flavonoid biosynthesis. In this study, we detected 4 differentially expressed CHS genes that were all upregulated (Table 1), and metabolic analysis showed that the abundance of naringenin chalcone, which is primarily catalyzed by CHS, was also higher in JIZI 2 than in JIZI 1 (Table S1) The biosynthesis of flavonoids is mostly regulated by transcription factors at the transcription level. The MYB-bHLH-WD40 complex, Zinc finger, MADS-box, NACs, and WRKY proteins have been identified to play a role in flavonoid biosynthesis in plants [34][35][36] . The R2R3-MYB transcription factors regulate the accumulation of flavonol in different parts of Arabidopsis seedlings 37 . MYB12 controls flavonol biosynthesis by controlling the expression of CHS, CHI, F3H and FLS 38 . In addition, transcription factor WRKY23 controls flavonol biosynthesis during Arabidopsis root development 39 . The expression levels of bHLH transcription factors were significant upregulation in red perilla than in green perilla 31 . The bHLH factor Myc-F3G which regulated the expression of anthocyanin genes was detected specifically in red perilla but not in green perilla 40 . In our study, we analyzed the transcriptome data and found that 86 important transcription factors including MYB, AP2/ERFs, WRKY, bHLH, MADS, WD40, and NACs showed significantly different expression levels (Table S4). These differentially expressed transcription factors might be candidate regulators of flavonoid biosynthesis in the leaves of Perilla frutescens.
In summary, we found the total flavonoid content was different in the Perilla frutescens varieties "JIZI 1" and "JIZI 2". We used metabolomics and transcriptomics to data reveal the flavonoid biosynthesis metabolic pathway.

Measurement of total flavonoid content. Approximately 2.5 g powder of perilla frutescens leaves was
used to measure the total flavonoid content by the Aluminum nitrate colorimetric method. In brief, 0.5 mL crude extract of leaves was mixed with 5.5 mL of 50% ethanol and l mL of 5% NaNO 2 solution. Then, 1 mL of 10% Al(N0 3 ) 3 solution was added after 6 min of incubation, and the mixture was incubated for another 6 min. Subsequently, 10 mL of 4% NaOH solution and 7 mL of 50% ethanol were added, and the final volume of the mixture solution was 25 mL. The mixture solution stood for 15 min, and then the absorbance was measured at a wavelength of 506 nm by an ultraviolet spectrophotometer (V-5100B, METASH, Shanghai, China). Rutin was used as a standard solution to prepare a calibration curve, and the results were expressed as rutin equivalent on a dry weight basis 41 .
Measurement of relative anthocyanin content. Perilla frutescens leaves (0.1 g) were ground with 1 ml of methanol (0.1% HCl) and were washed twice into 10 mL centrifuge tubes. The final volume of the samples was 5 mL (including methanol (0.1% HCl)). The tissue homogenates were oscillated for 30 s and centrifuged at 4 °C and 12,000 g for 10 min, and the absorbance of the supernatants was measured at a wavelength of 530 nm using an ultraviolet spectrophotometer (V-5100B, METASH). The relative anthocyanin content was calculated using the following formula: Identification and quantitative analysis of metabolites. Qualitative analysis of the metabolite data was performed base on the Metware Database (Metware Biotechnology Co., Ltd. Wuhan, China). A home-made software reference to MetDNA of Zhu lab and Masterview from AB Sciex were used to give the matching score of metabolites with the Metware Database. Most compounds of the database were standards. The metabolites for which no standards were available, peaks in the MS2T library, mainly the peaks that showed similar fragmentation patterns with the identified metabolites, were used to query the MS 2 spectral data taken from the literature or to search the public metabolite databases (such as MassBank, KNAPSAcK, HMDB, MoToDB, and METLIN) 45 . The quantitative analysis of metabolites used multiple reaction monitorin 46 . In the MRM mode, the precursor ions were fragmented after the induced ionization of the collision chamber to form many fragment ions, then selected a characteristic fragment ion by triple quadrupole filtering to eliminate the interference of non-target ions. After obtaining the metabolite mass spectrometry data of different samples, the peak area of all substance mass peaks was integrated, and the peaks of the same metabolite in different samples were integrated and corrected. An unsupervised PCA (principal component analysis), HCA (hierarchical cluster analysis), and OPLS-DA (partial least squares-discriminant analysis) were performed within R 47 . The unsupervised PCA was performed by statistics function prcomp within R (www.r-proje ct.org), and the data was unit variance scaled before unsupervised PCA. The HCA results of samples and metabolites were presented as heatmaps with dendrograms,

RNA sequencing (RNA-seq) data analysis and annotation.
To acquire high-quality reads, the raw reads in fastq format were processed through in-house Perl scripts. Clean reads were obtained from raw data by removing adaptor sequences, low-quality reads, and reads containing ploy-N. All downstream analyses were based on clean data with high quality. Transcriptome assembly was accomplished using Trinity software (ver-sion2.5.1) 49 50 . Analysis of the differentially expressed genes of the two groups was performed with the DESeq R package (1.10.1). DESeq provides statistical routines for determining differentially expressed genes using a model based on the negative binomial distribution. The results of all statistical tests were corrected by multiple testing using the Benjamini and Hochberg false discovery rate. Genes were determined to be significantly differentially expressed at an adjusted P value < 0.05 according to DESeq. GO enrichment analysis of the differentially expressed genes was implemented by the topGO R package based on the Kolmogorov-Smirnov test. Pathway analysis elucidated significant pathways of differentially expressed genes according to the KEGG database (https ://www.genom e.jp/kegg) [51][52][53] . We tested the statistical enrichment of differentially expressed genes in KEGG pathways using KOBAS software 54 . The whole set of raw data can be found in the national center for biotechnology information (NCBI) SRA database (accession number SRP270356).

QRT-PCR expression analysis of genes involved in anthocyanin biosynthesis. Total RNA of
Perilla frutescens leaves was reverse-transcribed according to the Quantscript Reverse Transcriptase Kit (TIAN-GEN Biotech, Beijing, China). cDNA was used as a template to measure gene expression. The specific primers involved in the anthocyanin biosynthesis genes and the Perilla frutescens actin gene (internal control) are listed in Table S5. Quantitative real-time polymerase chain reaction (qRT-PCR) was conducted by a real-time PCR ABI Prism 7500 system (software for 7500 and 7500 Fast Real-Time PCR Systems, V2.0.1, Foster City, CA, USA) using SYBR Premix Ex Taq II (TaKaRa Code No. RR820A, https ://www.takar a biomed.com.cn). The comparative CT method (2 −ΔΔCT method) was used to quantify gene expression 55 .
Statistical analysis. Statistical analysis was performed using Excel 2010 software (Microsoft Office, USA).
Data are presented as means ± standard deviations (SD). The levels of statistical significance were analyzed by the least significant difference (p < 0.05).