Transcriptome analysis of flavonoid biosynthesis in safflower flowers grown under different light intensities

Plant materials

The safflower flowers used in this study were obtained from the phytotron at the Chengdu University of Traditional Chinese Medicine in Chengdu, China. Prior to the experiment, the plants were grown at 25 °C during the day and at 18 °C at night; the relative humidity ranged from 50% to 60%. Safflower plants with the same growth status were selected and grown in a peat: vermiculite: perlite (3:1:1) mixture in a phytotron room. The 8 months old plants were divided into three groups (30 independent samples in each group) with equal size buds. Plants were illuminated with light-emitting diodes (LEDs) among three groups. The high light (HL) group was grown under 40,000 Lux. The second group, middle light (ML), was grown under 20,000 Lux. The third group, low light (LL), was grown under 10,000 Lux. Light intensity settings referred to light compensation point (LCP) and light saturation point (LSP) of safflower (Guan et al., 2007). Three different treatments lasted from budding to blooming. Blooming flowers were collected per group from the same location on the plants. These samples were immediately frozen at −80 °C. Three biological replicates were used in each treatment.

RNA extraction

Flowers were ground in liquid nitrogen, and the total RNA was isolated by TRIzol reagent (Invitrogen, Carlsbad, CA, USA). To remove DNA, an aliquot of total RNA was treated with DNase (Takara, Dalian, China) by using the standard protocol described by the manufacturer. The purity of the RNA samples according to the A260/A280 ratio was determined with a NanoDrop; the A260/A280 ratios of all samples were in the approximate range of 1.9–2.1. The integrity of the RNA samples was assessed with an Agilent 2100 Bioanalyzer, and samples with no sign of degradation were selected for further analysis.

Library construction and sequencing

A total of 9 samples, three per treatment, were randomly selected to constructed libraries. The first step in the workflow involved purifying the poly-A containing mRNA molecules using poly-T oligo-attached magnetic beads. Following purification, the mRNA was fragmented into small pieces through divalent cations under elevated temperature. The cleaved RNA fragments were copied into first strand cDNA using reverse transcriptase and random primers. This was followed by second strand cDNA synthesis employing DNA Polymerase I and RNase H. A single ‘A’ base was then added to these cDNA fragments, and an adapter was subsequently ligated to the cDNA fragments. The products were then purified and enriched with PCR amplification. Then, the PCR yield was quantified by Qubit, and samples were pooled together to make a single-strand circular DNA (ssDNA circle), which was the final library. DNA nanoballs (DNBs) were generated with the circular ssDNA by rolling circle replication (RCR) to enlarge the fluorescent signals during the sequencing process. The DNBs were loaded into the patterned nanoarrays, and pair-end reads of 100 bp were read with the BGISEQ-500 platform for the following data analysis. For this step, the BGISEQ-500 platform combined the DNA nanoball-based nanoarrays and stepwise sequencing employing the combinational probe-anchor synthesis sequencing method. Raw reads containing a low quality, adaptor-polluted and a high content of unknown bases (N) were removed before downstream analyses. Reads cleaning and filtering were performed by BGI internal software.

De novo assembly and functional annotation

Because published genome of safflower was not integrated, this project was done without a reference genome. Clean reads in each sample were assembled independently to obtain a reference sequence for subsequent analysis. De novo transcriptome assembly was completed by Trinity v2.0.6 (Grabherr et al., 2011) with clean reads, and then TGICL v2.0.6 (Pertea et al., 2003) was used on cluster transcripts to remove redundant transcripts and obtain Unigenes. The Unigenes were divided into two types: one type was a cluster in which the prefix was CL with the cluster ID behind it. In one cluster, there are several Unigenes which similarity between them is more than 70%. Another type was a singleton in which the prefix was Unigene. The N50, N70 and N90 length were used to determine the assembly continuity. The higher the length of N50, N70 and N90, the better the assembling continuity was. The raw data were uploaded to the NCBI and Figshare.

