Transcriptome Analysis Reveals Molecular Signatures of Luteoloside Accumulation in Senescing Leaves of Lonicera macranthoides

Lonicera macranthoides is an important medicinal plant widely used in traditional Chinese medicine. Luteoloside is a critical bioactive compound in L. macranthoides. To date, the molecular mechanisms underlying luteoloside biosynthesis are still largely unknown. In this work, high performance liquid chromatography (HPLC) was employed to determine the luteoloside contents in leaves, stems, and flowers at different developmental stages. Results showed that senescing leaves can accumulate large amounts of luteoloside, extremely higher than that in young and semi-lignified leaves and other tissues. RNA-Seq analysis identified that twenty-four differentially expressed unigenes (DEGs) associated with luteoloside biosynthesis were significantly up-regulated in senescing leaves, which are positively correlated with luteoloside accumulation. These DEGs include phenylalanine ammonia lyase 2, cinnamate 4-hydroxylase 2, thirteen 4-coumarate-CoA ligases, chalcone synthase 2, six flavonoid 3′-monooxygenase (F3′H) and two flavone 7-O-β-glucosyltransferase (UFGT) genes. Further analysis demonstrated that two F3′Hs (CL11828.Contig1 and CL11828.Contig2) and two UFGTs (Unigene2918 and Unigene97915) might play vital roles in luteoloside generation. Furthermore, several transcription factors (TFs) related to flavonoid biosynthesis including MYB, bHLH and WD40, were differentially expressed during leaf senescence. Among these TFs, MYB12, MYB75, bHLH113 and TTG1 were considered to be key factors involved in the regulation of luteoloside biosynthesis. These findings provide insights for elucidating the molecular signatures of luteoloside accumulation in L. macranthoides.


Introduction
Lonicera macranthoides, a member of the Caprifoliaceae family, is a medicinal plant primarily distributed in Southern China. It has been widely used as a critical raw material in traditional Chinese medicine for thousands of years, because it can effectively treat H1N1, respiratory syndrome, and hand-foot-and-mouth disease (China Pharmacopeia Commission, 2010). Furthermore, L. macranthoides extracts are utilized as important and indispensable ingredients in functional foods, beverages, wine and cosmetics [1,2]. A subset of secondary metabolites, including flavonoids, phenolic acid, volatile oil and saponins, is present in whole plants of L. macranthoides [3]. These components possess a wide range of pharmacological activities, such as antioxidant, anti-inflammatory, antibacterial, In the present study, luteoloside in different tissues and developmental stages was determined using high performance liquid chromatography (HPLC) and transcriptome-wide sequencing in L. macranthoides was performed using the Illumina HiSeq™ 2000 platform to explore the key genes encoding biosynthetic enzymes and regulator proteins associated with the luteoloside biosynthesis in L. macranthoides. These results provide genomic resources for future enhancement of luteoloside production by genetic strategies in L. macranthoides.

Luteoloside Contents in Different Tissues and Developmental Stages of L. macranthoides
Luteoloside contents in the leaves, stems and flowers were quantified by HPLC ( Figure S1). The contents differed significantly among the three tissues and different developmental stages ( Figure 1). Senescing leaves (SL) can accumulate very large amounts of luteoloside (up to 1.625 mg/g DW), which was approximately 35-and 160-fold higher than that in young leaves (YL) and semi-lignified leaves, respectively ( Figure 1A,B). There was no obvious change in luteoloside accumulation among the various developmental stages of stems ( Figure 1C,D). Luteoloside contents fluctuated throughout the flower development and the highest luteoloside level in flowers was found in white flower, followed by green flower buds in length of 10 mm ( Figure 1E,F). For all tissues examined, flowers displayed remarkably higher luteoloside contents than that in leaves and stems of various developmental stages except for senescing leaves. Values are means ± SD (n = 3). Duncan's multiple range test was used to analyze the significance and the different lower-case letters (e.g., a and b) indicate significant (p < 0.05) differences between samples.

