Transcriptome Analysis of the Biosynthesis of Anthocyanins in Begonia semperflorens under Low-Temperature and HighLight Conditions

Anthocyanins are considered a stress indicator due to their involvement in the response to many stresses including high light (HL) and low temperature (LT). With the development of transcriptomics, it is necessary to find the different and common points in the mechanisms of LT-induced and HL-induced anthocyanin biosynthesis. In the present study, we determined the transcriptomes of Begonia semperflorens leaves under three different conditions (normal growing conditions (CK), HL, and LT). To validate the differentially expressed genes (DEGs), we selected four core genes involved in anthocyanin biosynthesis to perform real-time reverse transcription-quantitative PCR (RT-qPCR), and then determined anthocyanin content. In total, 94,880 unigenes with a mean length of 635 bp were assembled. The N50 values of the transcripts and unigenes were 2286 bp and 1064 bp, respectively. The functional annotations of the unigenes were analysed against five protein databases. DEGs related to anthocyanin biosynthesis, transportation, and regulation were identified. We also analysed the enrichment pathway, and the differences in mechanisms of anthocyanin induced under low-temperature and high-light conditions are discussed in this paper. This study is the first to examine broad-scale gene expression in Begonia semperflorens. By identifying DEGs regulated by both LT and HL conditions, we found that anthocyanin accumulation was employed as a common strategy by Begonia seedlings in resisting LT and HL stress. By identifying DEGs regulated differently by LT and HL conditions, we found that Begonia seedlings also had some different strategies for resisting LT and HL stresses: anthocyanins were biosynthesized under HL condition, while lignin, proanthocyanidins, and anthocyanins were biosynthesized under LT condition.


Introduction
Plants have developed many strategies to resist stresses and compensate for their immobility.Products of secondary metabolism, such as flavonoids, play many important roles in resisting stress.Flavonoids constitute a diverse family, which includes compounds such as flavones, isoflavones, flavonols, catechins, and anthocyanins.Anthocyanins, a subgroup of flavonoids, are of interest due to their significant roles as pigments and antioxidants.For example, anthocyanins are known to attract seed dispersers and pollinators and to protect plants against high-light irradiation [1].Under stressful conditions, anthocyanins can act as scavengers of free radicals that are produced within cells [2][3][4].
Anthocyanins are induced by many environmental stresses in plants.Low temperature and high light are two common environmental conditions that induce anthocyanin biosynthesis.High excitation pressure is thought to be the common driving force of low temperature-induced and high light-induced anthocyanin biosynthesis in plants [14,15].However, there are still differences between acclimations to low temperature and high light.Adaptation to low temperature is a complex bioprocess and involves many different metabolic pathways.Low temperature could upregulate carbohydrate metabolism, redox adjustment, cell wall remodelling, and defence/detoxification but can downregulate photosynthesis-related protein expression [16].Plant responses to high light involve coordinated chloroplast and extra-chloroplast adjustment processes that typically occur in a chloroplast-centric manner [17].Due to the development of transcriptomics, it is now possible to investigate the differences and commonalities between these mechanisms.
B. semperflorens is a perennial, herbaceous flower of the family Begoniaceae and is widely planted in both tropical and subtropical regions for its colourful leaves and flowers.B. semperflorens 'Super Olympia' has green leaves that accumulate anthocyanins and become red when subjected to several environmental stresses (for example, low temperature, high light, nutrient deficiencies, and wounding).As such, B. semperflorens is an ideal plant for studying the roles of anthocyanins and the mechanisms of anthocyanin biosynthesis.Although generic background information on anthocyanins is available for classical model plants, published gene sequences of members of the Begoniaceae family are sufficiently rare that the genomic characteristics of this family are poorly known.Transcriptome analysis has been shown to be a precise and reliable approach for studying the genomic characteristics of plants for which genomic information is unavailable [18][19][20][21][22][23].In the present study, we determined the transcriptomes of B. semperflorens leaves under different conditions to investigate the biological mechanisms of anthocyanins.Leaves of B. semperflorens plants exposed to low temperatures (LT), high light (HL), and normal growing temperatures (CK) were analysed through transcriptome sequencing.
In total, 94,880 unigenes, with a mean length of 635 bp, were obtained.These unigenes were aligned and analysed with five protein databases, namely, the Non-Redundant protein (Nr), Swiss-Prot, Clusters of Orthologous Groups (COG), Gene Ontology (GO), and the Kyoto Encyclopedia of Genes and Genomes annotations (KEGG) databases.Differentially expressed genes (DEGs) were identified based on comparisons of LT vs. CK and HL vs. CK.DEGs involved in anthocyanin bioprocesses, including biosynthesis, transportation, and regulation, were also identified.To the best our knowledge, this research represents the first broad-scale transcriptome analysis of B. semperflorens.