Unigene sequences were aligned to the following databases with BLAST (v2.2.23): NT (https://www.ncbi.nlm.nih.gov/nucleotide), NR (https://www.ncbi.nlm.nih.gov/refseq/about/nonredundantproteins/), COG (http://www.ncbi.nlm.nih.gov/COG), KEGG (http://www.genome.jp/kegg) and SwissProt (http://www.ebi.ac.uk/swissprot/). Based on the NR annotation, GO functional annotation (http://geneontology.org) was obtained through the Blast2GO program (v2.5.0). InterPro annotations (http://www.ebi.ac.uk/interpro) was obtained by InterProScan 5.

Analysis of the differentially expressed genes

Clean reads were mapped to the Unigenes by Bowtie2 v2.2.5 (Langmead & Salzberg, 2012) and then the gene expression levels were calculated through RSEM v1.2.12 (Li & Dewey, 2011). Differential expression analysis was performed by DEseq2 as described by Love, Huber & Anders (2014). DEseq2 is a package implemented in R based on the negative binomial distribution. A fold change ≥2.00 and an adjusted P value ≤ 0.05 were set as the thresholds to establish significantly differential expression. Classification and functional enrichment analysis including Gene Ontology (GO) and KEGG were performed to identify which GO terms or metabolic pathways were significantly enriched in the Differentially expressed genes (DEGs).

Real-time PCR expression analysis of DEGs

Because of the pharmacological effects of safflower flavonoid, we focused on DEGs related to flavonoid biosynthesis in this study. Real-time PCR (RT-PCR) analysis was employed to confirm the results of the RNA sequencing, and expressions of selected DEGs were measured in three different treatments. Each sample replicated 3 times. Specific primers were designed by Primer Premier 5 software (Table S1). Gene expression under different conditions was measured with the CFX96™ Real-time System (Bio-Rad, Hercules, CA, USA) using SYBR Premix Ex Taq II (TaKaRa, Japan). The 25S rRNA gene obtained from Carthamus tinctorius L. was used as the reference gene to identify differences in each cDNA template. The RT-PCR cycling conditions were as follows: 95 °C for 3 min, followed by 40 cycles of 95 °C for 10 s and 61 °C for 30 s. The 2^−ΔΔCt method (Livak & Schmittgen, 2001) was employed to analyze relative gene expressions. ANOVA with a post-hoc Tukey-Test in SPSS (version 20) was used for data analysis. Three HCT genes, two FLS genes and two ANS genes related to flavonoid biosynthesis, namely, CtHCT1, CtHCT2, CtHCT3, CtFLS1, CtFLS2, CtANS1 and CtANS2, were successfully amplified and chosen to certify the results.

Flavonoid quantification

Total flavonoid content was determined by the colorimetric method as previously described (Yu et al., 2010) with some modifications. Powdered flower tissue (0.2 g) was dissolved in 10 ml of 50% (v/v) methanol and subjected to ultrasonic extraction for 90 min. Scanning in the wavelength range of 200–600 nm, the optimum absorption peak wavelength was 255 nm. Therefore, the absorbance of the solution was measured at 255 nm by SpectraMax iD3 (Molecular, USA). Each sample replicated 3 times. ANOVA with a post-hoc Tukey-Test in SPSS was used for data analysis.

To determine the HSYA content, 0.4 g of powdered flower tissue was dissolved in 50 ml of 25% (v/v) methanol and then subjected to ultrasonic extraction for 40 min. The samples were then diluted with methanol, filtered through a 0.45 μm filter membrane and analyzed by high performance liquid chromatography (HPLC) (Agilent 1200, Santa Clara, CA, USA). The HSYA content was analyzed with an Agilent C18 chromatographic column (4.6 mm × 250 mm, 5 μm) in conjunction with a mobile phase that consisted of methanol, acetonitrile and 0.7% (v/v) phosphoric acid (26:2:72) for elution. The flow rate was 0.8 ml min⁻¹, the injection volume was 10 μl and the UV detector was set such that λ = 407 nm. To determine the kaempferol content, 1.0 g of powdered flower tissue was dissolved in 25 ml of methanol and then subjected to heat reflux extraction at 95 °C for 30 min followed by hydrochloric acid hydrolysis for 30 min. The sample solutions were diluted with methanol, filtered through a 0.45 μm filter membrane and analyzed by HPLC. The kaempferol content was analyzed with an Agilent C18 chromatographic column in conjunction with a mobile phase consisting of methanol and 0.4% (v/v) phosphoric acid for elution. The flow rate was 1.0 ml min⁻¹, the injection volume was 20 μl and the UV detector was set such that λ = 360 nm and 254 nm. Each sample replicated 3 times. ANOVA with a post-hoc Tukey-Test in SPSS was used for data analysis.

Sequencing and functional annotation

By the use of the BGISEQ-500 high-throughput sequencing platform, a total of 99.16 Gb of clean data were obtained. The Q30 fraction of each sample was greater than or equal to 88.41%. Clean reads quality metrics were shown in Table S2. These results indicated that the quality of the sequencing data was sufficient to continue with the following analysis. 129,544 Unigenes were generated with a total length, mean length, N50 and GC content of 170,757,704 bp, 1,318 bp, 2,090 bp and 41.59%, respectively. The distribution of the Unigenes lengths was shown in Fig. 2. The Unigenes quality metrics were shown in Table 1.

Finally, a total of 78,179 Unigenes were annotated according to seven functional databases. Among all annotated Unigenes, 47,750 were annotated via the SwissProt database, 53,683 were annotated via the KEGG database, and 23,382 were annotated by the GO database (Table 2).

Analysis of DEGs under different light intensities

Based on their gene expression levels, DEGs between groups were identified using the methods described above. A total of 1,938 DEGs were identified among all Unigenes. Comparison of the DEGs in the three different light intensity conditions revealed 1,249, 242 and 698 DEGs when the following conditions were compared: HL vs LL, HL vs ML and ML vs LL, respectively. The distributions of the upregulated and downregulated Unigenes are shown in Fig. 3.

According to GO functional analysis and significant enrichment analysis of the DEGs, three comparisons showed a similar functional enrichment (Fig. 4). The “biological process” category was the most enriched, followed by the “cellular component” category. In the biological process category, the “metabolic processes”, “cellular processes” and “single-organism processes” terms were significantly enriched. In the cellular component category, “cell”, “cell part” and “membrane” were the three largest categories. In the molecular function category, the two largest categories were “binding” and “catalytic activity”.

According to KEGG functional enrichment analysis, the “Metabolism” pathway was the most enriched pathway in all three comparisons. A total of 110 relevant metabolic pathways were identified by KEGG pathway analysis (Fig. 4).

To further investigate the Unigenes involved in flavonoid biosynthesis under different light conditions, 13 DEGs related to flavonoid biosynthesis were screened and identified through GO and KEGG (Table S3). Their FPKM (Fragments Per Kilobase of transcript per Million fragments mapped) were showed in Fig. 5 and sequence data were listed in Supplemental Files.

RT-PCR analysis of DEGs involved in flavonoid biosynthesis

A total of 7 of 13 DEGs related to flavonoid biosynthesis were successfully confirmed by RT-PCR. CtHCT1 and CtHCT2 were significantly downregulated in the LL condition compared to their expressions in the other groups. However, CtHCT3 was downregulated in ML and upregulated in LL. Regarding the FLS and ANS genes, CtFLS1 and CtANS1 were upregulated in LL, but CtFLS2 and CtANS2 were downregulated in LL compared to their expressions in other treatment groups (Fig. 6). These results indicated that a decrease in light intensity had a remarkable impact on the expressions of flavonoid biosynthesis related genes in safflower flowers, and the decrease in light intensity to different degrees resulted in different effects on the expressions of homologous genes.

Analysis of flavonoid content under different light intensities

The total flavonoid content reached a maximum in the ML condition, followed by that in the HL condition, and the lowest flavonoid content was measured in the LL condition. However, the differences of total flavonoid contents among three conditions were not significant. The HSYA content showed little change between the ML and LL conditions but was significant decreased in the HL condition. However, kaempferol showed a different accumulation pattern in three conditions compared with HSYA; it was the highest in the ML condition and was distinctly lower in both the HL and LL conditions (Fig. 7). These results indicated that a reduction in the light intensity within a suitable range was advantageous for flavonoid accumulation in safflower. In addition, the kaempferol concentration showed a significant negative relationship with CtHCT3 expression (P < 0.01), but no significant relationships with other DEGs. The authors suggest that flavonoid biosynthesis was affected by the differential expression of HCT genes, especially CtHCT3, in safflower flowers grown under different light intensities.