RNA-Seq Analysis
To further explore the candidate genes involved in luteoloside biosynthesis of L. macranthoides, two libraries including SL and YL were established for RNA-Seq analysis. The raw data generated by sequencing for each library ranged from 82.4 to 87.5 million reads. After filtering, a total of 76.8 to 82 million clean reads were obtained, accounting for 93% of the raw reads, with an average GC % of 44.19% for all clean reads being (Supplementary Table S1).
These clean reads were assembled into 260,079 and 173,214 contigs for SL and YL libraries, respectively. After processing by using Trinity and TGICL software, the contigs of SL were further incorporated into 183,667 unigenes with an average length of 903 bps and an N50 length of 1501 bps. Moreover, a total of 122,824 unigenes were obtained in YL library, with an average length of 1131 bp and an N50 value of 1859 bps ( Table 1). The length distribution for all unigenes in L. macranthoides is shown in Figure S2. The numbers of unigenes with sequence length longer than 500 bp, 1000 bp, 2000 bp and 3000 bp accounted for 70.05% (111,272), 41.3% (65,614), 16.8% (26,673) and 5.7% (9053) of the total unigenes, respectively ( Figure S2). The assembled data for L. macranthoides in this study was improved compared to our earlier published transcriptome assemblies, with remarkable increases in overall unigene length and N50 value [11].

Differentially Expressed Unigenes (DEGs) in Senescing and Young Leaves
To improve the accuracy of the RNA-Seq data, two independent biological replicates were performed for each sample in this study. The Pearson and Spearman correlation coefficients were calculated to investigate the correlation of the gene expression data among the biological replicates and visualized with a color-coded diagram. The Pearson correlation coefficients were 0.99 and 0.83 between the two libraries of SL and YL, respectively. And the Spearman correlation coefficients were 0.79 and 0.95 ( Figure S4A,B). These results indicate strong correlation between variables and reliability of the RNA-Seq data.
The DEGs between senescing leaves and young leaves were identified (Supplementary Table S3), and the distribution of log-fold changes was visualized in the volcano plot ( Figure S4C). With respect to young leaves, a total of 17,020 unigenes were differentially expressed in senescing leaves, of which 13,662 unigenes were up-regulated and 3,358 unigenes were down-regulated ( Figure 2).

GO-and KEGG-Based Functional Classification of DEGs
To gain an insight into the functional characterization of DEGs, GO-based classification was performed by using the Blast2GO program v 3.0 [32]. Significantly enriched GO categories were identified (p < 0.05), 16 for cellular component ( Figure 3A), 9 for biological process ( Figure 3B), and 7 for molecular function ( Figure 3C). Several biological processes, such as cuticle development, epidermal and epithelial cell differentiation, long-chain fatty acid metabolic process and wax biosynthetic process, were differentially regulated at various developmental stages of leaves ( Figure 3B). Moreover, for the molecular function category, GO terms corresponding to oxidoreductase activity, hydrolase activity and monooxygenase activity were remarkably enriched ( Figure 3C). This result indicates the potential roles of these GO terms in luteoloside biosynthesis.
To better understand the gene function, we assigned 10,216 of the DEGs to one of 128 KEGG terms by using the KEGG database. Significantly altered biological pathways were identified during leaf senescence (Q value ≤ 0.05). Several metabolic pathways, including biosynthesis of secondary metabolites, phenylpropanoid biosynthesis, lipid and fatty acid metabolism, starch and sucrose metabolism, tryptophan metabolism, cutin and wax biosynthesis, and flavone and flavonol biosynthesis displayed outstanding enrichments ( Figure S5). This finding indicates a significant difference in primary and second metabolism between young and senescing leaves in L. macranthoides. A total of 339 and 67 DEGs were assigned to phenylpropanoid metabolic and flavone/flavonol biosynthetic pathways, respectively ( Figure S5), by which generate precursors for luteoloside biosynthesis.