Plant Material
Seedlings of B. semperflorens 'Super Olympia' were cultured in plastic pots containing a mixture of peat and vermiculite (7:3, v/v) under 25/15 • C (day/night) and 300 µmol•m −2 s −1 for 10 h, with 95% relative humidity for five days.The normal control, low-temperature-treated, and high-light-treated plants were designated CK, LT, and HL, respectively.CK plants were cultured under 25/15 • C (day/night) and 300 µmol•m −2 s −1 for 10 h and exposed to 95% relative humidity every day.LT plants were cultured under 15/5 • C (day/night) and 300 µmol•m −2 s −1 for 10 h and subjected to 95% relative humidity every day.HL plants were cultured under 25/15 • C (day/night) and 500 µmol•m −2 s −1 for 10 h and exposed to 95% relative humidity every day.After 12 days, the second pair of leaves from the growing apex of each plant was removed.For each set of treatment conditions, four parallel samples from four individual plants were prepared.All of the samples were harvested one at a time.The four parallel samples of each treatment were then mixed together as each treatment sample.A leaf sample of at least 0.3 g was prepared from each sample, immediately frozen in liquid nitrogen, and then stored at −80 • C for RNA extraction.The RNA-seq analysis was performed with a single biological replicate per experimental condition.In addition, another CK seedling was used as a biological replicate with the former one.RNA-seq was done and correlations in the RNA-seq data were analysed.

RNA Isolation, Library Preparation, and Sequencing
Total RNA was isolated from plant cells with a Quick RNA isolation kit (Bioteke Corporation, Beijing, China) according to the manufacturer's recommended protocol and then subjected to RNA purification using RNase-free DNase I (TaKaRa, Dalian, China) to remove the genomic DNA contamination.RNA degradation and contamination were monitored by 1% agarose gel electrophoresis.A NanoPhotometer ® spectrophotometer (Implen, Westlake Village, CA, USA) was employed to assess the RNA purity.The RNA concentration was measured using a Qubit ® RNA Assay Kit with a Qubit ® 2.0 Fluorometer (Life Technologies, Carlabad, CA, USA).The RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies, Palo Alto, CA, USA).For each sample, 3 µg of total RNA was pooled and used for subsequent RNA-Seq analysis.Sequencing libraries were generated using the NEB Next ® Ultra™ RNA Library Prep Kit from Illumina ® (New England Biolabs, Ipswich, MA, USA) following the manufacturer's recommendations, and index codes were added to each sample.Briefly, mRNA was purified from total RNA using poly-T oligo-attached magnetic beads.Fragmentation was performed using divalent cations under elevated temperature in NEBNext First-Strand Synthesis Reaction Buffer (5×).The first strand of cDNA was synthesized using a random hexamer primer and M-MuLV Reverse Transcriptase (RNase H-).Second-strand cDNA synthesis was subsequently performed using DNA Polymerase I and RNase H.The remaining overhangs were converted into blunt ends via exonuclease/polymerase treatment.After adenylation of the 3 ends of DNA fragments, NEB Next Adaptors with hairpin loops were ligated to prepare for hybridization.To preferentially select cDNA fragments with a length of 150-200 bp, the library was purified with the AMPure XP system (Beckman Coulter, Beverly, MA, USA).Then, 3 µL of USER Enzyme (NEB) was added to size-selected, adaptor-ligated cDNA, and the resulting mixtures were incubated at 37 • C for 15 min and then 95 • C for 5 min.PCR was performed with Phusion High-Fidelity DNA polymerase, Universal PCR primers and Index (X) Primers.Finally, the PCR products were purified using the AMPure XP system, and the library quality was assessed with an Agilent Bioanalyzer 2100 system.
Clustering of index-coded samples was performed on a cBot Cluster Generation System using a TruSeq PE Cluster Kit v3-cBot-HS (Illumina) according to the manufacturer's instructions.After cluster generation, the library preparations were sequenced on the Illumina HiSeq 2000 platform and paired-end reads were generated.

