Transcriptome profiling of anthocyanin-related genes reveals effects of light intensity on anthocyanin biosynthesis in red leaf lettuce

Plant materials

Lactuca sativa L. var. “Capitata” was grown in the green house of the Luoyang Normal University (China). Seedlings were grown under white light at an intensity 40 μmol m⁻² s⁻¹ and a 16/8 h day/night photoperiod, with the temperature maintained at 22 °C/20 °C. Light intensity was measured using a FGH-1 photosynthetic radiometer (Beijing Normal University Photoelectric Instrument Factory, Beijing, China). Forty days after germination, one group of seedlings was supplied with a light intensity of 100 μmol m⁻² s⁻¹, and the other group was used as a control. After three days of treatment, the two leaves at the top of plant were sampled, and the leaves from five individuals were pooled in each group. Samples were frozen in liquid nitrogen and stored at −80 °C. Quantitative (q) analysis of anthocyanin was performed as described by (Zhang et al., 2015), and cyanidin-3-O-glucoside (Cy3G) was used as standards for quantification. The soluble sugar content was detected using the anthrone method (Li, 2003).

cDNA library construction and sequencing

Total RNA was extracted using CTAB-LiCl method (Gambino, Perrone & Gribaudo, 2008), and genomic DNA contamination was removed using DNase I (TaKaRa, Dalian, China). RNA quality was verified by agarose gel electrophoresis and a Bioanalyzer 2100 (Agilent, CA, USA). For cDNA library construction, mRNA was enriched with oligo (dT) magnetic beads and then broken into smaller pieces using fragmentation buffer. The first strand cDNA was reversed-transcribed by random hexamers and small fragment as templates. This was followed by second strand cDNA synthesis using DNA Polymerase I and RNase H. The double strand cDNA was ligated to paired-end adapters, and suitable fragments were selected and enriched with PCR amplification. cDNA libraries were sequenced performed using an Illumina HiSeq2000 platform at BGI Co., Ltd. Sequence data were deposited in the NCBI database with accession numbers SRR5868088 and SRR5943715.

De novo assembly and gene annotation

