The lily Tango Series cultivar “Tiny Padhye” was grown in a greenhouse at the Chinese Academy of Agricultural Sciences (Beijing, China). Flowers sampling was performed at four different developmental stages (Figure 1A): Stage 1 (S1; bud length ~1.5 cm, and no anthocyanins visible on tepals), stage 2 (S2; anthocyanins visible on tepal bases), stage 3 (S3; 1 day before anthesis where the lower half of tepals was fully pigmented), and stage 4 (S4; 0 days post-anthesis). At each stage, samples were taken from the upper inner tepals of 20 flowers and pooled together as one biological sample; three independent biological replicates were collected for each stage. This process was repeated for the inner tepal bases. The collected samples were immediately frozen in liquid nitrogen and then stored at −80°C until use.

Frozen lily tepals were ground into powder in liquid nitrogen. Anthocyanins were extracted with a solvent mixture of trifluoroacetic acid, methanol, methanoic acid, and water (1:70:2:27, v:v:v:v). The extract was then analyzed by high-performance liquid chromatography (HPLC) on a Waters 2695 Alliance HPLC instrument (Waters Corporation, Milford, MA, USA) connected to a Waters 2998 photodiode array detector. The separation was performed on a Waters SunFire C18 column (150 × 4.6 mm I.D., 5 μm) with a C18 guard column at 30°C. Products were detected at 530 nm. Anthocyanins were eluted with a gradient mobile phase consisting of solvent A (HCOOH:H₂O = 0.1:100, v:v) and solvent B (HCOOH:CH₃CN = 0.1:100, v:v) at a flow rate of 1.0 ml·min⁻¹. The solvent gradient program is shown in Table S1. Six authentic standards [cyanidin (Cy), pelargonidin 3,5-diglucoside (Pg3,5dG), cyanidin 3-O-β-glucoside (Cy3G), cyanidin 3-O-β-rutinoside (Cy3R), and peonidin 3-O-β-glucoside (Pn3G)] were used to identify peaks. Chl content was determined by measuring the absorbance of the supernatant at 663 and 645 nm according to the protocol of Arnon (1949).

Total RNA was extracted from each sample using an RNAprep pure Plant Kit (Tiangen Biotech Co., Ltd., Beijing, China), according to the manufacturer's instructions. Total RNA was extracted from upper tepals at stage 2 (S2-U) and from basal tepals at the four different stages (S1-D, S2-D, S3-D, and S4-D). Three biological replicates were used for each sample. A total of 15 RNA-seq libraries were constructed from these samples using a ScriptSeq mRNA-Seq Library Preparation Kit (Epicentre Biotechnologies, Madison, WI, USA) according to the manufacturer's protocol. The libraries were sequenced to generate 150 paired-end raw reads on an Illumina Hiseq 4000 platform. After clean reads were generated by removing adapters, low-quality reads, and ambiguous reads from the raw data, transcriptome assembly was accomplished using Trinity software as previously described for de novo transcriptome assembly without a reference genome (Grabherr et al., 2011).

The assembled unigenes were searched against six public databases, including the NCBI Non-Redundant Protein Sequences (NR) database, NCBI Nucleotide Sequences (NT) database, Protein Family (PFAM) database, Swiss-Prot protein database, Gene Ontology (GO) database, Eukaryotic Ortholog Groups (KOG) database, and Kyoto Encyclopedia of Genes and Genomes (KEGG) database.