Identification of Candidate Genes Associated with Luteoloside Biosynthetic Pathways
To elucidate the key genes for luteoloside biosynthesis, DEGs involved in phenylpropanoid pathway were screened, including three PAL, one C4H, fourteen 4CL, eight CHS, two CHI, one FNS, eight F3 H and seven UFGT genes ( Figure 4). PAL2 (Unigene109136), C4H2 (CL11118.Contig2) and the majority of 4CL genes (13 out of 14) were significantly upregulated in senescing leaves compared with young leaves. We also found that CHS2 (CL19869.Contig1), six F3 H (CL11828.Contig1, CL11828.Contig2, CL7653.Contig1, Unigene102655, Unigene65437 and Unigene76746), two UFGT (Unigene2918 and Unigene97915) genes displayed dramatic increases in transcript levels during leaf senescence. Meanwhile, CHI and FNS paralogous genes appeared remarkable down-regulation in senescing leaves, compared with those in young leaves ( Table 2). The mRNA levels of ten unigenes possibly involved in luteoloside biosynthesis were validated by qRT-PCR, which displayed similar expression patterns to the results obtained by RNA-Seq ( Figure 5A-J).These DEGs, especially those functioning downstream in the metabolic pathway, may provide valuable clues to illustrate the luteoloside biosynthetic pathway. The high expression levels of CHS2, F3 H and UFGT genes in senescing leaves were consistent with the large amounts of luteoloside, indicating their fundamental involvement in luteoloside biosynthesis. Proposed pathways for luteoloside biosynthesis in Lonicera macranthoides were illustrated by RNA-Seq analysis. Luteolin, the precursor of luteoloside, is biosynthesized from the general flavonoid precursor: naringenin. Circle 1( 1 ) indicates that luteolin is biosynthesized directly from apigenin catalyzed by F3'H. Circle 2( 2 ) indicates that luteolin is generated directly from eriodictyol catalyzed by FNS. Expression profile for each gene was shown in colored blocks and each blocks represented the expression changes (represented by Log2Ratio) in senescing leaves with respect to young leaves. Red colors/green colors correspond to up-/down-regulation of these genes and Log2Ratio ≥ 1 is considered statistically significant. Details were showed in Table 2. Abbreviations: PAL, phenylalanine ammonia lyase; C4H, cinnamate 4-hydroxylase; 4CL, 4-hydroxycinnamoyl CoA ligase/4-coumarate-CoA ligase; CHS, chalcone synthase; CHI, chalcone isomerase; FNS, flavone synthase; F3H, flavonoid 3 -monooxygenase/flavonoid 3 -hydroxylase; UF7GT, flavone 7-O-β-glucosyltransferase. "Inf" and "−Inf" indicate that the Log2Ratio of SL to YL is infinity and negetive infinity, respectively.
To further characterize the molecular properties of the F3 H and UFGT genes, phylogenetic analyses of protein sequences were performed. Two F3 H genes, namely CL11828.Contig1 and CL11828.Contig2, were grouped into the same clade as the F3 H genes in Sesamum indicum, Solanum lycopersicum, Ziziphus jujube and Nicotiana tabacum and also shared high sequence homology to the F3 H genes in Arabidopsis lyrata and Oryza sativa. This result suggests the similar functions of these F3 H genes. However, the four remaining F3 H genes showed relatively lower levels of similarity to those genes in other species ( Figure S6A). Two UFGTs were classified into the same subcluster with UDP-glucose flavonoid 3-O-glucosyltransferase (UF3GTs) in several species and also displayed amino acid sequence similarities with UF7GTs in Arabidopsis and Scutellaria ( Figure S6B).