Quality Control, Transcriptome Assembly, and Gene Functional Annotation
Clean data (clean read pairs) were obtained by removing read pairs containing adapters, read pairs containing poly-N sequences, and low-quality read pairs from the raw data.The steps to filter low quality data were as follows: firstly, the sequences whose N percentages were more than 10% were filtered out; secondly the low quality reads from the whole were filtered out, while making sure that Q30 was above 85%; finally, for a single read with low quality bases whose Q values were lower than 10, if the percentage of the bases were above 50%, it was filtered out.The Q20 and Q30 values, GC contents, and sequence duplication levels were also calculated for the clean data.All of the subsequent analyses were performed with the high-quality, clean data.The left files (read1 files) for all libraries/samples were pooled into one large left.fqfile, and the right files (read2 files) were pooled into one large right.fqfile.Transcriptome assembly was carried out with the left.fqand right.fqfiles using Trinity [24].The default value for min_kmer_cov was set to 2, and all of the other parameters were set to default values.The filtered data were deposited in the National Center for Biotechnology Information (NCBI) Sequence Read Archive (SRA) under accession number SRP126322.These data and the data of the assembled contigs were then deposited in the Gene Expression Omnibus (GEO) repository under accession number GSE107805.The data of the replicated CK sample were also deposited in the SRA database under accession number SPR130727.BLASTX alignment [25] with an E-value ≤ 10 −5 was performed to implement functional annotation.Protein databases, including the Nr [26], Swiss-Prot [27], KEGG [28], COG [29], and GO [30] databases, were employed for Basic Local Alignment Search Tool (BLAST) searches and annotation.

Expression Annotation
Gene expression levels were determined for each sample using RNA-Seq by expectation maximization (RSEM) [31].The clean data were mapped back onto the assembled transcriptome, and the read count for each gene was obtained from the mapping results.Prior to differential gene expression analysis, for each sequenced library, the read counts were adjusted using the edgeR programme package by implementing one scaling normalized factor.Differential expression analysis was performed on pairs of samples using the DEGseq R package [32].The p-value was adjusted with the q-value [33].A q-value < 0.005 and an absolute log 2 value (fold change) > 1 were set as thresholds for significantly different expression.
The GO enrichment analysis of the DEGs was implemented using the GOseq R package based on the Wallenius' noncentral hypergeometric distribution [34], which adjusts for gene length bias.The KEGG [35] database is a resource for understanding the high-level functions and uses of biological systems (such as cells, organisms, and ecosystems) based on molecular-level information.The KEGG database is particularly useful for large-scale datasets generated via genome sequencing and other high-throughput technologies [28].We used KOBAS software [36] to test the statistical enrichment of DEGs in KEGG pathways.C, 500 µmol m −2 s −1 ).Leaves were harvested at 0, 1, 4, and 6 days after the start of each treatment and used for RT-qPCR.
We selected four anthocyanin-related genes (BsF3H, BsCHS, BsPAL, and BsANS) and designed primers based on these genes (see Supplementary Materials-Table S1).Gene expression was measured under different conditions with the StepOne and StepOnePlus Real-Time PCR Systems (Applied Biosystems, Foster City, CA, USA) using SYBR ® Premix Ex Taq™ II (TaKaRa, Dalian, China).We used the B. semperflorens 18S rRNA gene (GenBank Accession No. KJ959633) as an internal reference to identify differences among each cDNA template.The RT-qPCR cycling conditions were as follows: 95 • C for 10 s followed by 40 cycles of 95 • C for 10 s and 60 • C for 45 s.The relative transcript level of each target gene was estimated as the transcript ratio of 18S rRNA gene using the 2 −∆∆CT method [37] based on three biological replicates.Each RT-qPCR was repeated three times.The data and figures were generated using SPSS (v22.0,SPSS Inc., Chicago, IL, USA) and Origin (v8.0, OriginLab, Northampton, MA, USA) software.Means were analysed using a one-way analysis of variance implemented in DPS (v3.01,Ruifeng Information Technology Co., Ltd., Hangzhou, China).Multiple comparisons between means were conducted using Tukey's test at the p < 0.05 level.

Anthocyanin Quantification
The anthocyanin contents were determined according to the methods described by Mita et al. [38].Frozen leaf samples (0.3 g) were homogenized with 3 mL of 1% (v/v) HCl: MeOH and then maintained under dark conditions at 4 • C for one day.After centrifugation at 3500 g for 15 min, the supernatant was measured spectrophotometrically at 530 and 657 nm.One unit of anthocyanins equals one absorbance unit (A530-0.25 × A657) per ml of extraction solution.Means were analysed using a one-way analysis of variance implemented in DPS 3.01.Multiple comparisons between means were conducted using Tukey's test at the p < 0.05 level.

