Transcriptomic insights into the regulatory networks of chilling-induced early flower in tobacco (Nicotiana tabacum L.)

ABSTRACT

Appropriate timing of flowering is pivotal for tobacco, while chilling stress occurring at the seedling stage often undesirably leads to early flowering. However, the potential mechanism underlying chilling-induced early flowering remains unknown. Here, transcriptome sequencing was performed in tobacco with or without chilling at both seedling and budding stages. Chilling affected the expression of numerous genes at the seedling stage, while these dramatic expression changes were largely eliminated at the budding stage. A small number of genes related to metabolism, flower development, and stress tolerance continued to keep their altered expression patterns from the seedling stage to the budding stage. Many potential flowering-related genes involved in flowering pathways were identified and over half of them were differentially expressed. Functional analysis revealed that the down-regulation of NbXTH22 rendered tobacco less sensitive to chilling-induced early flowering. These results provide valuable resources for the investigation of flowering regulatory mechanisms and contribute to the genetic improvement of crops.

KEYWORDS:

• chilling stress
• tobacco
• transcriptome
• RNA-Seq
• early flowering

1. Introduction

Flowering marks the developmental transition from vegetative growth to reproductive growth and is one of the most important events in the plant life cycle. The initiation of floral development in plants is regulated to maximize their reproductive success by integrating endogenous signals with environmental cues, e.g. temperature, light, developmental stage, and phytohormones (Fornara et al. 2010). The genetic basis underlying flowering regulation has been extensively studied in model plants, such as Arabidopsis (Srikanth and Schmid 2011), rice (Hori et al. 2016), maize (Dong et al. 2012), soybean (Jung et al. 2012), leading to the wide identification of flowering-related genes and insightful understanding of floral induction regulatory networks. In Arabidopsis, six flowering pathways are considered to regulate the floral transition process, including vernalization, photoperiod, autonomous, gibberellic acid (GA), thermosensory, and ageing (Fornara et al. 2010). Vernalization, referring to the initiation or acceleration of flowering that occurred after the prolonged exposure of seeds or young seedlings to low temperature, is one of the main environmental cues leading to floral transition in nature (Kim et al. 2009).

In winter-annual Arabidopsis, FLOWERING LOCUS C (FLC) is a key regulator in the vernalization pathway (Sheldon et al. 1999), which inhibits flowering by repressing the expressions of floral promoters such as FLOWERING LOCUS T (FT) and SUPPRESSOR OF OVER-EXPRESSION OF CONSTANS 1 (SOC1). Analogously, in winter wheat and barley, VERNALIZATION2 (VRN2) is a major floral repressor of vernalization-mediated flowering (Yan et al. 2004). As for perennial tree peach (Prunus persica), the DORMANCY ASSOCIATED MADS-box (DAM) genes are considered alternatives to FLC and VRN2 in regulating vernalization-mediated flowering. Intriguingly, these DAM genes are a cluster of six MADS-box transcription factors with high similarities to Arabidopsis AGAMOUS-LIKE 24 (AGL24) and SHORT VEGETATIVE PHASE (SVP) genes (Sasaki et al. 2011). Arabidopsis SVP is induced by both vernalization and photoperiod and interacts with FLC and FLOWERING LOCUS M (FLM) to repress FT and flowering (Gregis et al. 2009; Tao et al. 2012). In contrast, the wheat SVP-group gene TaVRT2 is suppressed by vernalization (Kane et al. 2007). These studies reveal not only noticeable similarities in basal mechanisms and main regulatory genes of flowering across diverse plants, but also distinct differences in numbers and the roles of individual genes in the flowering regulatory pathways between some species.

Early flowering is a destructive problem that limits vegetative growth and reduces the yield and quality of economic products in tobacco. Low temperature (i.e. chilling/vernalization), an unfavorable environmental condition that can easily occur during the tobacco seedling stage, is a primary cause of early flowering. Unfortunately, the molecular mechanisms underlying tobacco’s early flowering induced by chilling have largely been unexplored, and a way to genetically manipulate the flowering time of tobacco independently of temperature control is still unavailable. With increasing concerns on the genetic mechanism and requirement for molecular insights, several flowering-related genes have been identified for tobacco in the past decades. Jang et al. (2002) isolated two tobacco MADS-box genes (NtMADS4 and NtMADS11), which were demonstrated to play a critical role in promoting flowering in tobacco. Other MADS-box genes (NsMADS2 and NsMADS3), FLOWERING PROMOTING FACTOR 1 (NtFPF1), NtFT5, SOC1, the RING FINGER PROTEIN (NtRCP1), and phytochrome genes (NaPHYA and NaPHYB2) were also successively identified and found to be associated with floral induction or altered flowering time of tobacco (Jang et al. 1999; Smykal et al. 2007; Harig et al. 2012; Wang et al. 2015; Fragoso et al. 2017; Beinecke et al. 2018).