Raw reads were filtered by discarding the following: reads with adaptors, reads with unknown nucleotides (Nts) larger than 5%, and low quality reads (more than 10% of bases with a Q score ≤ 20). The remaining clean reads from each sample were assembled using Trinity v2.0.6 software (Grabherr et al., 2011). TGICL v2.06 was then used to construct a non-redundant (Nr) unigene set from the two assembled datasets (Pertea et al., 2003). Unigenes were aligned to functional databases to obtain gene functions, mainly including the NCBI Nr, Nt, InterPro, SWISS-PROT, and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases by BLAST v2.2.23 with an E-value threshold <1e⁻⁵. Based on the Nr annotation, gene ontology (GO) annotation (http://www.geneontology.org) was performed using the BLAST2GO program v2.5.0 (Conesa et al., 2005).

Analysis of DEGs

Clean reads were mapped to unigenes using Bowtie2 v2.2.5 (Langmead & Salzberg, 2012), and gene expression level was calculated based on fragments per kb per million fragments (FPKM) method in RSEM v1.2.12 (Mortazavi et al., 2008). Significantly differentially expressed genes (DEGs) were scanned among samples under low and high light using EBSeq package v1.7.1 (Leng et al., 2013), with a threshold of an absolute log2 ratio ≥ 2 and a false discovery rate (FDR) significance score <0.001. Based on the KEGG and GO annotation, we classified DEGs and performed functional enrichment using phyper within R (R Development Core Team, 2008). The terms in which FDR not larger than 0.001 are defined as significant enrichment.

Prediction and sequence analysis of TF

ORF of each unigene was predicted using getorf EMBOSS:6.5.7.0 (Rice, Longden & Bleasby, 2000), and hmmsearch v3.0 was adopted to identify TF by aligning ORF to domain from PlantTFDB (Mistry et al., 2013). Neighbor-joining phylogenetic analysis was carried out using MEGA5 (Tamura et al., 2011).

Quantitative-PCR Analysis

Extraction of total RNA and elimination of genomic DNA contamination was performed as above, and 2 μg of total RNA was used as template and the first strand cDNA was synthesized using the PrimeScript™ II first strand cDNA Synthesis Kit (TAKARA, Dalian, China) with primer dT, according to the manufacturer’s instructions. Following 10 times dilution, the cDNA was used for quantitative PCR. qPCR was performed using the SYBR® Premix Ex Taq™ II kit (Tli RNaseH Plus; TAKARA, Dalian, China) according to the manufacturer’s instructions. Gene specific primers are shown in Table S1. The thermal-cycling conditions were as follows: an initial heat denaturing step at 95 °C for 3 min; then 40 cycles of 95 °C for 10 s, 55 °C for 20 s, and 72 °C for 20 s. Each sample was amplified in three independent replicates. Gene transcription levels were calculated using the 2^−ΔΔCT comparative threshold cycle (Ct) method (Livak & Schmittgen, 2001), and actin was used as an internal control to normalize the relative expression levels of the analyzed genes (Borowski et al., 2014).

Light intensity determines anthocyanin content

After three days irradiation at a high light intensity of 100 μmol m⁻² s⁻¹, the plant turned red, while the plant grown at low light intensity of 40 μmol m⁻² s⁻¹ still exhibited a pale green color (Figs. 1A and 1B). High-performance liquid chromatography analysis revealed that anthocyanins accumulated in the plants under high light at a level of 2.1 mg^−g, but were barely detected in those grown under low light (Fig. 1C). Our results indicated that anthocyanins were induced as light intensity increased to 100 μmol m⁻² s⁻¹, and which resulted in the red coloring in leaves of lettuce. Soluble sugar content was also increased with increasing light intensity (Fig. 1D). The content increased 2.9-fold when irradiated under 100 μmol m⁻² s⁻¹ in comparison to 40 μmol m⁻² s⁻¹.

Transcriptome sequencing and de novo assembly

After filtering raw data, we generated 59,774,792 clean reads, which included 8,966,218,800 nt from samples under light intensity of 100 μmol m⁻² s⁻¹, and 58,856,106 clean reads that included 8,828,415,900 nt from samples under light intensity of 100 μmol m⁻² s⁻¹. Trinity software was used to perform de novo assembly with the clean reads, generating 69,733 transcripts with an N50 of 1,500 bp in the high light sample, and 79,037 transcripts in the low light sample with an N50 of 1,460 bp. Redundancy in the transcripts was removed using Tgicl, finally resulting in the generation of 62,111 unigenes with an N50 of 1,681 bp.

A total of 24,129 unigene sequences (38.8%) had a length between 200 and 500 Nts, 13,180 unigenes (21.2%) were between 500 and 1000 nt in length, 16,058 unigenes (25.8%) were between 1,000 and 2,000 nt in length, and 8,744 unigenes (14.1%) were longer than 2,000 nt.

After assembly, gene functions were predicted by querying seven public databases, and a total of 48,435 unigenes (77.98%) were functionally annotated. Among them, 45,046 unigenes (72.52%) obtained hits in the Nr database, 31,786 obtained hits in the SWISS-PROT database, and 33,713 unigenes obtained hits in the InterPro database. TransDecoder v3.0.1 (https://transdecoder.github.io) was used to predict the ORFs. Overall, 22,623 genes were predicted to contain full-length ORFs, accounting for 36.42% of the assembled unigenes. Based on the Nr annotation, the top sequence matches obtained from BLASTX are shown in Fig. S1. The red leaf lettuce sequences had the highest similarity to Cynara cardunculus var. Scolymus sequences (57.02%), followed by Vitis vinifera (5.29%), Sesamum indicum (2.53%), Coffea canephora (1.95%), and Nicotiana sylvestris (1.57%).

Differentially expressed genes

A total of 3,929 unigenes were recovered as differentially expressed between samples grown under low and high light conditions. Compared with the low light sample, 1,377 unigenes were up-regulated and 2,552 unigenes were down-regulated in the high light sample (Fig. 2). To further understand the biological functions of the DEGs, we performed pathway analysis based on the KEGG database. A total of 2,891 unigenes were assigned to six categories including 132 KEGG pathways. KEGG pathway enrichment analysis was performed based on hypergeometric tests. The DEGs between the two samples were significantly enriched in 14 pathways (Fig. 3). Among them, four pathways including “phenylpropanoid biosynthesis (ko00940),” “flavonoid biosynthesis (ko00941),” “flavone and flavonol biosynthesis (ko00944),” and “anthocyanin biosynthesis (ko00942)” are closely associated with anthocyanin biosynthesis, and three glucide metabolic pathways including “other glycan degradation (ko00511),” “starch and sucrose metabolism (ko00500),” and “galactose metabolism (ko00052)” may be related to the synthesis of substrates.

We assigned 1,257 of the 3,929 DEGs to three main GO categories including “molecular functions,” “biological processes,” and “cellular components” (Fig. S2). Among them, 742 unigenes were grouped in the category “cellular components,” 992 unigenes in “molecular function,” and 884 unigenes in “biological processes.” Under the GO category “molecular functions,” the “polynucleotide adenylyltransferase activity” and “oxidoreductase activity” were the most highly enriched terms. Under the category “Molecular function,” the “flavonoid biosynthetic process” were the most highly enriched term.

Analysis of genes involved in anthocyanin biosynthesis and transport

Using gene annotation and phylogenetic analysis, we identified putative structural genes involved in anthocyanin synthesis and transport (Table 1). The accession numbers in GenBank were MF579543–MF579560. In terms of anthocyanin structural genes, a total of nine genes covered each step of the anthocyanin biosynthetic pathway. Among them, the CHS and 3GT gene family contained two members, while the other genes contained only one member. All nine candidate genes were up-regulated under high light, and their normalized transcript levels were two- to ninefold higher than under low light. Among the CHS and 3GT gene family members, Unigene12000_All and CL4808.Contig1_All displayed greater transcript abundance and more obvious growth at the transcription level. It is indicated that Unigene12000_All and CL4808.Contig1_All may be more insensitive to light intensity and contribute more to anthocyanin biosynthesis. In brief, all the putative structural genes were co-up-regulated and the anthocyanin pathway was active under high light conditions.

We also detected the differential expression of two anthocyanin transport genes, GST (Unigene10814_All) and MATE (Unigene12020_All), with transcript levels 7.1- and 1.1-fold higher under high light conditions. GST possessed higher transcript abundance and was significantly up-regulated, suggesting that anthocyanins might be primarily transported by GST.

Transcription factors regulating anthocyanin biosynthesis

We predicted a total of 291 TFs belonging to 39 families (Table S2). The MYB gene family represented the largest group containing 33 members, followed by the AP2-EREBP (30 members), MYB-related (26 members), and bHLH gene families (19 members). In model plants, TFs from the MYB, bHLH, and WD40 families regulate transcripts of anthocyanin structural genes. In the 33 MYBs filtered from the detected DEGs, 14 were up-regulated and 19 were down-regulated under high light conditions. Among the up-regulated genes, Unigene12430_All, Unigene12294_All, and Unigene23058_All formed a clade with anthocyanin regulatory genes from Arabidopsis, grape, and Antirrhinum, which have been shown to play a central role in regulating LBGs. Unigene24751_All and CL6440.Contig1_All were closely associated with MYB12 of Arabidopsis, which mainly regulates EBGs (Fig. 4A). Following treatment with high light, the transcription levels of MYBs were up-regulated 1.9- to 6.0-fold.

In our study, six bHLH TFs were up-regulated and 13 were down-regulated under high light conditions. Of these, Unigene13011_All was predicted as a anthocyanin regulatory gene and was up-regulated 2.6-fold under high light. Sequence alignment revealed that it contained the BOX18, BOX19, and HLH motifs, which were depicted as conserved in sub group III of the bHLH gene family. It grouped with bHLH in the phylogeny (Fig. 4B), which has been proven to regulated anthocyanin in apple, Petunia, and Arabidopsis, and was most closely related to DvIVS with 62.6% amino acid similarity.

HY5 is a member of the bZIP gene family and acts downstream of the light receptor network and directly affects the transcription of light-induced genes. It was previously verified to regulate anthocyanin structural genes and PAP1 expression by directly binding to their promoters (Shin, Park & Choi, 2007; Shin et al., 2013). In the bZIP family, 14 genes were up-regulated and 19 were down-regulated. Phylogenetic analysis placed Unigene19629_All in the same cluster as the HY5 proteins from different plant species (Fig. S3). Unigene19629_All showed 65.5% amino acid similarity with AtHY5 in Arabidopsis, the transcription level of which was up-regulated 1.7-fold under high irradiance.

Quantitative real-time-PCR validation of DEGs

To further validate the comparative transcriptome results, the transcript level variances of the 18 putative genes involved in anthocyanin synthesis, transport, and regulation between the low and high light conditions were detecting using quantitative real-time (qRT)-PCR analysis. The results indicated that the transcript levels of 17 genes were significantly up-regulated in the high light samples (Fig. 5), which is in agreement with the alterations in gene expression detected by the transcriptome analysis. However, the MYB gene CL6440.Contig1_All was not significantly up-regulated under light intensity of 100 μmol m⁻² s⁻¹. Overall, these results indicate that the transcriptomic profiling data correlate with the light intensity responses of red leaf lettuce.