RNA-Seq and Sequence Assembly
Leaves of B. semperflorens from the three treatment groups (i.e., CK, LT, and HL), with different leaf anthocyanin levels, were employed for library construction using Illumina sequencing technology (Figure 1).

RNA-Seq and Sequence Assembly
Leaves of B. semperflorens from the three treatment groups (i.e., CK, LT, and HL), with different leaf anthocyanin levels, were employed for library construction using Illumina sequencing technology (Figure 1).After quality checks, low-quality data were filtered out, resulting in 25,022,374, 23,025,183, and 22,506,321 clean read pairs in the LT, HL, and CK libraries, respectively.The Phred quality scores (>30, Q30) for the three samples were 92.10%, 94.16%, and 94.07%, respectively.The GC contents of the LT, HL, and CK libraries were 46.81%, 47.16%, and 47.49%, respectively.
Short clean read pairs were then assembled into 14,221,515 contigs in the total library, and these were further assembled into 205,601 transcripts and 94,880 unigenes.The calculating of N50 is as follows (take Contig N50, for example): Firstly, different lengths of contigs were achieved after splicing the read pairs.Adding all the lengths of Contigs, the total length of the Contigs could be obtained.Then all Contigs were sorted from long to short, such as Contig 1, Contig 2, Contig 3 … Contig 25.The Contigs were added in this order.When the length of the addition reached half of the total length of the Contigs, the last added length of the Contig was Contig N50.The N50 values of the transcripts and unigenes were 2286 bp and 1064 bp, respectively.The mean lengths of the transcripts and unigenes were 1297.91 bp and 635.25 bp, respectively (Table 1).After quality checks, low-quality data were filtered out, resulting in 25,022,374, 23,025,183, and 22,506,321 clean read pairs in the LT, HL, and CK libraries, respectively.The Phred quality scores (>30, Q30) for the three samples were 92.10%, 94.16%, and 94.07%, respectively.The GC contents of the LT, HL, and CK libraries were 46.81%, 47.16%, and 47.49%, respectively.
Short clean read pairs were then assembled into 14,221,515 contigs in the total library, and these were further assembled into 205,601 transcripts and 94,880 unigenes.The calculating of N50 is as follows (take Contig N50, for example): Firstly, different lengths of contigs were achieved after splicing the read pairs.Adding all the lengths of Contigs, the total length of the Contigs could be obtained.Then all Contigs were sorted from long to short, such as Contig 1, Contig 2, Contig 3 . . .Contig 25.The Contigs were added in this order.When the length of the addition reached half of the total length of the Contigs, the last added length of the Contig was Contig N50.The N50 values of the transcripts and unigenes were 2286 bp and 1064 bp, respectively.The mean lengths of the transcripts and unigenes were 1297.91 bp and 635.25 bp, respectively (Table 1).
Another CK seedling was used as a biological replicate with the former one.RNA-seq was done and the correlations of the RNA-seq data were analysed.The duplicate showed linear correlation with its corresponding sample.The r 2 value (r represents Pearson r values) was 0.9779.

Transcriptome Functional Annotation
Because genomes closely related to B. semperflorens were not available, the B. semperflorens unigenes were annotated as a transcriptome.All of the unigenes were aligned to the Nr protein, Swiss-Prot, KEGG, COG, and GO protein databases.A total of 42,027 unigenes were identified and annotated, with 41,161, 28,154, 17,165, 12,577, and 25,415 unigenes being identified in the Nr protein, Swiss-Prot, KEGG, COG, and GO protein databases, respectively (Table 2).

Differentially Expressed Genes
DEGs were screened using the criteria of a false discovery rate (FDR) ≤ 0.01 and an absolute value of the log 2 ratio ≥ 2. A total of 4643 DEGs (4.89% of unigenes) were identified in the comparison of LT vs. CK, and these included 2323 upregulated and 2320 downregulated genes.A total of 757 DEGs (0.80% of unigenes) were identified in the comparison of HL vs. CK; among these, 415 were upregulated, and 342 were downregulated (see Supplementary Materials-Table S2).In total, 4156 of the 4643 DEGs from the LT vs. CK comparison were identified and annotated in the five protein databases, whereas 663 of the 757 DEGs from the HL vs. CK comparison were identified and annotated in the five protein databases (see Supplementary Materials-Table S3).Common DEGs in the comparisons of HL vs. CK and LT vs. CK were retrieved and then analysed against the GO, COG, and KEGG databases.We annotated 1353 DEGs using the GO classification and assigned them to 36 subgroups (Figure 2).The results of the COG category assignment are shown in Figure 3.