Identification of DEGs Involved in MYB, bHLH and WD40 Transcription Factors
MYB, bHLH, and WD40 have been previously identified to regulate flavonoids biosynthesis [30]. RNA-seq data showed that thirty MYB, thirty-one bHLH and eighty-six WD40 TFs in senescing leaves displayed differential expression with respect to young leaves (Table 3 and Supplementary  Table S4). In contrast to young leaves, approximately half of the differentially expressed MYBs (16 out of 30) underwent significant down-regulation in senescing leaves, such as those orthologous to TT2, AtMYB12, PsMYB1, SlMYB46, VvMYB44, and VvMYB306. Meanwhile, the remaining differentially expressed MYBs, including those homologous to AtMYB48, VvMYB1R1, VvMYB75, AmMYB305, EsMYB13, and GhMYB9, showed dramatically increased transcripts in senescing leaves (Table 3). Moreover, out of the 30 bHLHs, 14 bHLHs, which include bHLH25, bHLH49, bHLH66, bHLH78, bHLH113, bHLH123 and bHLH143, were strongly up-regulated during leaf senescence (Table 3). In addition, among the DEGs involved in WD40 repeat-containing proteins, one transparent testa glabrous1 (TTG1) ortholog was remarkably upregulated in senescing leaves compared with that in young leaves (Supplementary Table S4). This gene contained five WD repeats and showed 37% amino acid identity with TTG1 in Arabidopsis (XP_020877597), 38% identity with TTG in Medicago (XP_003624554) and 37% identity with TTG1 in Malus (AHM88209). QRT-PCR analysis also showed that the transcripts of several unigenes, including those homologous to MYB, bHLH and TTG1, were in dramatic regulation during leaf senescence ( Figure 6A-I).

Discussion
Luteoloside belongs to a group of natural flavonoids isolated from Lonicera, that exerts human health benefits, including antiallergic, free radical scavenging and antioxidant, antihyperglycemic and anticholestatic activities [9]. Luteoloside accumulation is a precisely regulated process that varies considerably depending on plant species, developmental stages and different tissues. To date, fingerprint analysis of luteoloside has been carried out in Lonicera. The amounts of luteoloside are remarkably lower in flower buds of L. hypoglauca than inthose of L. japonica, whereas no obvious change is observed between the flower buds of L. japonica and L. macranthoides [8]. However, according to the description in the Chinese pharmacopoeia 2010 and the view reported by Wu et al. [7], luteolin and luteoloside are in lower abundance in flower buds of L. macranthoides than in those of L. japonica [33]. In the present study, luteoloside accumulation initially increased from stage 2 to stage 4 and then decreased in stage 5 during flower development in L. macranthoides ( Figure 1F), which is in accordance with the results of previous reports in L. japonica [1,7]. Moreover, in contrast to the findings of Yuan et al. in L. japonica [8], our results revealed that luteoloside displayed significantly higher level in flowers than that in stems and leaves at different growth stages except for senescing leaves (Figure 1). In general, luteoloside is relatively low in L. macranthoides (only 0.01 to 0.2 mg/g). Our study demonstrated that luteoloside was accumulated at extremely higher level in senescing leaves of L. macranthoides ( Figure 1B). This finding provides us with the view of significant medical value of the senescing leaves and new insights into the molecular basis underlying luteoloside biosynthesis.