Despite the intensive efforts on identifying the roles of flowering-related genes involved in tobacco floral transition, the mechanism underlying global regulatory networks at the transcriptome level is still poorly understood, especially for the chilling mediated changes that might be responsible for early flowering. To this end, a genome-wide transcriptional analysis was performed on tobacco at different development stages using RNA-Seq. We sequenced the cDNA libraries of chilled and non-chilled tobacco samples and compared their gene expression patterns at both the seedling and budding stages. We identified global differentially expressed genes (DEGs) involved in tobacco floral transition, specifically, mined potential flowering pathway genes. We also analyzed the inheritance of altered gene expression of chilled tobacco from the seedling stage to the budding stage. Functional analysis revealed that the down-regulation of NbXTH22, a homolog gene of NtXTH22 in N. benthamiana, can alleviate early flowering caused by chilling stress to some extent. Our findings are useful for further characterizing the molecular regulatory mechanisms underlying floral induction in tobacco.

2. Materials and methods

2.1. Plant materials and growth conditions

Tobacco (N. tabacum L., cv. NC82) was used for chilling treatment and transcriptional analysis. Plants were grown in a growth chamber with a normal growth condition (28°C, 12-h photoperiod) until four-leaf-stage, and then, were transferred to a growth chamber under a chilling condition (12°C, 12-h photoperiod) for 10 days. At the end of the treatment, leaves were randomly sampled from three chilled seedlings (CS) and three non-chilled seedlings (NS) at the same time. Leaves were sampled for the two groups, i.e. chilled plant at budding stage (CF) and non-chilled plant at budding stage (NF). The collected tissues of the four groups (three biological replicates per group) were frozen in liquid nitrogen immediately after harvest and stored at −80°C. For Virus-Induced Gene Silencing (VIGS) assay, three-week-old plants of N. benthamiana grown under normal conditions were used for Agrobacterium infiltration.

2.2. RNA isolation and transcriptome sequencing

Total RNA extraction was extracted using a TriZol reagent (Invitrogen) following the manufacturer’s instructions. The RNA quantity and quality were verified using a NanoDrop 2000 spectrophotometer (Thermo Fisher). For each sample, a paired-end library was prepared using the Illumina TruSeq RNA Sample Prep Kit (Illumina) and sequenced on a HiSeq 2000 platform (Illumina, San Diego, CA, USA). FastQC program was used to assess the quality of sequencing reads. Adaptor sequences were removed from the raw reads and only those with no unknown bases and with low-quality bases (Q ≤ 20) less than 50% were retained as clean reads.

2.3. Differential expression analysis