Pathway Enrichment Analysis
To analyse the effects of LT and HL on B. semperflorens, we employed the unigenes in BLAST searches of the five protein libraries.DEGs were screened using the criterion of a false discovery rate (FDR) ≤ 0.05.
As shown in Table 3, LT stress significantly (p ≤ 0.05) altered 27 metabolism pathways; the categories that contained the most DEGs were "carbon metabolism", "biosynthesis of amino acids", and "phenylpropanoid biosynthesis".HL stress significantly (p < 0.05) altered only 10 metabolic pathways, with "protein processing in endoplasmic reticulum", "ribosome", and "flavonoid biosynthesis" containing the greatest numbers of DEGs.These data suggest that HL stress promotes flavonoid biosynthesis through protein processing, whereas LT stress redirects carbon metabolism to phenylpropanoid biosynthesis through the biosynthesis of amino acids in Begonia seedlings.Although the LT and HL stresses had similar effects on anthocyanin biosynthesis, these likely occur through different mechanisms.The DEGs that were altered by both LT and HL stresses are annotated in Table 3.We identified 1017 DEGs involved in 28 pathways.However, the common DEGs, which comprised 89 DEGs in nine pathways, were almost identical to the DEGs regulated by HL stress alone, whereas 928 DEGs in 18 pathways were regulated specifically by LT stress.The categories "protein processing in endoplasmic reticulum," "ribosome", and "flavonoid biosynthesis" contained the highest numbers of DEGs (see Supplementary Materials-Table S4).
The above data suggested that the induction of anthocyanin biosynthesis under LT stress was more complex than that under HL stress.
The DEGs regulated by both LT and HL stresses are shown in Figure 4.All of the common DEGs, including PAL, CH4, 4CL, HCT, CAD, CHS, F3H, F3'H, DFR, ANS, 3GT, and GST, were upregulated by both LT and HL, although a difference in expression level of each common DEG was found between the LT and HL conditions.We therefore infer that the common products of phenylpropanoid and flavonoid biosynthesis regulated by both LT and HL include lignins and anthocyanins.
We also analysed the DEGs that were regulated differently in phenylpropanoid and flavonoid biosynthesis by LT and HL.As shown in Figure 5, the LT condition has complex effects on expression of genes involved in phenylpropanoid and flavonoid biosynthesis.For example, both the LT and HL conditions upregulated expression in some numbers of 4CL and CAD (Figure 4), other number(s) of 4CL and CAD were downregulated only by LT.The up-and downregulation in expression were also observed on COMT and PER, which were regulated only by LT; CCoAMOT and C3′H were upregulated by only LT and HL, respectively (Figure 5).Additionally, the LT condition upregulated the expression of CHS, CHI, F3H, and F3′H, genes belonging to the genes upstream of flavonoid biosynthesis.We therefore infer that the LT condition might have enhanced the biosynthesis of other products upstream of phenylpropanoid biosynthesis (such as chalcone and flavanone biosynthesis) in Begonia seedlings.Although the LT condition upregulated one of the ANSs (Figure 4), it downregulated another ANS that was not expressed under the HL condition (Figure 5).Interestingly, Anthocyanins are one type of the major products of phenylpropanoid biosynthesis.In the present study, F3 5 H was not identified in HL-and LT-treated seedlings (data not shown), whereas F3 H and DFR were identified in both HL-and LT-treated seedlings.
The DEGs regulated by both LT and HL stresses are shown in Figure 4.All of the common DEGs, including PAL, CH4, 4CL, HCT, CAD, CHS, F3H, F3'H, DFR, ANS, 3GT, and GST, were upregulated by both LT and HL, although a difference in expression level of each common DEG was found between the LT and HL conditions.We therefore infer that the common products of phenylpropanoid and flavonoid biosynthesis regulated by both LT and HL include lignins and anthocyanins.
We also analysed the DEGs that were regulated differently in phenylpropanoid and flavonoid biosynthesis by LT and HL.As shown in Figure 5, the LT condition has complex effects on expression of genes involved in phenylpropanoid and flavonoid biosynthesis.For example, both the LT and HL conditions upregulated expression in some numbers of 4CL and CAD (Figure 4), other number(s) of 4CL and CAD were downregulated only by LT.The up-and downregulation in expression were also observed on COMT and PER, which were regulated only by LT; CCoAMOT and C3 H were upregulated by only LT and HL, respectively (Figure 5).Additionally, the LT condition upregulated the expression of CHS, CHI, F3H, and F3 H, genes belonging to the genes upstream of flavonoid biosynthesis.We therefore infer that the LT condition might have enhanced the biosynthesis of other products upstream of phenylpropanoid biosynthesis (such as chalcone and flavanone biosynthesis) in Begonia seedlings.Although the LT condition upregulated one of the ANSs (Figure 4), it downregulated another ANS that was not expressed under the HL condition (Figure 5).Interestingly, LT and HL had different effects on proanthocyanidin (PA) biosynthesis.As monomers of PAs, catechin and epicatechin are biosynthesized via LAR and ANR, respectively [39].As shown in Figure 5, the LT condition upregulated both LAR and ANR, whereas the HL condition had no effect on LAR and downregulated ANR.Although anthocyanins were biosynthesized in both the LT-and HL-treated seedlings (Figure 1), the LT-treated seedlings might have produced more PAs than the HL seedlings.
LT and HL had different effects on proanthocyanidin (PA) biosynthesis.As monomers of PAs, catechin and epicatechin are biosynthesized via LAR and ANR, respectively [39].As shown in Figure 5, the LT condition upregulated both LAR and ANR, whereas the HL condition had no effect on LAR and downregulated ANR.Although anthocyanins were biosynthesized in both the LT-and HLtreated seedlings (Figure 1), the LT-treated seedlings might have produced more PAs than the HL seedlings.Lignin is an important compound produced by the phenylpropanoid biosynthesis pathway [40], and lignin biosynthesis occurs in response to many abiotic and biotic stresses [41].
Light itself has a positive effect on lignin biosynthesis [42], and high light might increase the lignin content [43] via increased expression of lignin-biosynthesis-genes (such as PAL, CAD, and PER) [40].In our research, the HL condition upregulated CAD and PER, which is expected to increase the lignin content; the expression of HCT and C3′H were also enhanced by the HL condition, which might indicate a change in the lignin composition from H lignin to G lignin.Lignin deposition in vascular bundles has been found to improve stress tolerance in Arabidopsis thaliana (L.) Heynh [44].
The effects of low temperature on lignin biosynthesis are complex.Low temperature (10 °C) increased lignin content in stored loquat fruit [45] and Increased the lignin content of ex vitro, but not in vitro, poplar seedlings [46].Annual average temperature is positively correlated with lignin content in Picea abies (L.) Karst.[47].Cold acclimation and freezing treatment did not modify the lignin content, but increased pectins content in rape leaves [48].The relationship between the activities of enzymes involved in lignin biosynthesis and lignin content in plants treated with low temperature is unclear [41].In our research, the LT condition had different, complex effects on genes involved in lignin biosynthesis.In addition to the genes upstream of lignin biosynthesis (for example, PAL and C4H), HCT, C3′H, and CCoAOMT were upregulated in Begonia seedlings under the LT and Lignin is an important compound produced by the phenylpropanoid biosynthesis pathway [40], and lignin biosynthesis occurs in response to many abiotic and biotic stresses [41].
Light itself has a positive effect on lignin biosynthesis [42], and high light might increase the lignin content [43] via increased expression of lignin-biosynthesis-genes (such as PAL, CAD, and PER) [40].In our research, the HL condition upregulated CAD and PER, which is expected to increase the lignin content; the expression of HCT and C3 H were also enhanced by the HL condition, which might indicate a change in the lignin composition from H lignin to G lignin.Lignin deposition in vascular bundles has been found to improve stress tolerance in Arabidopsis thaliana (L.) Heynh [44].
The effects of low temperature on lignin biosynthesis are complex.Low temperature (10 • C) increased lignin content in stored loquat fruit [45] and Increased the lignin content of ex vitro, but not in vitro, poplar seedlings [46].Annual average temperature is positively correlated with lignin content in Picea abies (L.) Karst.[47].Cold acclimation and freezing treatment did not modify the lignin content, but increased pectins content in rape leaves [48].The relationship between the activities of enzymes involved in lignin biosynthesis and lignin content in plants treated with low temperature is unclear [41].In our research, the LT condition had different, complex effects on genes involved in lignin biosynthesis.In addition to the genes upstream of lignin biosynthesis (for example, PAL and C4H), HCT, C3 H, and CCoAOMT were upregulated in Begonia seedlings under the LT and HL conditions, which might indicate a change in the lignin composition from H lignin to G and S lignin.
In addition, both up-and downregulation of 4CL, COMT, CAD, and PER were observed under the LT condition.LT condition therefore has more complex effects on lignin accumulation than HL condition in B. semperflorens.Despite the upregulation of genes involved in anthocyanin biosynthesis, and the increase of anthocyanin content under LT conditions (Figures 1 and 4), accumulation of anthocyanins is not associated with increase of lignin content.There may be different manners in which lignin and anthocyanin inductions occur under low temperature [48].Begonia seedlings biosynthesised anthocyanins during flavonoid biosynthesis, suggesting a possible role of anthocyanin accumulation in leaves under both LT and HL stress.Anthocyanins can perform photoprotective and antioxidative roles under many environmental stimuli [49,50].Although LT and HL stresses rarely occur simultaneously in nature, both are important inductors of anthocyanin biosynthesis in many plants.Anthocyanins were shown to resist LT and HL stresses in B. semperflorens leaves [51,52], which might explain why anthocyanins accumulated during flavonoid biosynthesis under both the LT and HL conditions in the present study.
PAs and anthocyanins are derived from the common flavonoid biosynthetic pathway and share the same precursors (leucoanthocyanidins), or are oxidized to anthocyanidins (Figure 4), suggesting a competition between PAs and anthocyanin biosynthesis.With the promotion of PAR, PAs might be produced by repressing IFS or redirecting anthocyanin flux to PA biosynthesis pathway in legume [53].Although anthocyanins and PAs share part of the same biosynthesis pathway, they have different synthesis patterns [54].In berries, PAs accumulate at the initial stage, whereas anthocyanins accumulate only during ripening [54].
In our research, the HL condition did not affect the expression of LAR, but it decreased ANR expression.We therefore infer that anthocyanin biosynthesis is the dominant flavonoid pathway.The LT condition upregulated LAR and ANR, suggesting it enhanced both anthocyanin and PA biosynthesis.
Similar to anthocyanidins, PAs are important products of the flavonoid pathway [53].Anthocyanidins and PAs are both involved in plant resistance and have significant potential to scavenge free radicals in both in vitro and in vivo models [55].However, anthocyanins can attenuate light to protect plants from photoinhibition.When plants are subjected to excessive light, decreasing light absorption becomes the immediate and cardinal task.This phenomenon might explain why the HL condition primarily upregulated genes involved in anthocyanin biosynthesis rather than those involved in PA biosynthesis in Begonia seedlings (Figure 5).Under the LT condition, plants first suffered the stress of low temperature, which will lead to oxidation from many other organelles besides chloroplasts.The Begonia seedlings therefore upregulated genes involved in both anthocyanin and PA biosynthesis under the LT condition (Figure 5).