Molecular Features Underlying Abundant Accumulation of Luteoloside during Leaf Senescence in L. macranthoides
Transcriptome analysis in leaves at different developmental stages was aimed at elucidating molecular mechanism underlying high accumulation of luteoloside in senescing leaves of L. macranthoides. Several genes encoding PAL, C4H and 4CL displayed higher transcript levels in senescing leaves (Figure 4), which might contribute to luteoloside accumulation and thus were considered to be critical members associated with luteoloside generation. These findings are in accordance with the previous observations that the activities and expressions of PAL, C4H and 4CL positively correlate with luteoloside biosynthesis [1,8]. However, the three enzymes are located upstream of the luteoloside metabolic pathway and as common enzymes regulating biosynthesis of various secondary metabolites, such as lignin [34] and flavonoids [35]. Therefore, other genes enabling irreversible commitments to luteoloside biosynthetic pathways need to be further investigated. LjCHS1 and LjCHI2 were proposed to be the important genes participating in luteolin biosynthesis [8], whereas the transcripts of two CHI genes exhibited negative correlation with luteoloside levels in our study (Figure 4), suggesting that the regulation of luteoloside accumulation differs among plant species.
Three enzymes including F3 H, FNS and UFGT, are proved to play crucial roles in luteoloside metabolic pathway [7,10,14,36]. In the present study, the enhanced expression of two F3 H homologs (CL11828.Contig1 and CL11828.Contig2) coincided precisely with the high level of luteoloside in senescing leaves ( Table 2), indicating that they are potential candidates modulating the biosynthesis of luteolin and luteoloside. CYP75B3, a F3 H gene in Oryza sativa, exhibits similar preference for naringenin and apigenin, in nearly catalytic efficiencies for these substrates [36]. The two F3 H genes showed close homology with CYP75B3 ( Figure S6A), suggesting that these genes were not altered to exhibit flavonoid-3-hydroxylation activity and might facilitate the biosynthesis of eriodictyol or luteolin from naringenin or apigenin, respectively. FNS is considered as the critical enzyme in the two routes of luteoloside metabolic pathway (Figure 4), which is implied by the accumulation of luteolin-7-O-glucoside (luteoloside) and apigenin-7-O-glucoside in FNSII-overexpression transgenic lines in Lonicera [7]. Wu et al. [7] also demonstrated that FNSII in L. japonica (LjFNSII-1.1) exhibits an apparently higher catalytic activity (approximate 4-fold) than that in L. macranthoides (LmFNSII-1.1). They concluded that the less abundance of flavones in flowers of L. macranthoides compared with that in L. japonica might be attributed to the weak catalytic efficiency of LmFNSII. However, in our dataset, a FNSII gene (Unigene2335) exhibiting 99% identity with LmFNSII-1.1 (deletion of a leucine at position 22 in LmFNSII-1.1) displayed lower transcript abundance in senescing leaves than in young leaves (Table 2). Thus, it seems not to be responsible for luteoloside generation. In addition, high levels of luteoloside in senescing leaves might result from the highly-expressed UFGTs (Unigene2918 and Unigene97915), homologs of AtUF7GT and SbUF7GT (Figure S6B), which can maximize the luteolin conversion. Hence, the extremely higher amounts of luteoloside in senescing leaves might be attributed to significantly higher transcript levels of F3 H in senescing leaves compared to that in young leaves. The increase in F3 H expression might provide sufficient eriodictyol to FNSII, and then large amounts of luteolin would accumulatedue to the relatively high background expression of FNSII ( Table 2). Given that the FNSII enzymes prefer eriodictyol as a substrate over naringenin in Lonicera [7], the role of FNSII in luteolin generation from the other routes (naringenin-apigenin-luteolin) could be neglected. On the basis of our results and the findings in the literature described above, we hypothesized that luteoloside is preferentially biosynthesized via the following alternative route: naringenin is hydroxylated by F3 H (CL11828.Contig1 and CL11828.Contig2) to produce eriodictyol and then FNSII converts eriodictyol to luteolin, which can be metabolized to luteoloside via UFGTs. Therefore, although the transcript levels of FNSII in young leaves were approximately 7-fold higher than that in senescing leaves (Table 2), lower levels of luteoloside accumulation were observed due to the lower expressions of F3 Hs in young leaves than that in senescing leaves ( Figure 1B).

