Transcriptional regulation of Lonicera japonica Thunb. during flower development as revealed by comprehensive analysis of transcription factors

Abstract

Background

Lonicera japonica Thunb. flower has been used for the treatment of various diseases for a long time and attracted many studies on its potential effects. Transcription factors (TFs) regulate extensive biological processes during plant development. As the restricted reports of L. japonica on TFs, our work was carried out to better understand the TFs’ regulatory roles under different developmental stages in L. japonica.

Results

In this study, 1316 TFs belonging to 52 families were identified from the transcriptomic data, and corresponding expression profiles during the L. japonica flower development were comprehensively analyzed. 917 (69.68%) TFs were differentially expressed. TFs in bHLH, ERF, MYB, bZIP, and NAC families exhibited obviously altered expression during flower growth. Based on the analysis of differentially expressed TFs (DETFs), TFs in MYB, WRKY, NAC and LSD families that involved in phenylpropanoids biosynthesis, senescence processes and antioxidant activity were detected. The expression of MYB114 exhibited a positive correlation with the contents of luteoloside; Positive correlation was observed among the expression of MYC12, chalcone synthase (CHS) and flavonol synthase (FLS), while negative correlation was observed between the expression of MYB44 and the synthases; The expression of LSD1 was highly correlated with the expression of SOD and the total antioxidant capacity, while the expression of LOL1 and LOL2 exhibited a negative correlation with them; Many TFs in NAC and WRKY families may be potentially involved in the senescence process regulated by hormones and reactive oxygen species (ROS). The expression of NAC19, NAC29, and NAC53 exhibited a positive correlation with the contents of ABA and H₂O₂, while the expression of WRKY53, WRKY54, and WRKY70 exhibited a negative correlation with the contents of JA, SA and ABA.

Conclusions

Our study provided a comprehensive characterization of the expression profiles of TFs during the developmental stages of L. japonica. In addition, we detected the key TFs that may play significant roles in controlling active components biosynthesis, antioxidant activity and flower senescence in L. japonica, thereby providing valuable insights into the molecular networks underlying L. japonica flower development.

Background

Lonicera japonica Thunb. (Caprifoliaceae), a medicinal plant, has long been used in traditional Chinese medicine for the treatment of various diseases, such as influenza, cold, fever, and infections [1]. Modern pharmacological studies have proved that the extracts of L. japonica exhibit therapeutic potency for many biological and pharmacological activities, such as anti-inflammatory [2], antiviral [3], antibacterial, antioxidant [4], hepatoprotective and anti-tumor [5, 6]. Additionally, recent studies indicated that polysaccharides from L. japonica flower buds show hypoglycemic and hypolipidemic effects on streptozotocin-induced diabetic rats, and neuroprotective effects on cerebral ischemia-reperfusion injuries in rats [7, 8]. These findings demonstrated that, L. japonica has received much interest in recent years in specialized pharmacological research studies.

L. japonica is a perennial, evergreen and twining vine which has double-tongued flowers that open white and fade to yellow [9]. The main components of L. japonica include essential oils, phenolic acids, flavone, triterpenoid saponins, iridoids, and inorganic elements, which are considered to be closely related to its pharmacological effects [5]. In L. japonica, the accumulation of some active components is various during floral development. For example, the contents of chlorogenic acid (CGA) and luteoloside, the standard chemicals for evaluating the chemical quality of L. japonica [10], are highest in slightly white alabastrum, while lower in other developmental stages [11]. Although studies have been conducted to analyze the biosynthesis of active compounds by transcriptome [12,13,14], the regulatory mechanism of secondary metabolic pathways and the physiological processes during different developmental stages of L. japonica remain largely unknown.

The expression levels of genes involved in secondary metabolism are often regulated through developmental, environmental or hormonal processes [15]. As sequence-specific DNA-binding proteins, transcription factors (TFs) function by interacting with regulatory regions and modulating the rate of transcriptional initiation in target genes [16]. Publications showed that several families of TFs play important roles in controlling the biosynthesis and accumulation of secondary metabolites [17]. In Arabidopsis MYB family, AtMYB12, AtMYB114, and AtMYB90 are involved in regulating anthocyanin biosynthesis via activation of the entire phenylpropanoid pathway [18, 19]. MYB and bHLH TFs combine with WD40 proteins to form MYB–bHLH–WDR protein complexes (MBW), which regulate flavonoid metabolism processes [20]. However, transcriptomic research that focused on TFs in L. japonica is still very limited.

Furthermore, TFs regulate physiological and developmental phenotypes by interacting with plant hormones [21]. Flower development is determined by the combined action of multiple pathways, which involve floral homeotic genes and hormone signaling molecules [22]. TFs and phytohormones including auxins, gibberellin (GA), abscisic acid (ABA), salicylic acid (SA) and jasmonic acid (JA), all play important roles in flower development [23]. For examples, the expression of TF NAC19 is strongly induced by ABA and slightly induced by JA, and NAC19 is involved in mediating programmed cell death (PCD) by reactive oxygen species (ROS) accumulation in plant cells [24]; WRKY53 is modulated by the JA and SA equilibrium in a complex TF signaling network that regulates plant senescence [25]. Yet relevant researches for exploring the physiological and developmental regulation during L.japonica growth remain scarce.