RT-qPCR Analysis of Four Core Genes Involved in Anthocyanin Biosynthesis
Due to the feedback from upregulated DEGs involved in anthocyanin biosynthesis, both LT and HL stresses induced significant anthocyanin accumulation in Begonia seedlings (Figure 1).
To verify the validity of the anthocyanin-biosynthesized genes we got from the data of transcriptomes, we selected four core anthocyanin biosynthesis genes (BsPAL, BsCHS, BsF3H, and BsANS) from the DEGs for RT-qPCR analysis of their expression under LT and different light environments.
As shown in Figure 6, the same light-intensity (300 µmol m -2 s -1 , CC) and low-temperature (15 • C/5 • C, LC) conditions resulted in significant upregulation of the expression of the four anthocyanin biosynthesis genes after 24 h (day 1).We also found that a significant quantity of anthocyanins accumulated during the second day (day 2).Low temperature and low light (100 µmol m -2 s -1 LD) did not result in upregulation of the four genes and did not increase anthocyanin accumulation.In contrast, low temperature and high light (500 µmol m -2 s -1 , LH) significantly induced the expression of the four genes, and the resulting increased expression was higher than that observed in seedlings grown under the LC conditions.These results indicate that both high-light and low-temperature conditions promote the expression of these four core genes involved in anthocyanin biosynthesis.Elevated anthocyanin concentrations were observed in leaves exposed to LC and LH conditions (Figure 7), which is consistent with the transcript levels of these four core genes.These results indicate that these four core genes are related to anthocyanin biosynthesis.Elevated anthocyanin concentrations were observed in leaves exposed to LC and LH conditions (Figure 7), which is consistent with the transcript levels of these four core genes.These results indicate that these four core genes are related to anthocyanin biosynthesis.Elevated anthocyanin concentrations were observed in leaves exposed to LC and LH conditions (Figure 7), which is consistent with the transcript levels of these four core genes.These results indicate that these four core genes are related to anthocyanin biosynthesis.