The newest tobacco genome (var. K326; PI552505) sequences were obtained from the Sol Genomics Network database (https://solgenomics.net/). Clean reads of three biological replicates for CS, NS, CF, and NF were analyzed, respectively. For each data set, paired-end reads were aligned to the tobacco reference genome using HISAT2 (Kim et al. 2015). Read counts were calculated using featureCounts implemented in subread with parameters ‘-Q 10 -B -C’ (Liao et al. 2014). Transcript abundance was estimated using the TPM (Transcripts per million reads) method. Pearson’s correlation coefficients were calculated to assess the agreement of expression data between biological replicates of each group. Principal component analysis (PCA) was conducted to evaluate the variance of gene expression patterns between groups using R software (https://www.r-project.org/). Cluster analysis was performed and hierarchical clustering heat maps were generated with the R package ‘pheatmap.’ Differentially expressed genes (DEGs) between groups were detected using DEGseq2 (Love et al. 2014) and genes with adjusted P < 0.05 were regarded as expressed differentially. Genes with FDR<0.05 were regarded as expressed differentially. Notably, raw P values were used for the detection of DEGs between CF and NF to avoid missing weak but effective information at this development stage.

2.4. Enrichment analysis

Functional descriptions and gene ontology (GO) terms were assigned to the DEGs based on the annotations of the tobacco reference genome (Edwards et al. 2017). Kyoto encyclopedia of genes and genomes (KEGG) pathways were obtained for the DEGs by the online KEGG Automatic Annotation Server (Moriya et al. 2007). Enrichment analysis for GO terms and KEGG pathways were completed using hypergeometric tests in the R package ‘clusterProfiler’ (Yu et al. 2012) followed by FDR correction.

2.5. Identification of flowering-related genes

Arabidopsis gene sequences were downloaded from the TAIR server (https://www.arabidopsis.org/tools/bulk/sequences/index.jsp). The flowering-related genes of Arabidopsis were identified according to the reported flowering regulatory networks (Fornara et al. 2010; Srikanth and Schmid 2011; Walworth et al. 2016; Song and Chen 2018). To analyze the flowering-related genes in tobacco, the tobacco gene sequences were first subjected to local BLASTp similarity search against the Arabidopsis genes. The resultant tobacco genes showed best-hits to the Arabidopsis genes and with E-values lower than 10⁻⁵ and minimal identities of 40 were considered homologous to the corresponding Arabidopsis genes.

2.6. Validation of DEGs by quantitative real-time PCR (qRT-PCR)

The reliability of DEGs identified through RNA-seq was evaluated through qRT-PCR of ten candidate genes. Primers were designed using primer express 3.0 software. The nucleotide sequence of the candidate genes is based on a tobacco database (Nicotiana tabacum v4.5, http://solgenomics.net/). The reverse transcription of RNA to cDNA was performed using SuperScript II reverse transcriptase (Invitrogen, Carlsbad, CA, USA) following the manufacturer’s protocol. The qRT-PCR analyses were performed on a LightCycler 96 real-time PCR system (Roche) using the SYBR® Green Realtime PCR Master Mix (Roche). Relative transcript abundance was calculated using the 2^–ΔΔCt method as previously described (Livak and Schmittgen 2001). The gene L25 was used as an internal control. Primers pairs are listed in Supplementary Table S1.

2.7. Plasmid constructions and virus-induced gene silencing

To construct the vector for the VIGS assay, approximately 300 bp PCR products were cloned using PrimeSTAR® GXL DNA Polymerase (Takara). Plasmid TRV2-LIC (hereafter designated as TRV2, Dong et al. 2007) was digested with Pst I, and then, purified DNA fragment for each gene was inserted into the linearized plasmid using In-Fusion® HD Cloning Kit (Clontech). After sequencing verification, TRV1, TRV2, and the constructed plasmids were introduced into Agrobacterium tumefaciens strain GV3101 (pMP90). Primer pairs used for plasmid constructions are listed in Supplementary Table 1.

The Agrobacterium strain GV3101 containing TRV1, TRV2, and TRV2 derivatives were grown in LB liquid medium supplemented with antibiotics (Kanamycin and Rifampin) at 28°C. After centrifugation, cells were re-suspended with infiltration buffer (10 mM MES, 10 mM MgCl₂, 250 mM acetosyringone), adjusted to a final absorbance (OD₆₀₀= 1.0), and incubated for 2–3 h at room temperature. N. benthamiana leaf infiltration was performed as previously described by Liu et al (2002). Briefly, each Agrobacterium strain containing TRV1 and TRV2 derivatives was well mixed in a 1: 1 ratio, and leaf infiltration was performed with a needleless syringe in three-week-old plants of N. benthamiana. qRT-PCR was performed to determine the expression level of the target gene in the VIGS plants. For chilling stress, WT, TRV2 (Mock), and TRV2-XTHs plants were subjected to low temperature (12°C), followed by transferring them to normal conditions. The expression level of each gene in its VIGS plants was monitored by qRT-PCR. Bud emergence time was recorded for each plant. Primers used in this assay are listed in Supplementary Table S1.

3. Results

3.1. Effects of chilling stress on the flowering in tobacco

It is well established that chilling stress at the seedling stage often leads to early flowering in tobacco. Four-leaf-stage seedlings of NC82 were subjected to chilling stress for 10 days at 12℃. After chilling treatment, the chilled and non-chilled seedlings were grown under normal conditions. The bud emergence of chilling seedlings was 9 d earlier than that in the non-chilling seedlings (Figure 1(A,B)). Furthermore, other phenotypic differences were also observed between chilled and non-chilled plants at the budding stage, with the chilled plants exhibiting significantly less leaf numbers (Figure 1(C)). Our data revealed that NC82 is a preferable cultivar to investigate the effect of chilling stress on the timing of flowering.

Figure 1. The switch to reproductive growth in chilled and non-chilled tobacco plants. (A) Flower bud emergence (arrowed) in plants subjected to chilling or non-chilling stress at the seedling stage; (B) The effect of chilling stress on the timing of bud emergence. Bud emerging day indicates days from the end of chilling treatment to bud emerging; (C) The effect of chilling treatment on the leaf number. Significant differences were determined by the Student’s t-test. **P < 0.01.

3.2. Quality validation of RNA-seq data

To investigate the transcriptomic changes that occur during the process of chilling-induced early flowering in tobacco, four-leaf-stage seedlings were treated with chilling stress, followed by growing under normal conditions (Figure 2(A)). Twelve separate RNA-seq libraries were constructed and sequenced using an illumina sequencing platform. The number of raw reads yielded by RNA-seq ranged from 60314020 to 86155494. After removing low-quality sequences, adaptor sequences, and sequence contaminants, a total of 0.84 billion clean reads (126.06 Gb) from 12 libraries was obtained. On average, the percentages of GC content, Q20, and Q30 were 43.54%, 96.23, and 90.23%, respectively. The proportion of total clean reads mapped to the tobacco reference genome ranged from 93.48% to 94.55% (Supplementary Table S2). Pearson correlation analysis revealed that these results exhibited a good correlation between the three biological replicates of different samples (Figure 2(B)).

Figure 2. Quality assessment of RNA-seq data. (A) Overview of the experimental design. (B) Pearson correlation between the 12 samples based on expression level. (C) Principal component analysis performed on the 12 samples based on the gene expression profile. (D) Verification of the RNA-seq expression by qRT-PCR. L25 was used as an internal control.

Principal component analysis (PCA) showed that three biological replications of each sample were clearly separated. Temporally, samples were well separated along the PC1 axis with NS and CS on the left, NF and CF on the right. We also found perfect separation between NS and CS, suggesting a dramatic transcriptional alteration in the chilled tobacco seedlings. Interestingly, NF and CF clustered near together (Figure 2(C)), suggesting that the differentiation of gene expression induced by chilling at the seedling stage had largely been eliminated at the budding stage. Additionally, the expression of 10 randomly selected genes differentially expressed under chilling stress was evaluated by qRT-PCR. The relative expression levels calculated from qRT-PCR correlated with the log2foldchange values determined by RNA-seq analysis. All 10 genes showed a similar expression pattern in the qRT-PCR and RNA-Seq analyses (Figure 2(D)). Collectively, these results demonstrated the reliability of the RNA-seq data and provided support for the subsequent analysis.

3.3. Comprehensive profiling of gene expression analysis at seedling and budding stages

The results showed that the gene expression profiles in chilled plants were significantly changed at the seedling stage, with 22, 828 genes being differentially expressed (Figure 3(A); Supplementary Table S3). Among these differentially expressed genes (DEGs), 11,600 genes were up-regulated and 11,228 genes were down-regulated (Figure 3(B)). A statistical test showed that the number of CF vs. NF DEGs was significantly lower than those of CS vs. NS DEGs (X² test, P < 0.05). Only 1195 DEGs were at the budding stage, with 979 up-regulated genes and 216 down-regulated genes in CF (Figure 3(C)). A further comparison found that 660 DEGs showed differential expression at both the seedling and budding stages (Figure 3(C), Table 1), among which 209 DEGs were consistently up- or down-regulated.

Figure 3. Gene expression change between chilled and non-chilled plants at the seedling and budding stages. (A) Heatmap visualization of global expression profile for each sample. Columns and rows in the heatmap represent samples and genes, respectively. Color scale indicates fold changes of gene expression. (B) The number of differentially expressed genes at the two development stages. (C) Venn diagrams show the common and specific DEGs between two developmental stages.

Table 1. Common and unique DEGs in chilled tobacco at the seedling and budding stages.

Catalogue	Expression pattern				Total
CS vs. NS	up	down	down	Up
CF vs. NF	up	down	up	Down
	155	54	397	54	660

3.4. Functional classification of DEGs

To understand the gene networks of annotated DEGs for chilled seedlings at the seedling stage, overrepresented GO terms and KEGG pathways were grouped and visualized. The GO enrichment analysis showed that the DEGs were overrepresented in 107 GO terms (FDR < 0.05), such as DNA replication, translation, amino acid metabolism, fatty acid biosynthesis, protein modification, RNA processing, methylation, photosynthesis, response to hormones, antioxidation, and DNA repair (Figure 4(A); Supplementary Table S4). Notably, when we analyzed separately the up-regulated DEGs, which were still enriched in multiple processes related to DNA replication, the regulation of transcription and translation, RNA processing, protein folding/methylation/dimerization, response to hormones, glucose and fatty acid metabolism, antioxidation and DNA repair, while down-regulated DEGs were enriched in GO terms associated with photosynthesis, protein glycosylation/ubiquitination, signal transduction, and cell cycle arrest. KEGG enrichment revealed that the DEGs were overrepresented in 44 pathways, including 9 pathways in the genetic information processing category, 4 pathways in the translation category, 4 pathways in the amino acid metabolism category, 4 pathways in the energy metabolism category, 2 pathways in carbohydrate metabolism, 2 pathways in nucleotide metabolism, and 1 pathway in the environmental adaptation category (Figure 4(B); Supplementary Table S5).

Figure 4. Enrichment analysis of differentially expressed genes in chilled tobacco at the seedling and budding stages. (A) and (C) show enriched GO terms of DEGs in CS vs. NS and CF vs. NF, respectively. (B) and (D) show the KEGG pathways for DEGs in CS vs. NS and CF vs. NF, respectively.

At the budding stage, DEGs in chilled plants were enriched in 57 GO terms (FDR < 0.05), which were mostly related to stress tolerance, energy metabolism, signal transduction, and sexual reproduction (Figure 4(C); Supplementary Table S4). KEGG analysis revealed that the DEGs were significantly enriched in 6 pathways (FDR < 0.05) which were all related to metabolism (e.g. lipid metabolism and amino acid metabolism) (Figure 4(D); Supplementary Table S5). Collectively, these results will facilitate our understanding of the response mechanism to chilling in tobacco at different development stages.

3.5. Effect of chilling on potential flowering-related genes

To comprehensively investigate the candidate genes implicated in flowering regulation in tobacco, a local BLAST similarity search was performed against the Arabidopsis genes, and putative tobacco flowering-related genes in tobacco were identified accordingly. About 458 tobacco gene sequences were found homologous to the 124 Arabidopsis flowering pathway genes. Of these potential flowering-related genes, 212 homologs showed differential expression in CS vs. NS, including 27 vernalization pathway genes, 106 photoperiod/circadian clock pathway genes, 5 autonomous pathway genes, 8 GA pathway genes, 3 ageing pathway genes, 4 integrators, and 16 genes involved in flower development or other related processes, according to the reported flowering pathways and regulatory networks (Supplementary Tables S6 and S7).

For identified flowering genes in the vernalization pathway, the putative flowering activators, including FIE1, MSI1, and CSTF64 were up-regulated, while repressors FRI and HUA2 were down-regulated (Figure 5). In the GA pathway, the DELLA gene GAI which represses flowering was down-regulated, while GID1B and GID1C that represses DELLA proteins to activate flowering were up-regulated. In the autonomous pathway, putative flowering activators, including MSI4/FVE, FY, and FLK were up-regulated. The main rhythm-related genes LHY, ELF4, TOC1, and transcription factor LUX, were all up-regulated in CS. Photoperiod-regulated flowering activators GI, ADO3/FKF1, CRY2, and transcription factor bHLH63/CIB1 were up-regulated while repressor TEM1 was down-regulated. The age-regulated flowering activator SPL5 was up-regulated. Flowering locus T (FT) is a primary integrator in the regulation of plant flowering. In tobacco (Nicotiana tabacum), five FT-like genes have been functionally identified, three of which (NtFT1–NtFT3) encode floral repressors, whereas NtFT4 and NtFT5 encode floral activators (Harig et al. 2012; Beinecke et al. 2018). We found that only NtFT5 (Nitab4.5_0001003g0220.1) was significantly up-regulated by chilling stress at the seedling stage (log2foldchange ≥1, P ≤ 0.05). As antagonistic FTs (NtFT1–NtFT3) and another floral activator NtFT4 were predominantly expressed under short-day conditions (Beinecke et al. 2018), and these genes showed no statistically significant difference between the chilled and non-chilled plants (Supplementary Table S8), suggesting that chilling might not affect this photoperiod-dependent expression pattern. Putative flowering repressor SVP, which responds to ambient temperature and represses FT to repress flowering, was down-regulated. The AP2/ERF and B3 domain-containing transcription factor RAV1, a putative negative regulator of flower development, was down-regulated. Several other genes involved in the flower development (e.g. DDL, EMF2, and LDL2) were up-regulated in chilled tobacco seedlings.

Figure 5. Response of major flowering-related genes in chilled tobacco at the seedling stage. Genes in blue, green, gray, purple, and cyan are ‘vernalization,’ ‘circadian clock/photoperiod,’ ‘autonomous,’ ‘GA,’ and ‘age’ pathway genes, respectively. Genes in black are the flowering pathway integrators. Frames in red, green, black, and dashed lines indicate up-regulated, down-regulated, unchanged, and not-identified genes, respectively. Frames in golden indicate genes with both up-regulated and down-regulated homologs. Arrows indicate positive regulation and bars indicate negative regulation. The relationships among the listed genes are drawn according to the diagram for A. thaliana by Fornara et al. (2010), Srikanth and Schmid (2011).

However, the putative vernalization-regulated flowering activators CLF, REF6, and FPF1 were down-regulated while the repressor PEP was up-regulated (Figure 5). A similar situation was also observed for photoperiod/circadian clock and autonomous pathway genes, in which the putative repressors PHYB, COP,1 and SPA were up-regulated while the activators FLD, SPL1, and SPL3 were down-regulated. Flowering pathway integrator AGL24 was down-regulated and LFY remained unchanged. In addition, homologous genes to several Arabidopsis key flowering pathway genes were not found in tobacco, including FL, CO, VRN2, CDF1, GID1A, ARP6, DNF, SMZ, SNZ, TEM2, and TOE3. These results suggested that although many of the known genetic flowering pathway genes shared a considerable degree of conservation in tobacco and other plants, a proportion of critical flowering genes might play different roles or have other function substitutes in tobacco.

3.6. Functional identification genes involved in chilling-induced early flower

Previous analysis revealed that 660 DEGs of CS vs. NS were also differentially expressed in CF vs. NF, among which 397 DEGs were down-regulated and up-regulated in CS and CF, respectively (Table 1, Supplementary Table S9). It is interesting to note that 13 genes that encoded xyloglucan endotransglucosylase/hydrolase (XTH) showed a contrary expression pattern in the chilled tobacco, with down-regulation at the seedling stage and up-regulation at the budding stage. Among these genes, the expression levels of eight genes were significantly changed (Figure 6(A), Supplementary Table S10), leading us to postulate that XTHs probably play roles in the process of chilling-induced early flowering. Three XTH genes (Nitab4.5_0003734g0030.1, henceforth referred as to NtXTH22; Nitab4.5_0003353g0040.1, NtXTH6; Nitab4.5_0002630g0010.1, NtXTH23), were selected for functional analysis. As NtXTH6, NtXTH22, and NtXTH23 shared 90.7%, 96.79%, and 89.86% identities with NbXTH6, NbXTH22, and NbXTH23, respectively, the VIGS assay was employed to identify the role of XTHs in the chilling-induced early flowering. When TRV2-PDS plants showed a bleaching phenotype, the expression of XTH in the corresponding TRV2-XTH plants was determined by qRT-PCR (Supplementary Figure S1). Under normal conditions, the flowering time of TRV2-NbXTH6, TRV2-NbXTH22, and TRV2-NbXTH23 plants showed no significant difference compared with WT and TRV2 plants. When plants were subjected to chilling stress, TRV2-NbXTH22 plants flowered later than TRV2 and WT plants (Figure 6(B)), and the expression level of NbXTH22 showed a certain negative correlation with flowering time (Supplementary Figure S1D). The flowering time of TRV2-NbXTH6 and TRV2-NbXTH23 plants was comparable with those in TRV2 and WT (Figure 6(B)). Taken together, our results indicated that the down-regulation of NbXTH22 can alleviate, to some extent, the effect of chilling-induced early flowering.

Figure 6. Functional analysis of NtXTHs in chilling-induced early flowering. (A) Heatmap visualization of the expression profile of NtXTHs. (B) Comparison of bud emerging time between WT, TRV2, TRV2-NtXTH6, TRV2-NtXTH22, TRV2-NtXTH23 plants. Plants with or without the downregulation of NtXTHs were subjected to chilling stress for 5 days, and then, returned to normal conditions. * indicates the significant difference, P < 0.05. Data were analyzed using the Student’s t-test.

4. Discussion

4.1. Potential flowering-related genes and regulatory networks in tobacco under chilling stress

In the vernalization pathway of Arabidopsis, the central flowering repressor gene FLC modulates downstream signaling, and all other genes in the vernalization pathway affect the flowering time by altering the expression of FLC (Srikanth and Schmid 2011). Although functional FLC orthologs are not identified in the tobacco genome (Sierro et al. 2014; Edwards et al. 2017), several genes potentially involved in the vernalization pathway were identified to be regulated in CS in the present study. Consistent with previous studies, the key regulator FRI, a plant-specific gene that activates the expression of FLC in Arabidopsis (Geraldo et al. 2009), was down-regulated in CS. The transcription factor HUA2, which functions as a repressor of flowering by enhancing the expression of several genes that delay flowering, was also down-regulated in CS. Meanwhile, the FIE1, MSI1, and CSTF64 genes required for the repression of FLC (De Lucia et al. 2008; Liu et al. 2010; Wood et al. 2016) were up-regulated in CS (Figure 5, Supplementary Table S6). Therefore, it is reasonable to infer that the flowering transition of tobacco derived by chilling stress could probably proceed through another factor that is functionally analogous to FLC.

Plants having a vernalization requirement adjust their flowering time according to the changes in both temperature and photoperiod. Components of the photoperiod pathway regulate flowering time in tight interaction with components of the circadian clock. We observed that vernalization regulated some genes involved in the circadian clock and photoperiod flowering pathway. GI, a component of the circadian clock and photoperiod pathway, is a key regulator of flowering time by increasing the expression of FT through miRNA172 (Sawa and Kay 2011). CO promotes flowering under long day conditions through the positive regulation of the floral pathway integrator FT. FKF1 is also a component of the circadian clock and the photoperiod flowering pathway promoter, which interacts with GI and activates CO (Ream et al. 2014). Intriguingly, similar to FLC as described above, both GI and FKF1 were increased in CS in this study (Figure 5, Supplementary Table S6), consistent with previous studies reporting the up-regulation of GI during cold treatment. Although CO orthologs are not identified in the tobacco genome (Edwards et al. 2017), the expression levels of several putative CONSTANS-LIKE (COL) genes were affected by chilling stress, with COL4 and COL9 being significantly down-regulated and up-regulated, respectively (Supplementary Table S7). Additionally, several other key circadian clock protein genes (e.g. TOC1, CCA1/LHY, and PRRs) and photoperiod pathway genes (e.g. PHYB, CRY1, CRY2, COP1, and SPA) involved in the regulation of CO were also regulated in CS. These results suggested connections between the photoperiod and vernalization mediated flowering pathways in tobacco, and the differences of certain key regulators between tobacco and other plants such as Arabidopsis.

4.2. Long-term memory of chilling-induced DEGs from the seedling stage to the budding stage

As described above, gene abundance between chilled and non-chilled tobacco changed dramatically at the seedling stage while tending to similar patterns when developed to the budding stage (Figure 2, Table 1), suggesting that most of the responses to vernalization occurred immediately after vernalization. Many flowering-related genes were deregulated in chilled seedlings (Supplementary Table S3), which may indicate that the required genetic prerequisites for floral initiation at a later stage during development are established at the seedling stage. Nevertheless, a proportion of genes were consistently regulated in their response to vernalization between the two development stages (Table 1). Of interest in the class of genes that are consistently up-regulated in chilled plants at both stages are those related to flowering, such as SPT, NAC07, DOG1, NPR5, and SPL9 (Alvarez and Smyth 1999; Hepworth et al. 2005; Xing et al. 2010; Pitaksaringkarn et al. 2014; Huo et al. 2016).

Except for genes involved in vernalization-induced flowering processes, many differentially expressed genes would be simply related to cold response. Consistent with the previous envisage that many cold-responsive genes would be rapidly regulated to vernalization, we did find an enrichment of stress response genes in CS (Supplementary Table S9). Noticeably, a proportion of these genes were similarly regulated (mostly up-regulated) at the two stages. Polyamines (PAs), mainly including putrescine, spermindine, and spermine, are important regulators for plant growth and development and closely related to stress resistance in plants. Previous studies demonstrated a positive relationship between PAs and cold tolerance (Guye et al. 1986, 1987; Kushad and Yelenosky 1987). ADC1 and SAMDC are two essential components that are required for the biosynthesis of PAs (Ge et al. 2005). These two genes were found consistently higher expressed, suggesting that PAs are essential for (cold) stress resistance in tobacco. Other consistently regulated genes response to stress include serine/threonine-protein kinase BSK5, hydrophobic proteins RCI2B, SNL6, NAC002, SCL14, heat shock factor proteins HSF24, AFP2, and PLP6. These genes might not only contribute to the rapid regulation of vernalization but also continue to play important roles in the development stage of plants, suggesting that plants might have developed intricate mechanisms to better cope with the period of and after chilling stress.

4.3. XTHs may plant essential roles in chilling-induced early flowering

Xyloglucan is a hemicellulose that is mainly present in the primary cell walls of land plants. The major role of XTH proteins in plants is to cleave and rearrange the backbone of xyloglucan, and thus, to regulate the composition and organization of the cell wall. Besides the reported role of XTHs in cell wall extensibility during primary root elongation, wood formation, hypocotyl growth, and vein differentiation (Matsui et al. 2005; Wu et al. 2005; Osato et al. 2006; Nishikubo et al. 2011), XTHs have been demonstrated to be involved in flower opening in Sandersonia, Alstroemeria, Dianthus caryophyllus, and Gerbera hybrid (O’Donoghue et al. 2002; Breeze et al. 2004; Laitinen et al. 2007; Harada et al. 2011). Tobacco (Nicotiana tabacum L.) harbors 56 XTH genes, which can be classified into two subfamilies (Wang et al. 2018). Among 397 DEGs that showed down-regulation at the seedling stage and up-regulation at the budding stage in chilled plants, it is interesting to note that the XTH gene family harbors the largest number of DEGs than the other gene families. We found that eight XTH genes were significantly down-regulated under chilling stress at the seedling stage and up-regulated at the budding stage in the chilled plant (Figure 6(A)). Further functional analysis revealed that the down-regulation of XTH22 can alleviate the early flowering phenotype caused by chilling stress (Figure 6(B)), suggesting the involvement of XTHs in the chilling-induced early flowering. To date, no literature reported the correlation between XTH activity and chilling-induced early flowering. In further study, the mechanism of XTH underlying chilling-induced early flowering remains to be elucidated.

5. Conclusion

In this study, we constructed a regulatory network of potential flowering-related genes under chilling stress, and revealed a long-term memory of chilling-induced genes from the seedling stage to the budding stage. We also identified NtXTH22 as a candidate gene to alleviate the negative effect caused by chilling stress at the seedling stage. To our knowledge, this is the first report on the genome-wide characterization of transcriptome profiles for chilling-induced early flowering in tobacco. Although additional experiments are necessary to validate the proposed genes and regulatory models, the present study provides valuable resources of candidate genes and is fundamental to further understand the complex networks that regulate the floral transition in tobacco.

Author contributions

Formal analysis: W.G., J.J. and Y.X.; investigation: G.X., Z.L. and H.Z.; validation: G.X., C.W., and Z.L.; data curation: S.D.; original draft preparation: G.X. and W.G.; editing: H.Z. and S.D.; funding acquisition: G.X. and H.Z. All authors have read and agreed to the published version of the manuscript.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability

The data used to support the findings of this study are available from the corresponding author on reasonable request. All transcriptome datasets generated in this study have been deposited in the Genome Sequence Archive (Genomics, Proteomics & Bioinformatics 2017) in BIG Data Center (accession number CRA002436, https://bigd.big.ac.cn/gsa).

Additional information

Funding

This research was funded by grants from the Zhengzhou Tobacco Research Institute (grant number 902018CA0280), the Key Scientific and Technological Project of Henan Province (grant number 212102110446), and the Project of Tobacco Genome (grant number 110202001026 (JY-09)).

Notes on contributors

Guoyun Xu

Guoyun Xu is an associate research fellow in the Zhengzhou Tobacco Research Institute of CNTC.

Wuxia Guo

Wuxia Guo is a research assistant in the South China Botanical Garden, Chinese Academy of Sciences.

Zunqiang Li

Zunqiang Li is a researcher in Mudanjiang Tobacco Science Research Institute.

Chen Wang

Chen Wang is a researcher in the Zhengzhou Tobacco Research Institute of CNTC.

Yalong Xu

Yalong Xu is a senior engineer in the Zhengzhou Tobacco Research Institute of CNTC.

Jingjing Jin

Jingjing Jin is an associate research fellow in the Zhengzhou Tobacco Research Institute of CNTC.

Huina Zhou

Huina Zhou is a professor in the Zhengzhou Tobacco Research Institute of CNTC.

Shulin Deng

Shulin Deng is a professor in the South China Botanical Garden, Chinese Academy of Sciences.