In this study, sample materials were collected from L. japonica flowers at representative five different developmental stages: the young alabastrum stage (S1), the green alabastrum stage (S2), the whole white alabastrum stage (S3), the silvery flower stage (S4), and the golden flower stage (S5) [26]. In order to investigate the regulation roles of TFs on secondary metabolism and development in L. japonica flower, the contents of CGA and luteoloside were measured, and TFs were identified and functionally categorized from RNA-seq data with bioinformatics techniques. What’s more, DETFs (differentially expressed TFs) were clustered based on their expression profiles, and functional annotation in each sub-cluster was performed. Correlation analysis and biological experiments were carried out to confirm the results of TFs analysis.

Materials and methods

Plant materials

Seeds of Lonicera japonica were purchased from Miaopu Seeds Limited Company (Hangzhou, China). More than 100 of L. japonica plants were grown in the Pingyi cultivation base (Shandong Province, China) without exposure to extreme drought, plant diseases, and insect pests. After being authenticated by Professor Lin Zhang at Zhejiang Sci-Tech University, China, flowers of L. japonica from five different developmental stages (the juvenile bud stage (S1), the third green stage (S2), the complete white stage (S3), the silver flowering stage (S4), and the gold flowering stage (S5)) were collected in late May. During sample collection, flowers from at least 10 plants were mixed and regarded as one biological replicate representing each stage, and three independent replicates were performed. Flower materials were frozen in liquid nitrogen immediately after collection and then stored at − 80 °C.

Determination of CGA and luteoloside contents by HPLC

L. japonica flower samples were ground in liquid nitrogen and lyophilized. Then, an accurately weighed powder sample (0.010 g) was suspended in 1 mL MeOH, ultrasonically extracted for 1 h, centrifuged at 13,000 g for 5 min and transferred the supernatant. Repeat the extraction step and combine the two supernatants. The extracting solution was concentrated and re-solubilized in 1 mL MeOH, which resulted in samples for CGA and luteoloside contents measuring.

To determine the concentration, the following procedures were performed on Waters 2695 Alliance HPLC system (USA), which equipped with a photodiode array detector, an online degasser and an autosampler for solvent delivery. An aliquot of sample (10 μL) was injected into Waters symmetry C18 column (250 mm × 4.6 mm I.D 5 μm) with a flow rate of 1 mL/min at 20 °C. The mobile phase was composed of 0.05% phosphoric acid in water (A) and 100% acetonitrile (B). The gradients were as follows: 0 min, 10% B; 50 min, 20% B; 60 min, 35% B; 70 min, 95% B. Spectra were measured at wavelength of 327 nm and 350 nm. Peak area was used for quantification.

Preparation of the standard calibration curves: Chlorogenic acid and luteoloside were purchased from Yuanye Biotechnology Co. (Shanghai, China). Calibration curves were constructed using the standard solutions diluted in MeOH at five different concentrations: CGA (0.06, 0.12, 0.18, 0.3, and 0.42 mg/mL), luteoloside (1.5, 3, 4.5, 9, and 18 μg/mL). The regression equations for CGA and luteoloside were y = 309,521x - 203,748 (r² = 0.9991) and y = 1697.2x - 19,711 (r² = 0.9989), respectively.

RNA isolation and library construction

RNA isolation and library construction were performed according to the method described by Yang et al. (2017) [27]. In brief, a portion (100 mg) of each sample was ground to powder in liquid nitrogen. Total RNA was isolated using RNeasy Plant Mini kit (Qiagen, Hilden, Germany). The quantity and quality of RNA were determined using Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA). The poly (A) mRNA was isolated from total RNA samples with Magnetic Oligo (dT) Beads (Illumina, San Diego, CA, USA) and used for mRNA-sequencing library construction. cDNA was synthesized using the fragmented mRNA with the incorporation of reverse transcriptase for end-repair, followed by a single ‘A’ base addition. mRNA-Sequencing Sample Preparation Kit (Illumina) was used to prepare the DNA fragments for ligation to the adapters. After purification, the cDNA fragments (200 ± 25 bp) were excised and retrieved. Polymerase chain reaction (PCR) was performed to selectively enrich and amplify the cDNA fragments. Libraries were prepared from a 150–200 bp size-selected fraction following adapter ligation and agarose gel separation.

Sequencing, de novo assembly and annotation

