Integrated metabolomics and transcriptome analysis on flavonoid biosynthesis in safflower (Carthamus tinctorius L.) under MeJA treatment

Transcriptome sequencing and differential transcript analysis

The transcriptomes of the mixed safflower samples were sequenced. In the untreated samples,a 6.50 G clean base was obtained (44,185,480 raw reads and 43,921,238 clean reads), and Q30 was 92.40%. In the MeJA-treated samples, a 6.17G clean base was obtained (42,087,290 raw reads and 41,720,790 clean reads), and Q30 was 91.89%. The different expressed genes were analysed with DESeq2. The total number of differentially expressed genes was 31,822 (20,741 upregulated genes and 11,081 downregulated genes). (Fig. 2).

Open in a separate window

Fig. 2

Volcano plot of differential expression genes. The abscissa represents the change of gene expression multiple (log₂ Fold Change), and the ordinate represents the significant level of differentially expressed genes (−log10 False Discovery Rate). The expression of red gene was up-regulated, that of green gene was down-regulated, and that of black gene was not significantly different

The differential flavonoid metabolite transcripts from each comparison group were screened by using the KEGG database. The experiments focused on annotations in the phenylpropane metabolic pathway, including flavonoid biosynthesis (ko00941), anthocyanin biosynthesis (ko00942), isoflavonoid biosynthesis (ko00943), and flavone and flavonol biosynthesis (ko00944). The results showed that there were 61 significantly different flavonoid biosynthesis transcripts between MeJA-treated and untreated materials (43 in ko00941, 8 in ko00942, 5 in ko00943, and 5 ko00944). Most of transcripts were involved in flavonoid biosynthesis (ko00941). The details could be viewed in Supplementary Table S3.

Integrated analysis of the transcriptome and metabolome on the KEGG pathway

The metabolic components were mapped onto the pathway of flavonoid metabolism (including ko00941, ko00942, ko00943 and ko00944) by combining metabolic components that had been detected by UHPLC-ESI-MS/MS, thus constructing the metabolic pathway map of integrative flavonoid biosynthesis in safflower. (Fig. 3) At the same time, the differentially annotated metabolites, together with the differentially annotated genes, were indicated on the integrated metabolic map. (Fig. 3) From the integrated metabolic map, it was shown that MeJA might upregulate the expression of the upstream genes in the flavonoid biosynthesis pathway (such as CHSs, CHIs, and HCTs) and might downregulate the expression of downstream genes (such as F3Ms, ANRs, and ANSs), thus promoting the biosynthesis of quinone chalcones, such as HSYA.

Open in a separate window

Fig. 3