Conclusions
In the present study, we employed RNA-Seq to study the transcriptome of anthocyanin biosynthesis and transport in B. semperflorens.By identifying DEGs regulated by both the LT and HL conditions, we found that anthocyanin accumulation was employed as a common strategy by Begonia seedlings in resisting LT and HL stress.By identifying DEGs regulated differently by the LT and HL conditions, we found that Begonia seedlings also had some different strategies for resisting LT and HL stresses: Begonia seedlings mainly biosynthesized anthocyanins under the HL condition, under which photoprotection is more necessary for plants than under LT condition; while Begonia seedlings mainly biosynthesized lignin, proanthocyanidins, and anthocyanins under the LT condition.

Figure 1 .
Figure 1.Phenotypes (a) and anthocyanin contents (b) in the leaves of B. semperflorens under different environmental conditions.Means denoted by the same letter (a, b, c) did not significantly differ at p < 0.05 according to Tukey's test.LT, low temperature-treated plants; HL, high light-treated plants; CK, control plants; FW: fresh weight.

Figure 1 .
Figure 1.Phenotypes (a) and anthocyanin contents (b) in the leaves of B. semperflorens under different environmental conditions.Means denoted by the same letter (a, b, c) did not significantly differ at p < 0.05 according to Tukey's test.LT, low temperature-treated plants; HL, high light-treated plants; CK, control plants; FW: fresh weight.