The mRNA-sequencing libraries were sequenced using the Hiseq 2000 sequencing platform (Illumina). Raw reads were filtered using BWA [28], with a quality threshold of 30. The Trinity program was used to assemble the clean reads and obtain non-redundant unigenes [29]. Annotation of the assembled unigenes was performed by searching against the NCBI non-redundant protein (NR), Swiss-Prot, TrEMBL, and Pfam databases. The search was conducted using BLAST with an E-value cut-off of 1e-5. Functional annotation by gene ontology terms (GO, http://www.geneontology.org) was carried out using the Blast2GO software [30].

Identification and expression analysis of TFs

The amino sequences of TFs in Arabidopsis thaliana and Coffea canephora were downloaded from the PlantTFDB database [31]. For TF identification in L. japonica, Blast+ was used for sequence alignment with the identity cut-off threshold of 30% and E-value of 1e-5 [32]. Candidates that contain DNA binding domains were recognized by GO annotation for the final TF identification. Differentially expressed TFs (DETFs) between samples were identified using the value of Fragments Per Kilobase of transcript Per Million fragments mapped (FPKM) with |log2(fold change)| > 1, p value ≤0.05 and q value ≤0.05 [33]. The Multiexperiment Viewer software (v4.9) was employed to exhibit the expression profiles of DETFs by clustering [34]. DETFs were clustered by K-means method. Using the total TFs in L. japonica as reference, GO enrichment analysis of DETFs was performed using Agrigo (http://bioinfo.cau.edu.cn/agriGO/) with hypergeometric test and FDR cut-off of 0.05 [35].

Measurement of endogenous hormones in L. japonica

Samples were prepared using the method of Pan et al. with minor modification [36]. Firstly, 50 mg sample was ground into fine power in liquid nitrogen. Then, 500 μL of 1-propanol/H2O/concentrated HCl (2: 1: 0.002, v: v: v) was added, and the mixture was agitated for 30 min at 4 °C. Next, the solution was centrifuged at 13,000 g for 5 min at 4 °C. 1 mL of CH₂Cl₂ was added to the isolated upper layer, and the new mixture was agitated for 30 min at 4 °C then centrifuged at 13,000 g for 5 min at 4 °C. The lower layer (~ 1 mL) was concentrated and re-solubilized in 0.2 mL of 50% MeOH, which resulted in the samples for hormones levels determination.

The following determination was performed on Agilent 6460 Triple Quadrupole LC-MS/MS (Agilent Technologies, Technologies, Palo Alto, CA, USA) with an electrospray ESI (Agilent Technologies). An aliquot (20 μL) was injected into Zobax XDB C18 column (2.1 mm × 150 mm × 3.5 μm, Agilent Technologies) with a flow rate of 0.3 mL/min. Peak abundance was used for ABA, JA, and SA quantification. The mobile phase was 0.1% formic acid in water (A) and 100% MeOH (B). The gradient was as follows: 0.00 min, 60% A (40% B); 1.5 min, 60% A (40% B); 8 min, 0% A (100% B); 8.5 min, 60% A (40% B).

Measurement of hydrogen peroxide (H₂O₂) and total antioxidant capacity in L. japonica

To measure the H₂O₂ content and the total antioxidant capacity, 100 mg sample was grounded with 0.9 mL deionized distilled H₂O (dd H₂O) on ice. The resulting mixture was centrifuged at 10,000 g for 10 min at 4 °C, and the upper layer was isolated for the subsequent process. Protein concentration was determined using the Bradford assay with bovine serum albumin as standard [37]. The Hydrogen Peroxide Assay Kit (Nanjing Jiancheng Bioengineering Institute, Nanjing, China) and the Total Antioxidant Capacity (T-AOC) Assay Kit (Nanjing Jiancheng Bioengineering Institute, Nanjing, China) were used for measuring the H₂O₂ content and the total antioxidant capacity.

Statistical analysis

Significant differences between contents of CGA, luteoloside, hormones and H₂O₂, and total antioxidant capacity were calculated using a one-way ANOVA analysis with a Turkey test and a significance level at α = 0.05 in SPSS software. All expression analyses were performed in three replicates. Reported values represent arithmetic averages of three replicates. Data was expressed as mean plus or minus standard deviation (mean ± SD).

Results

Contents of CGA and luteoloside in different developmental stages of L.japonica flower

The contents of CGA and luteoloside in S1, S2, S3, S4 and S5 of L.japonica flower were measured by HPLC (Fig. 1). The morphological photos corresponding to each stage were shown in Additional file 1: Figure S1. According to the results, the content of CGA was higher at S2 and S3, while was lowest at S4. The content of luteoloside at S3 was the highest, and about three times higher than other stages. Overall, the content of CGA exhibited an increase from S1 to S3, while reached its minimum at S4, and then increased again at S5; The content of luteoloside increased from S1 to S3, followed by a decrease from S3 to S5.

Fig. 1

Contents of CGA and luteoloside in five different developmental stages of the L. japonica flower. Data is shown as the mean ± SD from three independent experiment replicates. a, b, c indicate significant changes measured by one-way ANOVA analysis (α < 0.05)

Sequencing, assembly and annotation of L.japonica transcriptome

To acquire comprehensive TF expression profiles of L.japonica flowers at five stages of development, transcriptome sequencing and analysis were performed. Clean pair-end reads were obtained after removal of adaptor sequences and low-quality reads. The Trinity program was used for de novo assembly of all clean reads, which yielded a total of 43,689 unigenes with an average length of 1033 bp, N50 of 1612 bp, and GC content of 41.49% (Table 1). The assembled unigenes were annotated by common databases including the Swiss-Prot, TrEBML, NR, Pfam and GO, to which approximately 39.00, 65.30, 65.80, 53.20, and 54.60% of unigenes were mapped respectively (Table 2). Among the total unigenes, 15,585 (36%) were annotated by all the five databases. The number of unigenes annotated by only one database were 0, 19, 227, 0, 0 for Swiss-Prot, TrEBML, NR, Pfam and GO, individually (Additional file 1: Figure S2).

Table 1 Assembly summary of transcriptomic data

Table 2 Annotation results of the unigenes

Identification and classification of TFs

Many TFs have been found to play important roles in plant growth [38]. To detect the biological functions of TFs involved in the development of L. japonica flower, TFs expressing in five different developmental stages were identified. To ensure the accuracy of TF identification and classification, we aligned the unigenes of L. japonica with the amino acid sequences of TFs in two species including A. thaliana and C.canephora. A. thaliana is a model plant. C. canephora is one of the top species that has the most homologous genes with L.japonica in NR annotation, indicating that C.canephora may share a close polygenetic relationship with L. japonica (Additional file 1: Figure S3). Also, according to Flora of China, C. canephora and L.japonica both belong to Rubiales Family, which strengthened the close relationship between them [39]. The results showed that, 1391 unigenes were mapped to the amino acid sequences of A. thaliana’s TFs, and 1407 unigenes were mapped to the amino acid sequences of C. canephora’s TFs. As shown in Fig. 2a, we selected the 1316 unigenes of L. japonica which were mapped to both A. thaliana and C.canephora for further analysis.

Fig. 2

Identification and distribution of TFs. a Assembled L. japonica unigenes were aligned with the amino sequences of TFs in A. thaliana and C. canephora by Blast+ with an E-value cut-off of 1e-5. b Distribution of TFs in five developmental stages of L. japonica flower. The TFs identified were based on the intersection results of A. thaliana and C. canephora shown in (a), a total of 1316 TFs were detected. The five stages were the juvenile bud stage (S1), the third green stage (S2), the complete white stage (S3), the silver flowering stage (S4), and the gold flowering stage (S5)

Among the 1316 TFs, 1279, 1275, 1268, 1241, 1254 TFs were found in S1, S2, S3, S4 and S5, respectively. 1185 TFs were expressed in all five stages, representing 90.04% of the total TFs. S5 contained the most stage-specific TFs (12), which were from TF families including YABBY, GRAS, MIKC_MADS, ERF, HD-ZIP and SBP (Fig. 2b, Table 3). All the analyzed TFs (1316) were classified into 52 families. The five largest families represented were bHLH (108), ERF (95), MYB (89), bZIP (75), MYB_related (66) (Fig. 3a).

Table 3 S5 specific TFs

Fig. 3

Classification and differential expression of TFs in five stages. a Distribution of total TFs and differentially expressed TFs (DETFs) among different TF families. X-axis represents number of TFs. b Pairwise comparison of TFs in five developmental stages. Y-axis represents number of DETFs

Analysis of differentially expressed TFs (DETFs)

To acquire insights into functional and regulatory dynamics during flower development, pairwise differential analysis was conducted on the expression levels of TFs in the five developmental stages. Results showed that 917 TFs, classified into 48 families, were differentially expressed in at least one of the pairwise comparisons among the five flowering stages, representing 69.68% of the total TFs. Most DETFs (86) belonged to the bHLH family, followed by ERF (71), MYB (72), bZIP (49), NAC (46). In addition, expression levels of TFs from some families, such as TCP, CO-like, LSD, EIL were found to fluctuate significantly among the five stages. Lower portions of the TFs in MYB_related, C3H, SBP, FAR1, BBR-BPC and E2F/DP families were differentially expressed. Meanwhile, expression levels of TFs in HB-PHD, GeBP, STAT and VOZ families remained relatively stable throughout development (Fig. 3a).

To further decipher the expression profiles of TFs, the up- or down-regulated TFs in each pairwise comparison were identified (Fig. 3b). The comparison S1_VS_S4, S2_VS_S4, S1_VS_S5, S2_VS_S5 exhibited more DETFs, while S1_VS_S2, S4_VS_S5 showed fewer DETFs. This result indicated that S4 and S5 shared relatively similar situation during L. japonica flower development, so were S1 and S2. This observation was consistent with the cluster result shown in Additional file 1: Figure S4. Overall, an increasing number of up-regulated TFs was found along with the growth of L. japonica. For example, when expression of TFs in all stages were compared to that in S1, 72 DETFs in S1_VS_S2, 211 DETFs in S1_VS_S3, 245 DETFs in S1_VS_S4, and 277 DETFs in S1_VS_S5. This finding implied TFs might be involved in more activities in flowering development of L. japonica.

Time-course expression of DETFs

To clarify the time-course expression of DETFs in L. japonica, the 917 DETFs were grouped into six clusters based on their temporal expression profiles (Fig. 4). Cluster 1 contained 145 TFs that achieved their maximum expression at S2, followed by a decrease with further development. Cluster 2 contained 119 TFs that exhibited an increase in expression level from S1 to S3, followed by a decrease from S3 to S5. Cluster 3 contained 194 TFs that peaked at S1 and then steadily decreased across the five-time points. Cluster 4 was comprised of 116 TFs whose expression decreased sharply from S1 to S2, slowly at S3 and S4, and then increased slightly from S4 to S5. Cluster 5 contained 159 TFs that increased in expression from S1 to S4, and then exhibited a decrease in S5. The expression of TFs belonging to cluster 6 exhibited little changes between S1 and S2, followed by a surge from S3 to S5. Cluster 1 to 6 contained 13.0, 21.2, 15.8, 12.6, 17.3, and 20.1% of all DETFs, respectively. The largest TF family represented in cluster 1, 3 and 5 was bHLH and that in cluster 4 and 6 was NAC. Furthermore, several large TF families such as NAC, ERF, WRKY, MYB, bHLH were found to made up more than half of the TFs in cluster 4. In cluster 6, NAC, ERF, bHLH were the top three most abundant families (Additional file 1: Figure S5).

Fig. 4

Clusters of differential expression TFs (DETFs). DETFs were clustered and displayed in line chart using Mev. Six clusters were grouped. X-axis represents the five different developmental stages of L. japonica flower including S1, S2, S3, S4 and S5. Y-axis represents normalized value of the DETFs’ expression level. The middle black line in each cluster is the zero line and the red line indicates the average expression level. The number in the upper left corner represents the number of TFs in each cluster

GO enrichment analysis of time-course expression patterns

To further elucidate the functional mechanism of TFs in mediating flower development, GO enrichment analysis was performed on the six clusters of DETFs. Our results showed that, in cluster 1, there were no significantly enriched GO items; Many TFs in cluster 2 were associated with “protein self-association”; In cluster 3, the significant functions of TFs were “lipid binding” and “cotyledon development”; TFs in cluster 5 might play important roles in “ligase activity” (Additional file 2: Table S1). TFs in cluster 4 and cluster 6 were involved in multiple functions. Specifically, in cluster 4, the results of GO enrichment analysis indicated that TFs in this group were primarily associated with various stimuli, responses and biosynthetic processes, including JA stimulus, SA stimulus, ABA stimulus, carbohydrate stimulus, biotic stimulus, response to ROS, and organic acid biosynthetic process (Fig. 5). In cluster 6, most of the TFs were involved in the regulation of PCD, JA-mediated signaling pathway and ABA-mediated signaling pathway (Fig. 6). These findings demonstrated that the TFs in cluster 4 and cluster 6 might have essential roles in plant hormones, ROS and PCD during the growth and development of the L. japonica flower.

Fig. 5

Go enrichment of TFs in cluster 4. GO enrichment analysis was performed using Agrigo with FDR value cut-off of 0.05 and the total TFs in L. japonica flower were used as reference. Color of circle in the plot indicates FDR value and size of circle indicates the number of TFs in each GO_Term

Fig. 6

GO enrichment of TFs in cluster 6. GO enrichment analysis was performed using Agrigo with FDR value cut-off of 0.05 and the total TFs in L. japonica flower were used as reference. Color of circle in the plot indicates FDR value and size of circle indicates the number of TFs in each GO_Term

Endogenous hormones measurements

TFs and plant hormones have been shown to influence each other during plant growth in a complex relationship [21]. In order to better understand the interaction between endogenous hormones and TFs in the transcription activities during L. japonica flower development, the contents of hormones including JA, SA, and ABA were measured (Fig. 7). The results showed that both JA and SA were found at low concentrations in the early stages of flower development, and then a surge at S3 was observed, followed by a decrease from S4 to S5. The difference between the contents of JA and SA was that, in the last two periods, JA content decreased slightly while that of SA decreased dramatically. The content of ABA reached a plateau at S1 and S2, then steadily increased until a maximum was reached at S5. Moreover, the trend of ABA concentration throughout development was very similar to the expression trend of TFs observed in cluster 6 (Figs. 4 and 7).

Fig. 7

Contents of endogenous hormones including JA, SA, and ABA in different developmental stages of the L. japonica flower. Data is shown as mean ± SD from three independent experiment replicates. a, b, c indicate significant changes measured by one-way ANOVA analysis (α < 0.05)

Determination of hydrogen peroxide (H₂O₂) content and the total antioxidant capacity

In plants, various abiotic stresses can lead to the formation of ROS, which has been proposed as key inducers of developmental PCD [40]. H₂O₂ is one of the main members of ROS, and its effect on PCD has been well studied. Antioxidant system is important in the maintenance of balanced ROS levels [41]. In our research, the TFs in cluster 4 and 6 were found to have a close association with H₂O₂ and PCD (Fig. 6 and 7). To investigate the association among them, the content of H₂O₂ and the total antioxidant capacity were evaluated (Fig. 8). Our results showed both of them followed a similar trend, with low levels observed at the early stages of flower development, followed by high levels at S4 and S5. In addition, these trends were consistent with the trend of ABA’s content and the expression pattern of TFs in cluster 6 (Figs. 4 and 8).

Fig. 8

Content of hydrogen peroxide (H₂O₂) and the total antioxidant capacity in different developmental stages of the L. japonica flower. Data is shown as mean ± SD from three independent experiment replicates. a, b, c, d, e indicate significant changes measured by one-way ANOVA analysis (α < 0.05)

Correlation analysis of specific TFs

Based on the cluster and enrichment results, the key TFs involved in secondary metabolism and L. japonica flower development were detected. Figure 9 showed the expression patterns and their correlations. CHS and FLS have been reported as the key enzymes in the flavonol pathway [42]. Our results demonstrated that MYB114, MYB12 and MYB44 might be involved in flavonoid biosynthetic processes (Additional file 3: Table S2). In this study, the expressions of CHS and FLS exhibited similar time-course patterns to the expression of MYB12 from S1 to S5, while a reverse trend was observed when compared with the expression of MYB44 at S1, S2 and S3 (Fig. 9a). Additionally, NAC19, NAC29 and NAC53 might be involved in JA- and ABA-mediated signaling pathways and PCD (Table 4). The expressions of the three TFs in cluster 6 and the contents of ABA and H₂O₂ were very consistent, and JA content was also associated with the expression of these TFs (Fig. 9b). Moreover, the expression of LSD1 and SOD was highly correlated, while the expression of LOL1 and LOL2 displayed a negative correlation with them (Fig. 9c). In addition, the TF WRKY70 might be involved in ABA mediated signaling pathway, JA stimulus, and SA stimulus. WRKY53 might be involved in JA stimulus, and WRKY54 might be involved in JA stimulus, SA stimulus, and ROS (Table 4). These three TFs belonged to cluster 4 and showed similar expression pattern that exhibited high expression at S1 and followed by lower expression in later stages (Table 4).

Fig. 9

Correlations among TFs, genes, hormones, and H₂O₂. The expression of the genes, the contents of hormones and ROS were normalized by “Z-score”. a The expression patterns of MYB12, MYB44, CHS and FLS, which are involved in flavonoid biosynthesis. b The correlations among the expression of NAC19, NAC29 and NAC53, the contents of JA, ABA, and H₂O₂, which are involved in senescence processes. c The correlations among the expression of LSD1, LOL1, LOL2 and SOD. The number in each square represents the corresponding correlation coefficient. Color intensity is proportional to the correlation coefficient. Positive correlations are displayed in red and negative correlations in blue

Table 4 Enrichment results of TFs associated with hormones and PCD in cluster 4 and cluster 6

Discussion

Comprehensive analysis of TF expression profiles throughout L. japonica flower development

TFs are proteins that bind specific DNA sequences and regulate transcription. Recent literature demonstrated that TFs play significant roles in the growth and development of plants [43]. In A.thaliana genomes, nearly 6% of genes encode TFs [44]. Nevertheless, our knowledge on TFs in non-model plants are limited. In this study, TFs in five different developmental stages of L. japonica flower were investigated in order to better understand the roles of TFs L. japonica flower. To ensure the accuracy of identification and classification of TFs, the unigenes of L. japonica flower were aligned with the amino acid sequences of two other species including A.thaliana and C.canephora. Our results showed that 1316 TFs were identified and classified into 52 families (Fig. 2a, Fig. 3a). According to PlantTFDB, in A.thaliana, there were 2296 TFs (1717 loci) that belonging to 58 families. In C.canephora, which has a closer relationship to L. japonica, 1256 TFs (1256 loci) were identified and classified into 57 families [31].

Each stage of the L. japonica flower is characterized by different properties, which may be related to the regulatory activities of TFs. In our study, the accumulations of CGA and luteoloside, the most important active components in L. japonica, were variable at different stages (Fig. 1). Expression profiles of TFs were analyzed at the five different stages: S1, S2, S3, S4 and S5. Among the 1316 total TFs, 1279, 1275, 1268, 1241 and 1254 TFs were expressed in S1, S2, S3, S4 and S5, respectively (Fig. 2b). Additionally, 12 stage-specific TFs were found at S5. These 12 TFs included 3 ERF (ERF3), 3 SBP (SPL1), 2 YABBY (YAB2), 2 MIKC_MADS (AGL8), 1 GRAS (SCL15), 1 HD-ZIP(HAT7) (Table 3). SPL1 plays an important role in inflorescence [45], YAB2 is expressed in floral organ primordium and SCL15 represses embryonic traits in seeding acting [46, 47]. All of these TFs are related to the developmental processes in flowers, indicating that they may play important regulatory roles in S5.

Among the 52 families in the total 1316 TFs, the five largest families were bHLH (108), ERF (95), MYB (89), bZIP (75), and MYB-related (66) (Fig. 3a). All of these families are large TF families in plants and are involved in the regulation of numerous biological processes [48, 49]. As for the 917 DETFs, the five largest families were bHLH (86), ERF (71), MYB (72), bZIP (49), and NAC (46) (Fig. 3a). Compared with the five largest families of the total TFs, bHLH, ERF, MYB, bZIP were the constant, while NAC replaced MYB-related as the fifth largest family of the DETFs. The NAC family is the second-largest TF family in plant, and the TFs in NAC regulate diverse processes, such as floral development and auxin signaling [50]. This implies that the TFs in NAC family are active and may have significant roles in regulatory activities during the development of L. japonica flower.

The transcriptional regulation of active compounds biosynthesis by TFs in L. japonica flower

CGA and luteoloside are the standard chemicals for evaluating the chemical quality of L. japonica [10]. In this study, the content of CGA was higher at S2 and S3, and the content of luteoloside was highest at S3 (Fig. 1). CGA is a major member of soluble plant phenolics [51], and luteolin is an important plant flavonoid [52]. Both biosynthesis of CGA and luteolin derive from the phenylpropanoid pathway. Phenylpropanoids comprise an important class of plant secondary metabolites. A number of TFs have been reported to regulate branches of phenylpropanoid metabolism [53]. MYB–bHLH–WDR protein complexes (MBW) that consist of TFs of MYB and bHLH families and WD40 proteins regulate phenylpropanoids biosynthesis [54]. This process is highly conserved in higher plants [20]. In this study, most of the TFs involved in flavonoid pathway belonged to bHLH or MYB families, indicating that the important roles of these complexes in L. japonica flower (Additional file 3: Table S2, Additional file 1: Figure S5).

TFs MYB114 and WRKY44 are members of MBW. Previous study demonstrated that MYB114 increases accumulation of anthocyanin by activating the expression of anthocyanin biosynthesis genes [20]. WRKY44, also known as TTG2, controls proanthocyanidin biosynthesis [55]. In this study, MYB114 belonged to c1 cluster, the expression of it peaked at S3 and exhibited a high correlation with the content of luteoloside, while WRKY44 belonged to c2 cluster, the expression of it decreased from S1 to S5 (Figs. 1 and 4). These results indicate that MYB114 may play an important role in luteoloside biosynthesis and the regulatory mechanisms of MBW in L. japonica flower are complex.

Moreover, MYB12, also known as PFG1, is a member of MBW, too [56]. It has been reported that MYB12 is involved in phenylpropanoid pathway in Arabidopsis, tobacco, and tomato plants, acting as a very effective and positive regulator in the biosynthesis of caffeoyl quinic acids and flavonols [19]. CHS and FLS are key enzymes in the phenylpropanoid pathway. MYB12 specifically activates the flavonol pathway by inducing the expression of CHS and FLS [57]. MYB12 and CHS shows same transcription profiles in Asiatic hybrid lily during its flower bud development [58]. In tomato plants, AtMYB12 leads to the induction of the biosynthetic genes required for the production of flavonols and CGA [19]. In this study, the expression of MYB12, CHS and FLS was found to be highly consistent with each other (Fig. 9a). This result indicates that, in L. japonica, MYB12 may play a significant role in controlling the expression of CHS and FLS in different developmental processes. On the other hand, flavonoid biosynthesis is also controlled by negative regulation as the repression of JA-mediated defense by MYB44 contributes to the low accumulation [59]. The expression of MYB12, CHS, and FLS was maintained higher levels in S2, while the expression of MYB44 decreased in S2 (Fig. 9a). This result indicates that MYB44 may have a negative role in regulating the gene expression of CHS and FLS in L. japonica.

Besides, antioxidants have emerged as prophylactic and therapeutic agents for various diseases. Previous study demonstrated that extracts from the L. japonica flower exhibit antioxidative activity [6]. However, little is known about the transcriptional regulation concerning the antioxidant properties of L. japonica. LOL1, LOL2 and LSD1 belong to LSD TF family. LSD1 can positively regulate the accumulation of SOD, which is an enzyme capable of scavenging O2˙- and leading to H₂O₂ formation [60, 61]. On the contrary, LOL1 and LOL2 play negative roles in SOD accumulation [61]. In this study, the expression of LSD1 and SOD, the content of H₂O₂ and the total antioxidant capacity were positively correlated, while negatively correlating with the expression of LOL1 and LOL2 (Fig. 9c). Our results were in accordance with those from previous researches, indicating the regulatory roles of LSD1 LOL1, and LOL2 on antioxidation in L. japonica flower.

The transcription regulation of PCD by TFs in NAC and WRKY families in L. japonica flower

The mechanisms of plant developmental PCD have been extensively studied [62]. Flowers have a species-specific, limited life span with an irreversible program of senescence. To some degree, the terms ‘senescence’ and ‘PCD’ are equivalent for flowers [63]. NAC and WRKY TFs have been closely associated with senescence in several tissues such as Arabidopsis leaves, petals and siliques, and they interact with many regulators including plant hormones and ROS [17, 64]. Plant hormones like JA, SA, and ABA act as positive regulators of plant senescence [17, 65]. ROS have been proposed as key inducers of developmental PCD or senescence [40]. In L. japonica flower, many TFs known to involve in PCD or senescence were detected (Table 4). Moreover, our results support the complex relationships between TFs and many regulators like JA, SA, ABA and ROS that are involved in controlling PCD L. japonica flower.

In L. japonica, many TFs in cluster 6 were associated with PCD (Fig. 6). These TFs belonged to different families, including bZIP (4), NAC (2), MIKC_MADS (2), LSD (1), C2H2 (1), and TALE (1) (Table 4). These families have been reported to play regulatory roles in plant senescence [66,67,68,69,70]. In addition, the expression of these TFs was highest at S5, the last stage before senescence of L. japonica flower, indicating that these TFs play important roles in the process of senescence. As shown in Table 4, TFs in NAC family that associated with PCD included NAC19 and NAC29. The expression of NAC19 is strongly induced by ABA, slightly induced by JA, and is known to mediate PCD through ROS accumulation [24]. NAC29 can be induced by ABA and may have a positive role in precocious senescence [71]. In addition to NAC19 and NAC29, we also identified NAC53 in cluster 6. NAC53 promotes ROS production in leaves [72], which may further trigger PCD or senescence of plant. In L. japonica, high correlations among the expression of NAC19, NAC29, NAC53, the contents of ABA, JA, and H₂O₂ were observed (Fig. 9b). These results support that NAC19, NAC29, NAC53 may be important in the regulatory network involving PCD and hormones in the senescence process of L. japonica, a potential topic worthy of further analysis.

Furthermore, several WRKY TFs involved in regulation of leaf senescence and response to hormones were also identified in L. japonica flower (Table 4). WRKY70 and WRKY54 are regulated by SA and JA, and co-operate as negative regulators of senescence [73]. In this study, the expression of WRKY70 and WRKY54 were low at S3, S4 and S5, while the contents of SA and JA, the positive regulators of plant senescence, were high at S3, S4 and S5 (Figs. 4 and 7). This implied the negative regulation by WRKY 70 and WRKY 54 on senescence in L. japonica flower. Additionally, WRKY53 is modulated by the JA and SA equilibrium in a complex TF signaling network that regulates senescence [25, 74]. It has been reported that WRKY53, WRKY54, and WRKY70 might be involved in a regulatory network that integrates internal and environmental signals to modulate the onset and the development of leaf senescence [73]. In this study, our results indicate that this regulatory network may also have significant roles in L. japonica flower.

Taken together, the expression of TFs and the presence of different flower development regulators, including ABA, JA, SA and ROS, as well as their interactions, demonstrate complexity of the regulatory network involved in flower development. Further investigations on the various associations between these elements are essential to better understand the potential mechanisms involved in the medicinal effect of L. japonica flower.

Conclusion

TFs play significant roles during plant growth and development. To better understand the TFs’ regulatory roles in L. japonica flower, we aim to comprehensively characterize the expression profiles of TFs at different developmental stages of L. japonica flower. General analysis identified 1316 TFs that were classified into 52 families. The largest families included bHLH (108), ERF (95), MYB (89), bZIP (75), and MYB-related (66). S5, the last developmental stage investigated in this study, possessed the most stage-specific TFs. From the DETFs analysis, TFs from bHLH, ERF, MYB, bZIP, and NAC families exhibited obviously altered expression during the growth of flower. Moreover, cluster and GO enrichment analyses were performed on the DETFs, and 6 clusters were grouped. Many GO terms associated with JA, SA, ABA, ROS and PCD were enriched in cluster 4 and cluster 6. In addition, several significant TFs involved in the biosynthesis of active components, the antioxidant activity, and the development of the L. japonica flower were detected based on the enrichment results. In phenylpropanoids biosynthesis of L. japonica flower, MYB114, WRKY44, MYB12 and MYB44 may have regulatory roles. TFs in LSD family, including LOL1, LOL2, and LSD1, were found to exhibit potentially antagonistic effects on SOD accumulation and antioxidation in L. japonica flower. During the flower senescence processes, TFs in NAC and WRKY families, including NAC19, NAC29, NAC53, WRKY70, WRKY54, and WRKY53, may play significant roles in the regulation of hormones. Taken together, our results serve as valuable resource for further studies.