The metabolic pathway map of integrative flavonoid biosynthesis in safflower. The metabolic components were mapped into the pathway of flavonoid metabolism (including flavonoid biosynthesis (ko00941); Anthocyanin biosynthesis (ko00942); Isoflavonoid biosynthesis (ko00943); Flavone and flavonol biosynthesis (ko00944). The up arrow indicates an increase in content and the down arrow indicates a decrease in content. The reported cloned genes in safflower were tagged in the map

Expression analysis of differentially expressed flavonoid biosynthesis genes by real-time PCR

Real-time PCR was used to validate the transcriptome data. Ten differentially expressed genes were selected, including TRINITY_DN28401_c0_g1 (annotated as CHI), TRINITY_DN36537_c1_g1 (annotated as HCT), TRINITY_DN37574_c1_g2 (annotated as HCT), TRINITY_DN39063_c0_g1 (annotated as FLS), TRINITY_DN43344_c1_g1 (annotated as CHI), TRINITY_DN41734_c0_g5 (annotated as CHS), TRINITY_DN42700_c3_g2 (annotated as F3H), TRINITY_DN43120_c5_g1 (annotated as ANR), TRINITY_DN44094_c0_g1 (annotated as ANS), TRINITY_DN44565_c1_g1 (annotated as F3M). The primers used for this experiment could be found in Supplementary Table 4. Among them, HCT contains several transferase enzymes, which include anthranilate N-hydroxycinnamoyl/benzoyltransferase that catalyses the first committed reaction of phytoalexin biosynthesis. F3H can catalyse synthesis N-terminal flavanone 3-hydroxylase. F3M can catalyse flavonoid, NADPH, H⁺, and O₂ into 3′-hydroxyflavonoid and NADP⁺, and H₂O, which functions like F3H. FLS catalyzes the formation of flavonols from dihydroflavonols. CHS catalyses the reaction of one molecule of 4-coumaroyl-CoA and three molecules of malonyl-CoA to form tetrahydroxychalcone. Chalcone isomerase converts tetrahydroxychalcone into naringenin. ANR and ANS are all involved in the synthesis of anthocyanins. The results showed that the expressions of 6 genes, including TRINITY_DN28401_c0_g1 (CHI), TRINITY_DN36537_c1_g1 (HCT), TRINITY_DN37574_c1_g2 (HCT), TRINITY_DN39063_c0_g1 (FLS), TRINITY_DN43344_c1_g1 (CHI) and TRINITY_DN41734_c0_g5 (CHS), were significantly upregulated, and 4 genes, including TRINITY_DN42700_c3_g2 (F3H), TRINITY_DN43120_c5_g1 (ANR), TRINITY_DN44094_c0_g1 (ANS) and TRINITY_DN44565_c1_g1 (F3M), were significantly downregulated. (Fig. 4) From the integrated metabolic map, it could be seen that the expressions of upstream genes in the phenylpropane metabolic pathway were increased with MeJA treatment. And the results of gene expression also suggested that these genes might be involved in metabolic pathways.

Open in a separate window

Fig. 4

Real-time PCR expression of 10 genes from the flavonoid biosynthesis pathway. 6 genes significantly upregulated and 4 genes significantly down-regulated are included. Gene1 represents TRINITY_DN28401_c0_g1 (annotated as CHI), Gene2 represents TRINITY_DN36537_c1_g1 (annotated as HCT), Gene3 represents TRINITY_DN37574_c1_g2 (annotated as HCT), Gene4 represents TRINITY_DN39063_c0_g1 (annotated as FLS), Gene5 represents TRINITY_DN43344_c1_g1 (annotated as CHI), Gene6 represents TRINITY_DN41734_c0_g5 (annotated as CHS), Gene7 represents TRINITY_DN42700_c3_g2 (annotated as F3H), Gene8 represents TRINITY_DN43120_c5_g1 (annotated as ANR), Gene9 represents TRINITY_DN44094_c0_g1 (annotated as ANS), Gene10 represents TRINITY_DN44565_c1_g1 (annotated as F3M)

Sample preparation and extraction for metabolomic analysis

The freeze-dried flowers were crushed using a mixer mill (MM 400, Retsch) with a zirconia bead for 1.5 min at 30 Hz. One hundred milligrams of powder was weighed and extracted overnight at 4 °C with 1.0 mL 70% aqueous methanol. Following centrifugation at 10000 g for 10 min, the extracts were absorbed (CNWBOND Carbon-GCB SPE Cartridge, 250 mg, 3 mL; ANPEL, Shanghai, China, www.anpel.com.cn/cnw) and filtered (SCAA-104, 0.22 μm pore size; ANPEL, before LC-MS analysis.

UHPLC conditions

The sample extracts were analysed using an LC-ESI-MS/MS system (HPLC, Shim-pack UFLC SHIMADZU CBM30A system; MS, Applied Biosystems 6500 Q TRAP). The analytical conditions were as follows: UHPLC column, Waters ACQUITY UHPLC HSS T3 C18 (1.8 μm, 2.1 mm × 100 mm); solvent system, A water (0.04% acetic acid): B acetonitrile (0.04% acetic acid); gradient program, 0% B at 0 min, 95% B at 11.0 min, 95% B at 12.0 min,5% B at 12.1 min, and 5% B at 15.0 min; flow rate, 0.40 mL/min; temperature, 40 °C; and injection volume, 2 μL. The UHPLC effluent was connected to an ESI-triple quadrupole-linear ion trap (Q TRAP)-MS.

ESI-q trap-MS/MS

Mass spectrometry followed the method of Chen et al. [19]. LIT and triple quadrupole (QQQ) scans were acquired on a triple quadrupole-linear ion trap mass spectrometer (Q TRAP), API 6500 QTRAP LC/MS/MS system equipped with an ESI Turbo ion-spray interface, operating in positive ion mode and controlled by Analyst 1.6.3 software (AB Sciex, Waltham, MA, USA). The ESI source operation parameters were as follows: ion source, turbo spray; source temperature, 500 °C; ion spray voltage (IS), 5500 V; the ion source gas I (GSI), gas II (GSII), and curtain gas (CUR) were set at 55, 60, and 25.0 psi, respectively; and the collision gas (CAD) was high. Instrument tuning and mass calibration were performed with 10 and 100 μmol/L polypropylene glycol solutions in QQQ modes. QQQ scans were acquired as MRM experiments with collision gas (nitrogen) set to 5 psi. DP and CE for individual MRM transitions was performed with further DP and CE optimization. A specific set of MRM transitions was monitored for each time period according to the metabolites eluted within that period.

Qualitative and quantitative analysis of metabolites

The flavonoid identification and quantification in our study was made according to a method of scheduled multiple reaction monitoring (MRM), which has been previously described [19]. With this methods, Chen made a genome-wide association analyses between genetic and biochemical, and reported 840 metabolites, from which 277 were identified or annotated [20]. Mainly based on this research, a database was built. All the metabolites in our study were detected or annotated by this database. The database and methods have been wildly used to detected or annotate the metabolites, such as the research of Wang [36] and Liu [37]. This methods is also called as wild-targetd analysis. As HSYA was not in the database, we use its standard sample (No.: MUST-15072815, Chengdu Manster Biotechnology Co., LTD, China) to build a library and add it to the self-built database. The metabolites of the samples were qualitatively and quantitatively analysed by mass spectrometry. The characteristic ions of each substance were filtered by the triple quadrupole, and the signal strength of the characteristic ions was obtained in the detector. Chromatographic peaks were analyzed with Multi Quant software 3.0.3. The integrated area peak of each compound was used for the PCA and OPLS-DA analysis..

Differential metabolite analysis

Principal component analysis (PCA) can effectively extract main variance information and was used in many other research [29, 38]. In addition, orthogonal partial least squares-discriminant analysis (OPLS-DA) was used for variables with less correlation. The experiment had biological duplication; thus, the fold change and VIP value of the OPLS-DA model were combined to screen differential metabolites. PCA was analyzed by R software built-in functions (www.r-project.org/). The parameter: scale = True. After conversion of the original data by log2, the data was centralized (Mean Centering) and analyzed by OPLSR. Anal of Metabo Analyst in R software. The main steps of our study can be referred to a previouly research [39]. Screening criteria: 1. Metabolites must exhibit fold changes≥2 and fold changes≤0.5. If a metabolite in the experimental groupwas more than 2 times or less than 0.5 times that of the metabolite in the control group, the difference was considered significant. 2. If there is biological duplication in the sample grouping, the metabolites of VIP ≥1 are selected on the basis of the above. The VIP value indicates the influence intensity of the difference between groups of corresponding metabolites in the classification and discrimination of each group of samples in the model.

RNA sequencing and annotation

RNA isolation and purification and cDNA library construction and sequencing were performed as previously described [40]. All tissues were ground on dry ice, and total RNA was prepared by using TRIzol reagent (Invitrogen, CA, USA). To remove DNA, an aliquot of total RNA was treated with DNase (Takara, Dalian, China). RNA quantity and quality were determined by using a NanoDrop 2000 spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA) and an Agilent 2100 Bioanalyser (Agilent Technologies, CA, USA), respectively. mRNA was isolated from total RNA using magnetic beads with oligo (dT) primer; cDNA was synthesized using a cDNA synthesis kit (TaKaRa, Dalian, China) and linking the sequencing adapter to both ends. The library preparations were sequenced on an Illumina HiSeq 4000 platform, and the unigene sequences obtained from our laboratory transcriptome database by RSEM (RNAseq by Expectation-Maximization) software were integrated for annotation. The whole set of transcript data can be found in the National Center for Biotechnology Information (NCBI) SRA database (SRR10011980 and SRR10011979).

Screening of differential genes

DESeq2 [41, 42] was used for differential expression analysis among sample groups. DESeq2 requires the input of unregulated read counting data of genes rather than RPKM (Reads Per Kilobase Million), FPKM (Fragments Per Kilobase Million) and other standardized data. Read counts of genes are the expected_count outputs calculated using RSEM (RNAseq by Expectation-Maximization). The expected count is generally lower than the read number. After discrepancy analysis, multiple hypothesis tests are needed to correct the hypothesis test probability (P value) with the Benjamini-Hochberg procedure to obtain the false discovery rate (FDR) when |log₂fold change| ≥ 1and FDR < 0.05.

Real-time PCR analysis

Real-time PCR analysis was performed as previously described [14]. Total RNA was isolated using an RNA extraction kit (Invitrogen, CA, USA), and reverse transcription was carried out using the prime script reagent kit (Takara, Dalian, China). The primers used to amplify the screened genes (CHS, CHI, HCT, etc.) by real-time PCR were designed by Primer 5.0, and parts of the safflower 28S coding region were used as an internal reference gene. The primer details are listed in Supplementary Table 4. All of the primers were tested for their specificity by agarose gel electrophoresis. RT-PCR was performed using an SYBR prime script RT-PCR kit (Takara) with three replicates, and the cycling conditions were set according to the manual. The Bio-Rad CFX96 real-time PCR detection system (Hercules, CA, USA) was used in our experiment.

We acknowledge Dr. Hu Zeng and Yuanyuan Qin (Metware Biotechnology Co., Ltd. Wuhan, China) for the data analysis, American Journal Experts (AJE) for the language editing (verification code: B631-386C-E812-C6A9-351B), and Renchuan Yao for providing us the plant materials.