Transcription Factors as Potential Regulators of Luteoloside Accumulation in Leaves of L. macranthoides
Flavonoid biosynthesis is transcriptionally regulated by a multitude of transcription factors, of which R2R3-MYB, bHLH and WD40 repeat proteins are proven to be vital [37]. Numerous studies have demonstrated that MYB12/11/111 can modulate flavonol synthesis by activating the early biosynthetic genes including PAL, C4H, CHS, CHI and F3 H [16,17,[38][39][40][41]. A recent study in E. sagittatum showed that EsMYBF1, a homolog of AtMYB12 and VvMYBF1, functions as an activator regulating the flavonol synthesis. Ectopic expression of EsMYBF1 resulted in the enhanced accumulation of kaempferol and quercetin via upregulating the expression of CHS, CHI, F3H and FLS genes but a decline in the content of anthocyanin via downregulating the transcripts of DFR and ANS genes, suggesting that EsMYBF1 is a flavonol-specific regulator [18]. Moreover, our results indicated that two putative MYB12 genes in L. macranthoides (Unigene36582 and Unigene52460), which are homologs of AtMYB12, were more abundantly expressed in young leaves (Table 3), in which the phenylpropanoid products, chlorogenic acid and rutin, were largely accumulated (unpublished data) but extremely low content of luteoloside was observed ( Figure 1B). These findings suggest that the two MYB12 homologs might regulate the phenylpropanoid and flavonoid biosynthetic pathways in L. macranthoides. The disparate accumulation patterns of the related main metabolites and derivatives upon both up-and down-regulation of MYB12 might be preferentially ascribed to a flux shift in the metabolic pathway rather than a direct outcome of transcriptional activation or repression in leaves, which is manifested by the reports that overexpression of AtMYB12 in tomato activates the caffeoylquinic acid biosynthetic pathway, while down-regulation of SlMYB12 also leads to the accumulation of caffeic acid derivatives [42]. Moreover, several novel MYBs were observed to correlate with the accumulation pattern of luteoloside during leaf senescence, in particular the positive-correlated genes, such as those homologs to AtMYB48, VvMYB1R1, VvMYB75, AmMYB305, EsMYB13 and GhMYB9 (Table 3), indicating their involvement in the regulation of luteoloside biosynthesis in L. macranthoides leaves.
bHLH, the key component involved in the regulation of flavonoid biosynthesis, was observed to activate late biosynthetic genes (i.e., F3 H, ANR, DFR, and UFGT) through formation of highly dynamic MYB/bHLH/WD40 (MBW) complexes [43][44][45]. In Arabidopsis, three bHLH activators including TRANSPARENT TESTA 8(TT8), enhancer of glabra3 (EGL3), glabra3 (GL3), were reported to act in the transcriptional regulation of anthocyanin and proanthocyanidins (PAs) production via interacting directly or indirectly with Mybs [30,46]. Other bHLH proteins, such as MdbHLH3 in apple and BobHLH1 in cauliflower, interact with MdMYB10 and BoMYB2 respectively, were observed to confer anthocyanin accumulation [21,47]. In the current study, a subset of bHLHs, such as two AtbHLH113 homologs (CL3655.Contig4 and Unigene11884) and four AtbHLH78 homologs, displayed higher mRNA levels in senescing leaves than that in young leaves, which is positively correlated with the accumulation of luteoloside, suggesting their involvements in the regulation of luteoloside generation. This result is in accordance with the findings in Arabidopsis that AtbHLH113 is predicted to interact with PAP1/MYB75 modulating anthocyanin biosynthesis [48].
WD40 repeat protein is another pivotal factor of MBW complexes [46]. In this study, a TTG1 ortholog (Unigene108470) (Supplementary Table S4) in L. macranthoides was isolated and displayed a positive correlation with luteoloside production. In terms of the roles of TTG1 in regulating flavonoid biosynthesis, we considered that Unigene108470 might be the candidate factor in controlling luteoloside biosynthesis. The transcription factors, particularly R2R3 MYB TFs, may activate distinct sets of structural genes of flavonoid biosynthesis [49]. In our dataset, several novel candidate genes participating in the transcriptional regulation of flavonoids biosynthesis in L. macranthoides were identified. Further work will be required to determine whether these genes significantly induce or limit luteoloside biosynthesis in plants and if so, their molecular mechanisms have to be explored via genetic and biochemical approaches in vitro and in vivo.

Plant Materials
Lonicera macranthoides Hand-Mazz (cv. Yu Lei 1#) used in the experiments was planted in a greenhouse of Chongqing University of Arts and Sciences (Chongqing, China). Plant tissues including leaves, stems and flowers at different developmental stages were collected from 5-year-old seedlings. Leaves and stems at three developmental stages including young leaves and stems, semi-lignified leaves and stems, senescing leaves and stems were obtained and illustrated in Figure 1A, and C. Flowers at five developmental stages were harvest in June. These stages are illustrated in Figure 1E.

Determination of Luteoloside Contents by High Performance Liquid Chromatography (HPLC)
Tissue samples were subjected to lyophilization and homogenization with a grinding miller. The powder was then passed through a 40-mesh (420 µm) standard sieve before extraction. A total of 0.25 g homogenized samples was extracted with 20 mL of ethanol (70%, v/v) by ultrasonication for 30 min. The extracts were cooled to room temperature and centrifuged at 4000 rpm for 10 min. Afterward, the supernatant was filtered through a 0.22 µM microfiltration membrane for luteoloside analysis by HPLC.
HPLC was performed on a Shimadzu LC-20A HPLC analyzer (Shimadzu Corporation, Kyoto, Japan), equipped with a LC-20AT pump, a SIL-20A auto sampler, a CBM-20A system controller, a SPD-M20A diode array detector (DAD)and a CTO-20A column oven. Separation was carried out on a Shimadzu Shim-Pack VP-ODS C18 column (5 µm, 250 × 4.6 mm) using gradient elution. Mobile phase A was 2% formic acid in deionized water, while phase B was 2% formic acid in methanol (80:20, v/v). A 15% B linear gradient was conducted from 0.01 min to 8.00 min, followed by a 15~50% B linear gradient from 8.01 min to 25.00 min, and finally a 50% B linear gradient from 25.01 min to 40.00 min, at a flow rate of 0.7 mL/min. The column and the detector were operated at 35 • C. A volume of 20 µL supernatant was injected and the detection wavelength was monitored at 360 nm. Luteoloside standard (>98%) was purchased from SIGMA (Sigma-Aldrich, St. Louis, MO, USA). The luteoloside contents were analyzed in triplicate and calculated based on peak area measurements. Statistical significance was performed with SPSS using Duncan's new multiple range test.

RNA Extraction and Illumina Sequencing
Total RNA from young and senescing leaves was isolated and purified using QIAGEN RNeasy Plant Mini kit and RNase-free DNase set (QIAGEN, Hilden, Germany) according to manufacturer's instruction. RNA-Seq was performed at Hangzhou 1GENE Technology Co., Ltd. (Hangzhou, China). Illumina TruSeq™ RNA Sample Prep Kit (Cat# RS-122-2001) (Illumina, San Diego, CA, USA) was employed to construct cDNA library according to the manufacturer's protocol. Briefly, Poly(A) mRNA was enriched from 5 µg total RNA using oligo (dT) magnetic beads and fragmented in a thermomixer. The short fragments were used as templates to synthesize first-and second-strand cDNAs. Afterward, cDNA was end-repaired, A-tailed and ligated with Illumina-specific adaptors. Finally, the fragments were size selected and PCR amplified to generate the cDNA library. RNA-sequencing was performed using the Illumina HiSeq™ 2000 platform. Each sample was analyzed in two biological replicates in this RNA-Seq experiment.

De Novo Transcriptome Assembly and Functional Annotation
The raw reads were filtered by removing adaptors and low-quality reads. After that, the clean reads were generated. Trinity (version 2.0.6) (http://trinityrnaseq.github.io/) and TGICL software (version r2013-04-11) (https://sourceforge.net/projects/tgicl/) were employed to assemble our clean data. The clean reads from the SL and YL libraries were first processed independently. To obtain complete reference sequences, the clean reads from all the samples were mixed and assembled again.
To assess the quality of de novo assembly, the length distribution of assembled contigs and unigenes were collected. To obtain functional annotation of a given unigene, the sequence was aligned against protein sequence database entries including those in the Nr, Swiss-Prot, KEGGand COG databases using BLASTx (E-value < 0.00001). The unigenes were also aligned against the Nt database using BLASTn and we declared sequences similar if the corresponding E-value for the alignment was less than 10 −5 .

Screening of Differentially Expressed Unigenes
Gene expression was calculated using the number of reads aligned to a single gene and the total number of reads aligned to reference sequences in the reads per kb per million reads (RPKM) method [50]. R package DEGseq was employed to identify DEGs with random sampling model [51]. A p-value can denote its expression difference between the two libraries, while false discovery rates (FDRs) were used to determine the threshold of p value. We set "FDR < 0.001 and the absolute value of log2Ratio > 1" as the threshold to judge the significance of gene expression difference according to described method in the literature [52].

Validation of RNA-Seq Data by Quantitative Real-Time PCR (qRT-PCR)
To validate the accuracy of the gene expression levels of DEGs obtained from the RNA-Seq analysis, 19 genes possibly associated with luteoloside synthesis were randomly selected and subjected to qPCR detection. Gene-specific primers for the selected genes were designed using an online primer design software (https://www.genscript.com/ssl-bin/app/primer) (Supplementary Table S5) and a melting curve analysis was used to confirm specificity. qRT-PCR was performed using a Fast SYBR Mixture (CWBIO, Beijing, China) on a Bio-Rad CFX connect real-time PCR detection system using of 95 • C incubation for 10 min, then 40 cycles of 95 • C for 15 s and 60 • C for 60 s. For all qPCR experiments, three biological replicates were performed. Relative expression levels were calculated based on the 2 −∆∆Ct method using tubulin as a reference gene.

Conclusions
In this study, we demonstrated senescing leaves in L. macranthoides can accumulate very large amounts of luteoloside, which was extremely higher than those in other leaf samples and organs including stems and flowers. Transcriptome analysis of senescing leaves and young leaves screened a subset of candidate genes associated with luteoloside biosynthesis. The elevated mRNA levels of twenty-four unigenes including PAL2, C4H2, thirteen 4CLs, CHS2, six F3 Hs and two UFGTs were considered to contribute to luteoloside accumulation in senescing leaves with respect to young leaves. Thus, we hypothesized that luteoloside may be preferentially biosynthesized using the following route: naringenin-eriodictyol-luteolin. Furthermore, several unigenes encoding MYB, bHLH and WD40 TFs, such as MYB12, MYB75, bHLH113 and TTG1, involved in flavonoid biosynthesis, were coexpressed with luteoloside biosynthetic unigenes in our dataset. Their roles in regulating luteoloside biosynthesis in L. macranthoides will be characterized via genetic and biochemical approaches in vitro and in vivo in our future work.
Supplementary Materials: Supplementary materials can be found at http://www.mdpi.com/1422-0067/19/ 4/1012/s1. Table S1. Overview of RNA-Seq data; Table S2. Statistical results for functional annotation of transcriptome in Lonicera macranthoides; Table S3. Differentially expressed genes (DEGs) identified by transcriptomic approach during leaf senescence in Lonicera macranthoides; Table S4. DEGs involved in WD40 transcription factors during leaf senescence in Lonicera macranthoides; Table S5. Primer sequences used for qRT-PCR. Figure S1. HPLC chromatogram of luteoloside with detection at 360 nm; Figure S2. Length distribution of All-Unigene sequence; Figure S3. Statistical analysis results of unigene sequences against the NCBI non-redundant protein (Nr) database; Figure S4. The Pearson and Spearman correlation coefficients of gene expression data between the biological replicates. Figure S5. Identification of significantly altered KEGG pathways during leaf senescence in Lonicera macranthoides. The left Y-axis indicates the top 20 enrichment KEGG pathways between young and senescing leaves. The different dot sizes indicate the number of DEGs corresponding to each KEGG pathway; Figure S6. Phylogenetic analysis of putative F3 Hs and UFGTs in Lonicera macranthoides. The phylogenetic tree was constructed using MEGA software (version 5.1) based on Neighbor-Joining method. Values above the branches are bootstrap percentages (1000 replicates Acknowledgments: This work was supported by the National Natural Science Foundation of China (31770338) and The Chongqing Natural Science Foundation (cstc2017jcyjBX0026).
Author Contributions: Zexiong Chen, Ning Tang and Zhengguo Li designed this work. Zexiong Chen, Ning Tang and Guohua Liu performed the experiment. Zexiong Chen prepared the original draft the manuscript. Ning Tang and Zhengguo Li revised the manuscript. All authors read and approved the final manuscript.

Conflicts of Interest:
The authors declare no conflict of interest. The funding sponsors had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; and in the decision to publish the results.