The differentially expressed genes (DEGs) between pairs of samples were identified and filtered with the R package DESeq (Anders and Huber, 2010). The DEGs between two samples were determined based on |log2 (foldchange)| > 1 and padj < 0.05. The heatmap displays of the Trimmed Mean of M-values (TMM) normalized against the Fragment Per Kilobase of transcripts per Million mapped reads (FPKM) were performed created using the R package pheatmap (https://cran.r-project.org/src/contrib/Archive/pheatmap/) (Kolde, 2012).

A total of 10 unigenes related to anthocyanin biosynthesis and nine unigenes related to Chl metabolism were chosen for qRT-PCR analyses. qRT-PCR analyses were performed using SYBR® Premix Ex Taq™ II (Tli RNaseH Plus) (Takara, Dalian, China) and a Bio-Rad iQ5 Gradient RT-PCR system with the following reaction conditions: Denaturation at 95°C for 30 s and 40 cycles of amplification (95°C for 5 s, 60°C for 30 s). The lily Actin gene was used as an internal control for normalization. Relative expression levels of target genes were calculated using the 2^−ΔΔCT method (Livak and Schmittgen, 2001) against the internal control. The gene-specific primers are shown in Table S2. Experiments were performed with three independent biological replicates and three technical replicates.

To identify TFs, assembled unigenes were searched against the Plant Transcription Factor Database PlnTFDB (http://plntfdb.bio.uni-potsdam.de/v3.0/downloads.php) using BLASTX with an E-value cut-off of ≤ 10⁻⁵.

To investigate the co-expressed gene networks in tepal development with a particular focus on the transcriptional architecture of anthocyanin biosynthesis and Chl metabolism, we constructed gene co-expression networks from the DEGs using the WGCNA package (Langfelder and Horvath, 2008). The networks were visualized using Cytoscape _v.3.0.0 (Hollender et al., 2014).

During bicolor tepal development in Asiatic “Tiny Padhye,” the upper tepals underwent a rapid color change from green to white, whereas the tepal bases changed from green to purple (Figure 1A). To investigate these physiological changes, the anthocyanin, and Chl contents of upper and lower tepals were assessed at four different stages of tepal development.

Using HPLC, a single anthocyanin (cyanidin 3-O-β-rutinoside) was detected in the pigmented tepal bases of “Tiny Padhye” (Figure S1). No anthocyanins were detected in the upper tepals at any stage of tepal development (Figure 1B). The anthocyanin levels in tepal bases increased dramatically from S1 to S3, and then gradually decreased at S4 (Figure 1B). As shown in Figure 1C, Chls accumulated at the early developmental stages (S1 and S2) in the upper and lower tepal parts before decreasing to extremely low levels at the late stages (S3 and S4).

The transcriptomes of the 15 “Tiny Padhye” samples were separately obtained using Illumina technology. Overall, ~772.62 million 150-nt paired-end raw reads were generated (Table 1). All raw reads have been deposited in the NCBI Sequence Reads Archive (SRA) under the accession number SRP093907. After removing low-quality sequences, adapters, and ambiguous reads, we obtained a total of 713.12 million high-quality clean reads (Table 1). These clean reads were assembled into 400,850 transcripts with a mean length of 668 nt and an N50 length of 1,043 nt (Table 1). These transcripts were then assembled into 295,787 unigenes with an average length of 544 nt and an N50 length of 718 nt (Table 1).

To annotate the transcriptome with putative functions, the assembled unigenes were searched against four public databases: NR, NT, PFAM, and SWISS-PROT. Among them, 66,531, 25,508, 47,404, and 40,271 unigenes were annotated to the NR database, NT, PFAM, and SWISS-PROT databases, respectively (Table S3). To further illustrate the main biological functions of the transcriptome, GO, KOG, and KEGG pathway analyses were performed. GO analysis provides a description of gene products in terms of their associated Biological Process (BP), Cellular Component (CC), and Molecular Function (MF) (Berardini et al., 2004) (Figure 2A). A total of 48,927 unigenes were categorized into 47 major functional groups. Cellular process (GO:0009987), cell (GO:0005623), and binding (GO:0005488) were the most highly represented GO terms in BP, CC, and MF, respectively (Table S3 and Figure 2A).

In addition to the GO analysis, KOG analysis was performed to further evaluate the function of the assembled unigenes. A total of 17,714 annotated unigenes were grouped into 26 KOG classifications (Table S3 and Figure 2B). The three most abundant groups represented in the transcriptome were group J (translation, ribosomal structure, and biogenesis), group O (post-translational modification, protein turnover, and chaperones) and group R (general function prediction only) (Figure 2B).

In the KEGG pathway analysis, a total of 18,675 unigenes were mapped to 278 KEGG pathways (Table S3). The pathways with the highest unigene representations were ribosome (ko03010, 860 unigenes), followed by carbon metabolism (ko01200, 828 unigenes), and biosynthesis of amino acids (ko01230, 723 unigenes). Notably, 61 unigenes were involved inanthocyanin biosynthesis and 106 unigenes were involved in Chl metabolism (Table 2).

The DEGs between upper tepals and tepal bases at different developmental stages were identified and filtered with the R package DESeq (Anders and Huber, 2010) (Tables S4–S7). A total of 16,182 DEGs were detected in at least one pair-wise comparison (S2-D vs. S1-D, S2-D vs. S2-U, S3-D vs. S2-D, and S4-D vs. S3-D) (Figure 3A). Among the four comparisons, the largest number of DEGs was between the S2-D and S3-D libraries (9,958), with 4,775 down-regulated and 5,183 up-regulated unigenes (Figure 3B). The smallest number of DEGs was between the S2-D and S2-U libraries (2,036), with 671 down-regulated and 1,365 up-regulated unigenes (Figure 3B). Among the DEGs, 112 were significantly differentially expressed in all pair-wise comparisons (Figure 3A).

During bicolor tepal development in Asiatic “Tiny Padhye,” anthocyanins were heavily accumulated in the tepal bases but less so in the upper tepals. Therefore, we investigated genes encoding enzymes involved in anthocyanin biosynthesis. In this study, 61 unigenes encoding eight putative enzymes involved in anthocyanin biosynthesis were identified from the transcriptome (Table 2), and 21 of them were identified as DEGs (Table 2). As shown in Figure 4, most DEGs including LhCHS1 (c4453117_g1), LhCHS2 (c8851122_g2), LhCHS3 (c1209245_g1), LhCHS4 (c1209245_g2), LhCHI (c2301271_g2), LhF3H (c6307121_g1), LhF3′H (c9329113_g1), LhDFR (c4316912_g1), LhUFGT1 (c3956911_g1), and Lh3RT (c8092117_g1) showed significantly higher expression in tepal bases at S2 and S3 than in the other samples, and extremely low expression in tepal bases at S1 and S4 and in upper tepals at S2. The results of qRT-PCR analyses not only confirmed these results, but also showed that anthocyanin structural genes had extremely low expression levels in upper tepals at all four stages (Figure 5).

Degreening is a significant developmental change in the bicolor tepals of Asiatic “Tiny Padhye.” Therefore, we studied unigenes involved in Chl metabolism in detail. In this study, 106 candidate genes related to Chl metabolism were identified, and 28 of them were identified as DEGs (Table 2). These genes showed two distinct expression patterns (Figure 6). All of the DEGs related to Chl biosynthesis, except for LhDVR (c118209_g1), showed significantly higher expression in tepal bases at S1 and S2 and in the upper tepals at S2 than in the other stages. Their expression levels in tepal bases at S3 and S4 were extremely low (Figure 6). Conversely, most DEGs involved in Chl degradation showed the opposite pattern (Figure 6). The transcript levels of LhPPH1 (c119711_g1), LhPPH2 (c119711_g2), LhPaO (c120481_g1), LhRCCR (c111130_g1), and LhSGR (c124306_g4) were low at the early stages (S1 and S2) in the lower tepal parts, but high at the late stages (S3 and S4) (Figure 6). The results of qRT-PCR analyses not only confirmed these results but also showed that these DEGs shared similar expression patterns in upper parts during all four stages (Figure 7).

To further investigate candidate genes related to anthocyanin biosynthesis and Chl metabolism during bicolor tepal development, co-expression gene network modules were constructed using WGCNA. The co-expression network constructed based on the 16,182 DEGs revealed 26 modules (Figure 8). In the dendrogram, each branch represents a module and each leaf constitutes one gene.

The Red module contained 12 out of 21 anthocyanin-related DEGs (Table 3), indicating that the 405 unigenes in the Red module play important roles in anthocyanin accumulation in lily bicolor tepals. The majority of the Chl biosynthesis-related DEGs were in the Yellow2 (9/19) module and more than half of the Chl degradation-related DEGs were in the Turquoise1 module (4/7) (Table 3). Therefore, the unigenes in the Yellow2 and the Turquoise1 modules were potentially involved in regulating Chl metabolism in lily tepals.

To provide a system view of the regulatory network responsible for controlling anthocyanin and Chl levels, networks were extracted from the Red, Yellow2, and Turquoise1 modules. To search for regulatory genes potentially involved in anthocyanin biosynthesis in lily tepals, we constructed a subnetwork from the Red module using 12 anthocyanin biosynthesis-related genes as the seed nodes. As a result, 12 anthocyanin genes and 21 TFs belonged to this subnetwork (Figure 9A and Table 4). These TFs were classified into 10 putative TF families, with the four largest TF families being the MYB (three unigenes), WRKY (three unigenes), Dof (three unigenes), and ERF (three unigenes) families (Table 4 and Table S8). The expression profiles of these TFs potentially related to anthocyanin biosynthesis were hierarchically clustered and plotted in a heatmap (Figure 9D). The expression patterns of all 21 TFs, except for c123181_g3 and c114287_g2, were positively correlated with those of anthocyanin biosynthetic genes.

Using these Chl biosynthesis-related genes as bait, a huge subnetwork was extracted from the Yellow2 module. This subnetwork contained nine Chl biosynthesis-related genes and 78 TFs belonging to 25 TF families (Figure 9B and Table 4). Among the TF families, bHLH (21 genes) was the largest group, followed by MYB (six genes), C2H2 (five genes), and C3H (five genes) (Table 4 and Table S9). Among the TFs, 60 had expression patterns that were positively correlated with Chl content and the remainder had expression patterns that were negatively correlated with Chl content (Figure 9E).

In the Turquoise1 module, four Chl degradation-related unigenes were highly interconnected and 232 TFs, belonging to 35 different families, had interactions with these Chl degradation-related unigenes (Figure 9C and Table 4). The main TF family was the WRKY family (65 unigenes), followed by the NAC TF family (29 unigenes), and the bHLH family (16 unigenes) (Table 4 and Table S10). Notably, these TFs also showed two distinct expression patterns (Figure 9F).

To date, two previous studies have identified anthocyanin biosynthetic and regulatory genes in lily flowers based on transcriptome analysis by RNA-seq (Zhang et al., 2015; Suzuki et al., 2016). However, those studies were limited by the small number of samples used and the small scale of transcriptome data obtained. For example, to investigate the molecular mechanisms responsible for bicolor tepal development in lilies, Suzuki et al. (2016) analyzed the global transcription of pigmented and non-pigmented tepal parts from “Lollypop” at stage 3 with two replicates. This resulted in ~50 million raw reads and a total of 39,426 unigenes through de novo assembly. However, as that study did not provide a global view of transcriptome dynamics over the key tepal developmental stages of anthocyanin biosynthesis, it was hard to perform a co-expression network analysis to identify new candidate target genes. Our transcriptome profiling and network analysis differ from these prior studies in several ways. First, we not only analyzed the global transcription of pigmented and non-pigmented tepal parts from “Tiny Padhye” at stage 2 but also constructed a high-resolution transcriptome atlas of lily inner tepal development using a time series of tepal base samples taken from S1 to S4. Thus, ~772.62 million raw reads were generated and 295,787 unigenes were assembled. This large-scale transcriptome analysis serves as a valuable resource for analyzing gene function on a global scale, and for elucidating the developmental processes of anthocyanin biosynthesis and Chl metabolism. Second, all the genes involved in Chl metabolism were identified in Lilium spp. using the KEGG analysis, and their expression profiles during tepal development were analyzed to identify DEGs. Additionally, using WGCNA, a number of candidate TFs involved in anthocyanin biosynthesis and Chl metabolism were identified in lilies. Overall, with multiple time-points and a co-expression analysis, we can say that the present work is the first dynamic transcriptomic study of lily flower color development.

In this study, we identified 61 unigenes encoding eight candidate enzymes in the anthocyanin biosynthetic pathway from the lily transcriptome. Among these unigenes, 21 were DEGs and most of them were expressed at high levels coordinately in tepal bases, while showing extremely low expression levels in the upper tepals at each stage. Similar expression patterns have been reported for “Lollypop” tepals (Suzuki et al., 2016). These results indicate that the bicolor trait of “Tiny Padhye” is caused by the transcriptional regulation of anthocyanin biosynthesis genes rather than the PTGS of CHS genes.

Many studies have demonstrated the crucial role of TFs in regulating anthocyanin biosynthesis in plants. In this study, we identified 10 TF families comprising 21 TF unigenes that represented potentially important regulators of anthocyanin biosynthesis in Lilium spp. The R2R3-MYB family plays a key role in regulating the spatiotemporal expression of anthocyanin biosynthetic genes in plants (Gonzalez et al., 2008; Zhao and Tao, 2015). This family is further classified into several subgroups, with the members of subgroup 6 positively regulating anthocyanin biosynthesis (Dubos et al., 2010). In this study, the sequence of c8752811_g1, designated as a subgroup 6 R2R3-MYB, was found to be the same as that of LhMYB12-Lat, which activates anthocyanin accumulation in lily tepals (Yamagishi et al., 2014). This result shows that the WGCNA is a powerful method of identifying TFs relevant to a specific pathway. Several anthocyanin repressors have been characterized, including small R3-MYBs, AtMYBL2, and subgroup 4 R2R3-MYBs. Some R3-MYBs [CAPRICE (CPC) from A. thaliana, and MYBx from petunia (P. hybrida)] contain a single MYB DNA-binding domain (DBD) and a conserved amino acid motif ([DE]Lx2[RK]x3Lx6Lx3R) required to bind bHLH partners, and are found to regulate anthocyanin biosynthesis by competing with R2R3-MYB activators for interaction with bHLH partners (Zhu et al., 2009; Albert et al., 2014). AtMYBL2 from Arabidopsis lacks an intact R2-MYB repeat, and contains a repression motif (TLLLFR) to actively repress transcription (Dubos et al., 2008; Matsui et al., 2008). Subgroup 4 R2R3-MYBs [Fa/Fc-MYB1 from strawberry (Fragaria spp.) and Ph MYB27 from petunia (P. hybrida)] contain an ERF-associated amphiphilic repression (EAR) motif that is essential for the repressive function (Aharoni et al., 2001; Salvatierra et al., 2013; Albert et al., 2014). The unigene c120204_g1(LhMYB3) was designated as a subgroup 4 R2R3-MYB, suggesting that this MYB suppressor might be involved in bicolor development. We also identified one unigene (c115761_g1), LhMYBP, as a homolog of ZmP. In maize, ZmP activates the promoters of CHS, CHI, F3H, and FLS, but it does not control anthocyanin accumulation (Grotewold et al., 1994; Mehrtens et al., 2005). Thus, the function of LhMYBP needs to be further studied. In our study, three putative WRKY TFs were identified including the unigene (c117585_g2) LhWRKY44 that strongly attracted our attention. AtTTG2 (WRKY44) and its homologs have different functions in different species. In Arabidopsis, AtTTG2 is involved in trichome initiation, proanthocyanidin accumulation in the seed coat, and seed development (Johnson, 2002; Ishida et al., 2007; Dilkes et al., 2008). However, in Brassica napus, Bna.TTG2 genes, homologs of AtTTG2 (WRKY44), are involved in the salt stress response (Li et al., 2015). In P. hybrida, PH3, which is homologous to AtTTG2 (WRKY44), regulates color patterns by acidifying the central vacuole where anthocyanins are stored in epidermal petal cells (Verweij et al., 2016). In the present study, the expression pattern of LhWRKY44 was parallel to those of anthocyanin biosynthetic genes, but its function is still unknown in Lilium spp. and requires further research.

We identified 106 candidate structural genes involved in Chl metabolism in the transcriptome dataset, and 28 of these unigenes were identified as DEGs. The expression of most DEGs in Chl biosynthesis declined at S3 and S4, correlating strongly with Chl content. The expression pattern of most DEGs involved in Chl degradation was opposite to that of DEGs involved in Chl biosynthesis. Overall, based on the gene expression and Chl content analysis, our results indicate that high-level expression of Chl degradation genes and low-level expression of Chl biosynthetic genes lead to the absence of Chl from “Tiny Padhye” tepals after flowering (Figure 10).

TFs are crucial regulatory proteins that mediate transcriptional regulation. Many previous studies have identified the TFs responsible for controlling Chl metabolism in leaves (Lim, 2003; Lin et al., 2015; Wen et al., 2015). Three TFs, ANAC046, ETHYLENE INSENSITIVE3, and ORE1, positively regulate Chl degradation by binding to the promoter regions of NYC, SGR1, and PaO (Qiu et al., 2015; Oda-Yamamizo et al., 2016). FLU, a negative regulator of Chl biosynthesis in A. thaliana, may mediate its regulatory effect through interaction with enzymes (GSA and CHlH) involved in Chl biosynthesis (Meskauskiene et al., 2001). However, the molecular mechanisms regulating Chl metabolism in lily petals are still unknown. In this study, we identified 78 TFs in the Chl biosynthesis co-expression network. These TFs were categorized into 25 putative TF families (Table S9). The main TF family was the bHLH family (21 unigenes), followed by the MYB family (six unigenes), the C2H2 family (five unigenes), and the C3H family (five unigenes). These results indicate that there is a complicated regulatory network controlling Chl biosynthesis in lily tepals.

Several TF families, especially the NAC and WRKY families, have members that control Chl degradation in leaves (Lim, 2003; Wen et al., 2015). Here, we found 94 TFs belonging to these families in the Chl degradation co-expression network (Table 4; Table S10). This result indicated that, in both leaves and tepals, some members of the NAC and WRKY families have highly conserved functions in regulating Chl degradation. An additional 26 TF families (bHLH, C2H2, C3H, and others) were also detected in this study.

We further analyzed the expression patterns of these TFs, and observed two distinct expression patterns for TFs involved in both Chl biosynthesis and Chl degradation. The expression profiles of some TFs were positively correlated with Chl content, while those of other TFs were negatively correlated with Chl content, suggesting that both transcriptional activation and repression are involved in regulating Chl biosynthesis and Chl degradation in lily tepals.

In this study, we presented a dynamic transcriptome landscape of lily bicolor tepal development using RNA-seq technology. Based on de novo transcriptome analysis and functional annotation of the transcriptome dataset, we identified transcripts encoding most of the known enzymes involved in the anthocyanin biosynthetic and Chl metabolic pathways. Our results indicate that the white upper tepals are caused by a low rate of Chl biosynthesis, a high rate of Chl degradation, and a low rate of anthocyanin biosynthesis. The purple tepal bases are due to a high rate of anthocyanin biosynthesis, a low rate of Chl biosynthesis, and a high rate of Chl degradation (Figure 10). We identified regulatory genes that are likely to be key regulators of anthocyanin biosynthesis and Chl metabolism using WGNCA. Taken together, these results pave the way for the greater understanding of the molecular basis of lily bicolor tepal development and will allow for the identification of candidate genes associated with anthocyanin biosynthesis and Chl metabolism in Lilium spp.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.