Figure 2 .
Figure 2. Gene Ontology (GO) classification of common differentially expressed genes (DEGs) that were detected in both comparisons of LT vs. CK and high light HL vs. CK in B. semperflorens leaves.LT, low temperature-treated plants; HL, high light-treated plants; CK, control plants.

Figure 3 .
Figure 3. Clusters of orthologous groups (COG) classification of common differentially expressed genes (DEGs) that were detected in both comparisons of LT vs. CK and HL vs. CK in B. semperflorens leaves.LT, low temperature-treated plants; HL, high light-treated plants; CK, control plants.

Figure 2 .
Figure 2. Gene Ontology (GO) classification of common differentially expressed genes (DEGs) that were detected in both comparisons of LT vs. CK and high light HL vs. CK in B. semperflorens leaves.LT, low temperature-treated plants; HL, high light-treated plants; CK, control plants.

Figure 2 .
Figure 2. Gene Ontology (GO) classification of common differentially expressed genes (DEGs) that were detected in both comparisons of LT vs. CK and high light HL vs. CK in B. semperflorens leaves.LT, low temperature-treated plants; HL, high light-treated plants; CK, control plants.

Figure 3 .
Figure 3. Clusters of orthologous groups (COG) classification of common differentially expressed genes (DEGs) that were detected in both comparisons of LT vs. CK and HL vs. CK in B. semperflorens leaves.LT, low temperature-treated plants; HL, high light-treated plants; CK, control plants.

Figure 3 .
Figure 3. Clusters of orthologous groups (COG) classification of common differentially expressed genes (DEGs) that were detected in both comparisons of LT vs. CK and HL vs. CK in B. semperflorens leaves.LT, low temperature-treated plants; HL, high light-treated plants; CK, control plants.

Figure 5 .
Figure 5. Schematic of physiological and the expression of differentially expressed genes (DEGs) regulated only by LT or HL condition in phenylpropanoid and flavonoid biosynthesis in B. semperflorens.

Figure 5 .
Figure 5. Schematic of physiological and the expression of differentially expressed genes (DEGs) regulated only by LT or HL condition in phenylpropanoid and flavonoid biosynthesis in B. semperflorens.

Table 1 .
Summary of the sequence assembly.

Table 1 .
Summary of the sequence assembly.

Table 2 .
Annotation of the unigenes against different databases.COG, Clusters of Orthologous Groups; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; KOG, EuKaryotic Orthologous Groups; Pfam, Protein Families; Nr, Non-Redundant protein.

Table 3 .
Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis of differentially expressed genes (DEGs) detected in the comparison of LT vs. CK and HL vs. CK in B. semperflorens leaves.LT, low temperature-treated plants; HL, high light-treated plants; CK, control plants.