Skip to main content
Download PDF

Transcriptome-based gene regulatory network analyses of differential cold tolerance of two tobacco cultivars

BMC Plant Biology volume 22, Article number: 369 (2022) Cite this article

Abstract

Background

Cold is one of the main abiotic stresses that severely affect plant growth and development, and crop productivity as well. Transcriptional changes during cold stress have already been intensively studied in various plant species. However, the gene networks involved in the regulation of differential cold tolerance between tobacco varieties with contrasting cold resistance are quite limited.

Results

Here, we conducted multiple time-point transcriptomic analyses using Tai tobacco (TT, cold susceptibility) and Yan tobacco (YT, cold resistance) with contrasting cold responses. We identified similar DEGs in both cultivars after comparing with the corresponding control (without cold treatment), which were mainly involved in response to abiotic stimuli, metabolic processes, kinase activities. Through comparison of the two cultivars at each time point, in contrast to TT, YT had higher expression levels of the genes responsible for environmental stresses. By applying Weighted Gene Co-Expression Network Analysis (WGCNA), we identified two main modules: the pink module was similar while the brown module was distinct between the two cultivars. Moreover, we obtained 100 hub genes, including 11 important transcription factors (TFs) potentially involved in cold stress, 3 key TFs in the brown module and 8 key TFs in the pink module. More importantly, according to the genetic regulatory networks (GRNs) between TFs and other genes or TFs by using GENIE3, we identified 3 TFs (ABI3/VP1, ARR-B and WRKY) mainly functioning in differential cold responses between two cultivars, and 3 key TFs (GRAS, AP2-EREBP and C2H2) primarily involved in cold responses.

Conclusion

Collectively, our study provides valuable resources for transcriptome- based gene network studies of cold responses in tobacco. It helps to reveal how key cold responsive TFs or other genes are regulated through network. It also helps to identify the potential key cold responsive genes for the genetic manipulation of tobacco cultivars with enhanced cold tolerance in the future.

Peer Review reports

Background

Sessile plants are constantly subjected to various abiotic clues during the lifecycle such as drought, salinity, heat and cold. These stresses dramatically limit their normal growth and development, and final productivity as well [1]. Naturally occurring cold stress usually causes complex symptoms including growth retardation, leaf chlorosis, wilt and droop, chilling/frost/freeze injury and final death of tissues or plants [2,3,4]. Therefore, cold stress is considered as one of the most destructive threats that adversely affect plant lifecycle [5, 6]. Plants have to adjust multilayer transient changes, including physiological, biochemical, metabolic and molecular alterations, to adapt cold stress [4, 7,8,9,10,11,12].

Gene expression change is one of adaptive molecular mechanisms for plants responding to cold stress. In particular, rapid up-regulation of cold-responsive genes, such as C-REPEAT BINDING FACTOR (CBF) genes, COR (Cold Regulated) genes, and other key regulons can enhance cold tolerance [7, 12,13,14,15,16,17,18]. With the advent of high throughput sequencing, transcriptomic studies have been widely applied in many plant species, including model plants [19], crops, horticulture plants [20,21,22,23], vegetables [24,25,26,27,28], trees [29, 30] and tobacco [31, 32]. These studies provide insights into how transcriptional changes occur during cold stress.

Tobacco (Nicotiana tabacum), originally, is a tropical and subtropical plant. It currently becomes commercially popular economic crops grown worldwide. Moreover, tobacco leaves and BY-2 cells provide an excellent system for studies of functional genomics and metabolism [33, 34]. However, tobacco is overall highly susceptible to cold stress [35,36,37,38]. Mechanisms underlying cold responses in tobacco have already been investigated at physiological, biochemical, transcriptional, and metabolic levels. Cold acclimation helps tobacco to cope with cold stress through improving antioxidant enzymatic activities, osmotic adjustment potentials and photochemical efficiency [39], these changes are partially caused by reprogramming gene expression and metabolites [31]. Cold treatment can change accumulation levels of certain functional proteins involved in photosynthesis, protein/RNA procession and redox etc. [40]. NtbHLH123 acts as a positive regulator to activate expression of cold responsive genes, NtCBF (C-repeat binding factor) and alleviate damage of ROS (reactive oxygen species), thereby enhancing cold resistance [41]. Moreover, functional studies showed that tobacco plants with overexpression of chloroplast ω − 3 fatty acid desaturase gene [36], TaWRKY10 [42], Capsella bursa-pastoris CE53 and CBF genes (CblCE53 and CbCBF) [43] exhibit enhanced cold tolerance, indicating that these genes play positive roles in tobacco cold responses. Thus, elucidating the underlying mechanisms responsible for cold responses, especially identifying key cold-responsive genes, can facilitate the development of genetically modified tobacco varieties with cold tolerance. However, the gene networks related to cold responses are still understudied in tobacco.

Results

Identification of differentially expressed genes between two cultivars under cold treatment

Yan tobacco (YT) and Tai tobacco (TT) are two important tobacco cultivars for tobacco production, but exhibit contrasting cold responses. YT is cold resistant whereas TT is susceptible to cold stress. To investigate transcriptional changes between two cultivars under cold stress, we conducted RNA-seq with 6 time points before (0 h, CK) and after cold treatment at 4 °C (0.5 h, 1 h, 2 h, 4 h and 8 h), three biological replicates per sample (Additional file 2: Table S1). Correlation analyses showed that replicated data were generally well correlated within the same time point in YT (Additional file 1: Fig. S1A) or TT (Additional file 1: Fig. S1B).

To examine time-point-related transcriptional changes in YT and TT after cold treatment, we conducted a Principal Component Analysis (PCA), and observed that PC2 was quite different between TT and YT (Fig. 1A). YT and TT had a similar trend of time-point-related expression changes after cold treatment, slight expression changes occurred before 1 h, whereas dramatic expression changes occurred after 2 h (Fig. 1A). We then identified differentially expressed genes (DEGs) occurred in TY and TT at each time point after exposure to cold treatment compared with CK (Additional file 2: Table S2). We found that TT had more DGEs, especially for down-regulated genes, than YT at the same time point, suggesting that more genes are dynamic in TT relative to YT during cold stress. To assess cold-induced transcriptional changes in YT or TT, we counted total cold-induced DEGs, one of which differentially expressed at least at one time point, in YT (5622 DEGs) and TT (6799 DEGs). We then standardized the expression of DEGs by removing some extreme expression values and generated a heatmap with 10 subclusters to visualize time point related changes in YT (5602 DEGs, Fig. 1B) or TT (6774 DEGs, Fig. 1C) and between YT and TT. Time point related similar and distinct changes were observed between YT and TT in response to cold stimulus. After conducting pair-wise comparisons, we found that 1169 and 2418 genes were down- and up-regulated in TT and YT, 337 genes (165 genes up-regulated in TT but down-regulated in YT, and 172 genes down-regulated in TT but up-regulated in YT) exhibited opposing expression changes between TT and TY, and the rest of up- or down-regulated genes were specific to TT or YT (Fig. 1D).

Fig. 1
figure 1

Characterization of differentially expressed genes before and after cold treatment in TT or YT. A Principal Component Analysis (PCA) of gene expression profiles in TT and YT before and after cold treatment. Each point represents one treatment for gene expression profiling. B, C Heatmaps showing genes differentially expressed at each time point in YT (B) or TT (C) before (0 h) and after cold treatment (0.5, 1, 2, 4, 8 h). 5602 DEGs in YT and 6774 DEGs in TT were analyzed by using k-means clustering. Each heatmap was generated through k-means clustering algorithm with the Z-score result of gene expression in Cluster 3.0 software. The color key representing the standardized gene expression levels from high (red) to low (blue). D Venn plots showing comparisons of genes that were up-regulated or down-regulated after cold treatment between YT and TT

Full size image

To assess if DEGs are related to any significant biological relevance, we conducted GO term enrichment analyses using up- and down-regulated genes. We found that up-regulated genes had much more GO terms than down-regulated genes (Additional file 1: Fig. S2). Moreover, we found that GO terms differentially occurred in up-regulated genes specific to TT or YT. GO term functions in signal/molecular transduce activity, responses to stimuli, receptor/binding activity and (cellular) protein/primary metabolic processes were significantly enriched in the up-regulated genes in TT instead of YT. In contrast, functions in thylakoid, photosynthesis, membrane, extracellular region and generation of precursor metabolites and energy were overrepresented in the up-regulated genes in YT but not in TT (Additional file 1: Fig. S2). Furthermore, GO terms differentially occurred in genes up-regulated in TT but down-regulated in YT and genes down-regulated in TT but up-regulated in YT (Additional file 1: Fig. S2).

To further assess cold-induced transcriptional changes between YT and TT, we identified DEGs occurred between YT and TT at each time point (Fig. 2A). We noticed that there were 857 YT-higher genes and 750 TT-higher genes before cold treatment (CK), reflecting genotype-related gene expression changes. As shown in Fig. 2B, 4,448 DEGs between YT and TT were grouped into 10 subclusters for showing genotype-related gene expression changes at each time point. For example, DEGs in C6 and C7 were highly expressed in YT and TT, respectively, at almost each time point. To assess if the DEGs in each cluster have any significantly functional enrichment, we conducted GO term enrichment analyses. We observed variations in GO terms among 10 subclusters (Additional file 1: Fig. S3). DEGs in C10 exhibited highly enriched GO functions related to metabolic processes, membrane and thylakoid; GO functions of DEGs in C6 and C9 were overrepresented in response to stresses. Taken together, all above analyses indicate DEGs display variations in the expression profiles and GO terms functions at each time point between TT and YT, thereby resulting in changes in response to cold treatment between TT and YT.

Fig. 2
figure 2

Characterization of differentially expressed genes between YT and TT in each time point. A Box plots showing the number of DEGs (TY or TT-higher) occurred at each time point. B Heatmap showing expression profiles of DEGs at each time point for YT vs TT. 4448 DEGs were analyzed by hierarchical clustering. The heatmap was generated through k-means clustering algorithm with the row-wise Z-score result of gene expression in Cluster 3.0 software. 4448 DEGs were grouped into 10 distinct subclusters according to gene expression profiles. The color key representing the standardized gene expression levels from high (red) to low (blue). C Venn plots showing comparison of genes that were up-regulated or down-regulated before (CK) or after cold treatment

Full size image

Compared with CK, we identified 1443 YT-higher and 1561 TT-higher genes under cold treatment (Fig. 2C). According to GO term enrichment analyses, we observed distinct GO terms associated with up- or down-regulated genes in CK and cold treatment (Additional file 1: Fig. S4). In general, up-regulated genes in CK were mainly responsible for some fundamental biological functions such as DNA/protein/nucleotide binding activities and kinase/hydrolase/pyrophosphatase activities, whereas cold-induced up- or down-regulated genes were mainly enriched in some biologically regulatory functions, responses to stresses, photosynthesis-related functions and reproduction. Collectively, transcriptomic analyses show that TT and YT have some specific DEGs functioning in distinct biological processes during cold treatment, which possibly cause differential cold responses between YT and TT.

Identification of differentially expressed TFs and protein kinases under cold treatment

It has been documented that transcription factors (TFs) [44, 45] and various types of protein kinases (PKs) [46,47,48] act as key regulators in response to abiotic/biotic stresses in plants. To assess how TFs and PKs are differentially expressed in TT/YT or between TT and YT under cold treatment, we classified all DEGs in TT/YT or between TT and YT into TFs, PKs and other genes (Additional file 1: Fig. S5). Compared with TT, YT had a similar number of down-regulated TFs and PKs but have less up-regulated TFs and PKs. Moreover, we observed variations in percentage of down- and up-regulated TFs in TT/YT or between TT and YT (Additional file 1: Fig. S6A). For instance, WRKY, MYB, NAC and AP2-EREBP TFs had relatively higher percentage than other up-regulated TFs in TT and YT. GRAS, C2H2 and bHLH TFs had similar percentage in both down- and up-regulated TFs in TT and YT. Similarly, we also observed variations in percentage of down- and up-regulated PKs in TT/YT or between TT and YT (Additional file 1: Fig. S6B). More Leucine Rich Repeat (LRR) Kinases XI&XII were down-regulated than other PKs in TT and YT, while more LRR Kinase VII and III were down-regulated in TT relative to YT.

To examine time-point-dependent expression changes of TFs and PKs under cold treatment, we conducted pairwise comparisons using differentially expressed TFs and PKs in TT and YT (Additional file 1: Fig. S7). For 609 TFs/384 PKs in TT (Additional file 1: Fig. S7A, B) and 526 TFs/303 PKs in YT (Additional file 1: Fig. S7C, D), which were differentially expressed at least at one time point after cold treatment, 89 TFs/79 PKs in TT and 47 TFs/19 PKs in YT were differentially expressed across all time points (Additional file 2: Table S3), other TFs and PKs were differentially expressed at one time point under cold treatment. Similarly, we conducted similar analyses using 347 TFs (Fig. 3A) and 217 PKs (Fig. 3B) differentially expressed at least at one time point between TT and YT under cold treatment. Compared with TT, we found that 14 TFs/6PKs were differentially expressed across all time points under cold treatment (Additional file 2: Table S3).

Fig. 3
figure 3

TFs or PKs differentially expressed between YT and TT under cold treatment. A, B Venn plots showing TFs (A) or PKs (B) differentially expressed at each time point between YT and TT after cold treatment (0.5 h, 1 h, 2 h, 4 h, 8 h). C GO term enrichment analyses of TFs or PKs that differentially expressed at least at one time point between YT and TT before (0 h) and after cold treatment (0.5 h, 1 h, 2 h, 4 h, 8 h)

Full size image

After conducting GO term enrichment analyses, we found that 347 PKs had significantly enriched GO terms with functions in transferase activities, (cellular) protein metabolic processes, macromolecular/cellular protein modifications, kinase activities; as expected, GO terms with functions of 217 TFs were mainly overrepresented in transcription factor activities, nucleic acid/DNA binding, regulation of gene expression and metabolic processes and biological regulations (Fig. 3C). Based on the above analyses, we found that TFs differentially expressed in the two samples after cold treatment mainly included AP2-EREBP, BHLH, WRKY, MYB TF family, and protein kinases mainly included Leucine Rich repeat, Plant External Response Like Kinase and Calcium Dependent Protein Kinase. In addition to some of TFs or protein kinases already reported in cold responses of different plant species [41, 49], we found that some new genes such as Nitab4.5_0002105g0060 (C2H2) (Additional file 2: Table S3) were potentially in relation to cold responses, which are necessary to be functionally validated in tobacco or other plants.

Weighted gene co-expression network analysis (WGCNA)

To further unveil if hub genes or TFs are involved in cold responses through interacting with other genes or TFs in TT and YT, we conducted Weighted Gene Co-Expression Network Analysis (WGCNA) [50] using 9583 DEGs that were differentially expressed at least at one time point in TT or YT after cold treatment. In particular, WGCNA specifically focused on two well correlated biological replicates (Additional file 1: Fig. S8). We obtained 17 co-expression modules visualized using different colors according to the WGCNA package function (Fig. 4A). After extracting the number of genes and TFs in each module, we observed variations in the number of genes or TFs per module, ranging from 38 genes and 2 TFs in the grey60 module to 2258 genes and 233 TFs in the turquoise module with an averaged number as 488 genes and 44 TFs (Additional file 1: Fig. S9). In addition, we constructed module-time point correlation heatmap to clearly visualize the correlation between each module and each time point (Fig. 4B). We also examined the overall expression trend of co-expressed genes in each module at each time point (Additional file 1: Fig. S10). We found that the overall expression levels of genes in the brown module at each time was higher in YT than in TT; Genes in the pink module tended to be up-regulated at the first time point of cold treatment (0.5 h) in YT and TT.

Fig. 4
figure 4

Identification of co-expression gene modules using differentially expressed genes. A Hierarchical cluster tree showing co-expression modules identified by WGCNA. Each leaf in the tree represents one gene. The major tree branches constitute 17 modules labelled by different colors. The x-axis represents genes, and the y-axis represents the co-expression distance. Modules were identified using dynamic tree cutting by dividing the dendrogram at significant branch points. The horizontal bar immediately below the dendrogram represents modules indicated with different colors. B Module-time point correlation. Each row represents a module. Each column corresponds to a specific sample. The color of each cell at the row–column intersection indicates the correlation coefficient between the module and sample. Red indicates a high degree of correlation between a specific module and the sample

Full size image

Identification of key genes involved in cold responses in tobacco

We connected the 17 modules with each time point post cold treatment in TT and YT. According to module-time point association, we found that the correlation coefficient between genes in the pink module and time points was very low before cold treatment (0 h), while became very high at 0.5 h after cold treatment in TT. A similar trend occurred in YT, but the change was relatively weak (Figs. 4B and 5A). Meanwhile, the pink module containing 471 genes, including 80 TFs, was the module with the most TF relative to the other modules (Additional file 1: Fig. S9). To assess roles of genes in the pink module in cold responses, we conducted GO term enrichment analyses, and observed that they were mainly involved in various biological and metabolic processes (Additional file 1: Fig. S11A). GO terms related to the regulation of gene expression, metabolic processes, cellular processes and biological processes and DNA binding were enriched for genes after cold treatment, which is consistent with the possibility that plants adjust changes responding to cold treatment.

Fig. 5
figure 5

Analyses of the pink and brown module. A Module gene expression pattern of the brown module, red indicates high expression level and green indicates low expression level. B The correlation network of the pink module. A gene network is constructed by using WGCNA, in which each node represents a gene; the connecting line (edge) between genes represents the co-expression correlation. The genes with edge weights > 0.1 and sorting the top 50 by the number of edges are visualized by Cytoscape. The size and color of each circle represent the number of edges. C Module gene expression pattern of the brown module, red indicates high expression level and green indicates low expression level. D The correlation network of the pink module. A gene network is constructed by using WGCNA, in which each node represents a gene; the connecting line (edge) between genes represents the co-expression correlation. The genes with edge weights > 0.1 and sorting the top 50 by the number of edges are visualized by Cytoscape. The size and color of each circle represent the number of edges

Full size image

To further assess if genes in the pink module are involved in cold responses through interacting with other genes. We generated a co-expression network containing genes with edge weights > 0.1 and the top 50 number of edges, which were derived from DEGs in the pink module. We obtained the top 50 hub genes in the pink module, which contain 8 TFs and 42 other genes (Fig. 5B, Additional file 2: Table S4). Specific descriptions of hub genes in the network are listed in Additional file 2: Table S4. We indeed found some important cold-related and stress-related genes, such as NtCBF [41] that play important regulatory roles in coping with cold stress [18], other AP2-EREBP, GRAS and C2H2 TF family.

Identification of genes responsible for differential cold responses between TT and YT

After a closer examination, we found that only the brown module exhibited a high association with all time-points in YT instead of TT (Fig. 4B), suggesting that this module may contain genes responsible for cold tolerance in YT. Consistently, we found that genes in the brown module were mainly up-regulated at all time-points in YT as compared to TT (Fig. 5C). To assess if co-expression genes in the brown module are possibly involved in differential cold responses between TT and YT, we performed GO term enrichment analyses, and found that they were mainly involved in responses to stress or stimuli (Additional file 1: Fig. S11B).

Similarly, we conducted a co-expression network containing genes with edge weights > 0.1 and the top 50 number of edges, which were derived from DEGs in the brown module (Fig. 5D). We found that the co-expression network contained some hub genes like RAP74, UDP-Glycosyltransferase, WRKY21, ABI3/VP1 (ABSCISIC ACID INSENSITIVE 3/VIVIPAROUS 1, RAV), PSP (3-phosphoserine phosphatase) and ARR12 (type A response regulator 12) genes (Additional file 2: Table S5). Some of these genes have been found to be involved in stress responses or signal transduction in plants [51,52,53,54,55,56]. Moreover, AtRAP74 is involved in osmatic stress and ABA signaling through interacting with CTD phosphatase-like 3 (CPL3) and CPL4 [57].

Taken together, all above analyses indicate that differential expression of genes in the brown module may cause differential responses to cold treatment between TT and YT.

Roles of hub transcription factors (TFs) in cold stress

Based on WGCNA, we identified some important hub genes that were co-expressed in the same module with each other. We speculate if they have certain regulatory relationships with each other, and they play important roles in the regulation of DEGs under cold treatment. By using tree-based methods [58], we generated inferring regulatory networks to study the regulatory relationship between TFs and DEGs by using 11 hub TFs in the pink or brown module as regulators. We obtained genetic regulatory networks (GRNs) containing 11 TFs through the R package GENIE3 (Additional file 1: Fig. S12, Additional file 2: Table S7). We found that 11 TFs had possible interaction relationships with the aforementioned hub genes, but also with genes within other modules (Additional file 1: Fig. S12), indicating that 11 TFs play vital regulatory roles in responding to cold stress in TT or YT. To visualize the regulatory relationship among these genes more clearly, we exemplified the network only containing the hub genes/TFs of the pink or brown module identified in WGCNA (Fig. 6). We observed clearly regulatory relationships occurred between hub TFs and hub genes. For instance, in the brown module, the three hub TFs (ABI3/VP1, ARR-B and WRKY) exhibited regulatory relationships with hub genes or with each other (Fig. 6). It has been documented that WRKY TF family is involved in cold stress in several plants [59], such as Arabidopsis [60] and eggplant [61]. Moreover, regulatory relationships between TFs and hub genes also occurred in the pink module, including cold-responsive TFs NtCBF2, NtCBF3 and NtCBF5, which were well characterized in tobacco [41]. Some of TFs (GRAS, C2H2 and AP2-EREBP) associated with stress/stimulus responses were also detected in the network and exhibited regulatory relationships with other hub genes and DEGs (Fig. 6; Additional file 1: Fig. S12). In Arabidopsis, GRAS gene is stress-inducible and functions in enhanced cold tolerance via regulating the expression of other related genes [62]. Thus, according to GRN, we identified genes already known to be associated with cold stress in tobacco, such as NtCBF genes; more importantly, we also identified other genes that were involved in cold responses.

Fig. 6
figure 6

Genetic regulatory networks (GRNs) between hub genes. The circle represents gene, the diamond represents TF, and the connecting line is divided into source and arrow. The gene pointed by the arrow represents the gene targeted by source transcription factor. If both ends of the connecting line are arrows, the two TFs regulate each other

Full size image

To verify the accuracy of DEGs identified by RNA-seq data and hub genes identified through the GRNs, we conducted qRT-PCR assay for 13 randomly selected DEGs, including 6 hub genes /TFs in the pink module (Fig. 7A, B, C, D, E, F), and 7 hub genes/TFs in the brown module (Fig. 7G, H, I; Additional file 1: Fig. S14 A, B, C, D). We found that 9 genes (ca. 69%), including 6 (100%) genes in the pink module, 3 (ca. 43%) genes in the brown modules, exhibited a similar trend of expression change under cold treatment as compared to the results from RNA-seq (Fig. 7). For instance, Nitab4.5_0001083g0010 (NtCBF2) was up-regulated at all time points examined (0.5 h, 1 h, 2 h, 4 h, 8 h) as compared to the CK (0 h) (Fig. 7A), Nitab4.5_0001088g0190 was up-regulated in YT relative to TT in each time point (Fig. 7G), while the remaining 4 genes were distinct from the RNA-seq results (Additional file 1: Fig. S14). Inconsistency between RNA-seq and qRT-PCR assay were most likely caused by the variations in different batch of biological samples with cold treatment. Even though we followed the same procedures as we did for RNA-seq for growing the material and cold treatment, it is hard to obtain 100% reproducibility at gene expression levels as before. We also cannot exclude the possibility that parts of DEGs identified by RNA-seq could be false positive, which are possibly not directly caused by cold treatment.

Fig. 7
figure 7

qRT-PCR assay of some key genes identified in the network. A dual y-axis plot illustrating in parallel comparisons of each individual gene examined detected by RNA-seq (expression fold change compared to 0 h of TT, blue) and qRT-PCR (relative expression levels, red). The right coordinate (y axis) with blue represents the expression fold change of RNA-seq, the left coordinate (y axis) with red represents the relative expression levels of qRT-PCR. The transcription levels of each gene in each time point were normalized relative to the internal control (NtActin7). Relative expression levels of genes examined were calculated and expressed as 2-ΔΔCT [63] relative to the expression levels of the corresponding NtActin7, which were set as 1.0. The mean expression levels were calculated from three biological replicates. Error bars are standard deviations of three biological replicates. Significance test was determined using Student’s t test. “ns”: not significant, “*”: p < 0.05, “**”: p < 0.01, “***”: p < 0.001

Full size image

To further confirm the role of NtCBF2 in enhancing cold tolerance, we generated CRISPR-Cas9 based NtCBF2 knockout (KO) transgenic line with K326 background (Fig. 8A, Additional file 1: Fig. S13). After conducting cold treatment, we found that the transgenic line with loss of NtCBF2 gene function exhibited two droop leaves at 5 h post cold treatment, by contrast, the cold damage was even severe for extended cold treatment, since the transgenic line exhibited four droop leaves at 8 h post cold treatment as compared to the WT grown in the normal conditions. This result confirmed that NtCBF2 gene can function in cold responses and tolerance, indicating the accuracy of constructed TF-centered regulatory network related to cold responses in tobacco. Furthermore, we compared the time course related differential expression pattern of qRT-PCR results for NtCBF2 between TT and YT. We found that fold changes of its expression levels between YT and TT were increased at 0.5 h, 4 h and 8 h post cold treatment as compared to CK (0 h), while the cold-induced expression trend of NtCBF2 was opposite between YT and TT at 0.5 h relative to 4 h/8 h (Fig. 8B). After cold treatment, NtCBF2 was more expressed at 0.5 h, but much less expressed at 4 h and 8 h in YT than in TT (Fig. 8B). To assess if there exist any sequence variations resulting in differential expression of NtCBF2 under cold treatment between two varieties, we conducted RT-PCR to obtain full CDS sequences, and conducted PCR to amplify the promoter region (1.5 kb upstream of ATG) from each variety for sequencing. We did not observe any sequence polymorphisms in the cDNA and promoter regions of NtCBF2 between YT and TT (Additional file 1: Fig. S15). The underlying regulatory mechanisms need to be further investigated. The possible reasons related to cold inducible expression change of this gene between YT and TT were discussed in the discussion section.

Fig. 8
figure 8

Experimental validation of NtCBF2. A CRISPR-Cas9 based validation of NtCBF2 gene functioning in cold tolerance or responses. B The time course related differential expression of NtCBF2 between TT and YT

Full size image

Discussion

To cope with a variety of environmental stresses adversely affecting plant normal growth and development, sessile plants have to adjust changes at multiple levels to adapt to these stresses for better survival. They include morphological, physiological, biochemical, metabolic, genetic and epigenetic changes [64,65,66,67]. Global transcriptional changes in response to cold stress have been intensively investigated in various plant species [22, 23, 32, 68,69,70]. However, transcriptomic network studies related to cold treatment are quite limited in tobacco. In this study, we conducted comprehensive multiple time-point transcriptomic analyses between cold-susceptible tobacco cultivar TT and cold-tolerance tobacco cultivar YT. We observed genes up- or down-regulated at each time point in each cultivar or between two cultivars (Additional file 2: Table S2), including cold-responsive genes or TFs such as protein kinases and NtCBF genes (Additional file 2: Table S6). Similar comparative transcriptomic studies related to cold treatment have also been conducted in two tobacco varieties with contrasting cold sensitivity [71], due to tobacco genome sequences unavailable at that time, 3557 cold induced DEGs were identified using the assembled all unique genes and 1840 DEGs were functionally annotated when using Potato (Solanum tuberosum) genome as reference. For a parallel comparison between two studies, we reanalyzed 1840 DEGs as we did by using tobacco (Nicotiana Tabacum) genome sequence as reference for gene annotations. We found that only 940 DEGs best matched in Nicotiana Tabacum genome. According to the overlapping comparison, we found that 354 out of 940 (ca. 38%) DEGs were also present in 9583 DEGs and 3 of them, Nitab. 4.5_0000532G0070, Nitab. 4.5_0001244G0020 and Nitab. 4.5_0002394G0010 (NtCBF5), were present in 100 hub genes in our study. Inconsistency between two cold related studies in tobacco is most likely caused by using different reference genomes for identification of DEGs and subsequent gene annotation.

Moreover, our study also provided evidence showing transcriptional changes between TT and YT, which possibly explains transcriptome-based differential responses to cold stress between TT and YT. It provided evidence showing common mechanisms involved in cold responses between TT and YT. NtCBF genes, cold-responsive genes, were up-regulated in both TT and YT. In particular, TT and YT share similar co-expression gene networks (the pink module, Fig. 5B). We also observed genotype-related gene expression changes under normal growth condition (CK). There were 857 up-regulated genes and 750 down-regulated genes in YT compared to TT before cold treatment (Fig. 2A). Moreover, YT displayed less gene expression changes than TT under cold treatment (Additional file 2: Table S2), thus YT can spend less energy to remodel chromatin structure for reprogramming gene transcription than TT. We assume that those DEGs before cold induction may help to enhance YT to tolerate cold stress at initial stages in response to cold stress. In addition, we found that TT-specific up-regulated genes including gene upregulated in TT but down-regulated in YT exhibited more GO term enrichment associated with various biological and metabolic processes and stress responses. Plants need to consume more energy to complete those processes, thus have to use less energy to cope with cold stress, therefore resulting in sensitivity to cold treatment. Eukaryotic cells need to spend energy from oxidation of metabolic fuels to accomplish all kinds of energy-dependent biological processes [72]. Peroxisome with ATP supply is essential for seedling development in Arabidopsis [73]. By contrast, YT-specific up-regulated genes including gene up-regulated in YT but down-regulated in TT mainly function in photosynthesis, membrane functions, and production of energy. Thus, YT can produce more energy and maintain membrane stability to facilitate the plants to endure cold injury. It has been reported that plants with cold acclimation have increased photosynthetic capacity [74,75,76,77,78]. Remodeling of membrane proteins and lipids and maintenance of membrane fluidity and permeability are one of the principal adaptation strategies for plants accommodating to cold stress [79,80,81].

In addition, we observed differential co-expression network between TT and YT. Genes in the brown module with GO terms enriched in stress responses were up-regulated in YT instead of TT (Fig. 5C, Additional file 1: Fig. S11B). More importantly, we identified 3 TF (ABI3/VP1, ARR-B and WRKY) mainly functioning in differential cold responses between two cultivars, and 3 key TFs (GRAS, AP2-EREBP and C2H2) primarily involved in cold responses (Fig. 6). Some of them have been well documented to be involved in cold or stress responses in other plant species. For example, C2H2 zinc finger gene family (ZAT12 and ZAT10) have been reported to be involved in stress responses [82, 83]. Tobacco GRAS are homologous to SCL13 and PAT1 in Arabidopsis, which are involved in phytochrome A (PhyA) signal transduction and B (PhyB)-related environmental temperature responses, respectively [84,85,86]. Ethylene responsive element binding factors (ERFs) are involved in various abiotic stresses, such as cold, drought and salinity [87]. COBL7 gene has been reported to function in cellulose biosynthesis in Arabidopsis and rice [88, 89]. The involvement of AP2-EREBP families in stress responses have been documented in a couple of plants, such as Arabidopsis [90], rice [91], maize [92], cotton [93], Brachypodium Distachyon [94], Ammopiptanthus nanus [95]. CERK1 has been reported to mediate roles of chitin in salt stress in Arabidopsis [96]. AtWRKY21 and ARR 12 have been reported to act as negative regulators during drought stress in Arabidopsis [52, 97]. Salinity stress can increase expression of PSP in Brassica juncea [51]. Ser/Thr phosphatase type 2A (PP2A) plays crucial roles in the regulation of adaptive biotic and abiotic stress responses in plants [98]. Arabidopsis transgenic plants with overexpression of AB13/VP1 gens (MtRAV3) from M. truncatula increase tolerance to cold treatment [99]. Rice transgenic plants with overexpression of GSA1 encoding a UDP-glucosyltransferase can enhance tolerance to heat, NaCl and PEG stresses [56]. Our study also showed that proteinase kinases were differentially expressed under cold stress, suggesting that protein kinases can potentially act as key regulators in response to cold treatment in tobacco. Protein kinases have been reported to play vital roles in cascading signal transduction during abiotic stresses in plants [100, 101]. Maize plants with constitutive expression of tobacco mitogen-activated protein kinase kinase kinase (Nicotiana PK1) exhibited enhanced freezing resistance [102]. Tobacco plants with overexpression of Malus domestica calcium-dependent protein kinases (MdCPK1a) gene had low levels of ROS accumulation, resulting in increased cold resistance [103]. The roles of protein kinases in cold responses in tobacco need to be further investigated.

In addition, it is worth noting that NtCBF2 was selected through WGCNA using DEGs between TT and YT, but was only functionally validated in K326 background instead of either cultivar examined. Higher expression of NtCBF2 at the early stage of cold treatment may facilitate YT to enhance cold resistance than TT. This was evidenced by the qRT-PCR results that NtCBF2 was more expressed at 0 h (before cold treatment) and 0.5 h post cold treatment but was much less expressed at 4 h and 8 h after cold treatment in YT as compared to TT (Fig. 8B). Thus, YT may activate cold responses earlier, thereby being subjected to less cold damage than TT at the beginning of cold treatment. In addition to NtCBF2, YT had more cold responsive genes with higher expression than TT after cold treatment (Fig. 5). Thus, differential cold responses between two varieties are possibly caused by differential expression of some of key cold responsive genes or potential difference in functions of their downstream targets of cold responsive TFs.

Given that no sequence polymorphisms occurred in the cDNA and proximal promoter regions of NtCBF2 between YT and TT (Additional file 1: Fig. S15), differential transcription of NtCBF2 between YT and TT could be caused by several possibilities: 1) differential expression of its upstream regulators such as TFs; 2) sequence variations occurred in the distal cis-regulatory elements (CREs) affecting the binding of the regulatory protein or protein complex; 3) occurrence of distinct epigenetic mechanisms, such as DNA methylation, histone modifications and chromatin openness. Further evidence needs to be provided. Epigenetic characterization like DNA methylation and histone modifications may help to reveal possible mechanisms responsible for differential expression of cold responsive genes like NtCBF2 under cold treatment between YT and TT.

Conclusions

Collectively, through comprehensive multiple time-point transcriptomic and network analyses, we identified key genes and transcription factors that contribute to differential cold tolerance of two tobacco cultivars. Thus, our study provides valuable transcriptomic resources for studying cold adaptation mechanisms in plants. It helps to identify key cold responsive genes for genetic modifications of tobacco cultivars for improving cold/freezing tolerance.

Materials and methods

Plant growth and cold treatment

Two cultivated tobacco (Nicotiana tabacum L.) varieties with contrasting cold responses, Taiyan8 tobacco (TT, cold susceptibility) and Yanyan97 tobacco (YT, cold resistance) provided by Drs. Risheng Hu and Yangyang Li from Hunan Tobacco Research Institute, were used in this study. The seeds were sterilized in 75% ethanol for 10 min at room temperature (RT) and washed with sterile water for five times. Sterilized seeds were germinated and grown in 1/2 MS medium in an illuminated incubator (16/8 h day/night light cycle, 26/24 °C day/night temperature setting with 70% relative humidity). Tobacco seedlings at the five-leaf stage were used for cold treatment at 4 °C. Leaves were collected from CK and cold-treated plants at different time points (0, 0.5 h, 1.0 h, 2.0 h, 4.0 h and 8.0 h). Collected samples were immediately frozen in liquid nitrogen for RNA preparation. Each sample was biologically triplicated for transcriptomic analyses.

RNA-seq library preparation and data analyses

Total RNA was prepared from tobacco leaves using TRIzol Reagent (ThermoFisher Scientific, Cat # 15596026). After complete removal of genomic DNA contamination using DNase treatment, total 10 μg purified RNA was used for mRNA enrichment using poly(T)-conjugated magnetic beads. mRNA-seq libraries were constructed using an Illumina TruSeq RNA Sample Prep Kit following the manufacturer’s manual. All libraries were sequenced on Illumina Hiseq4500 platform.

The raw sequence data were processed to remove low-quality reads using FastQC (v0.11.5) and Trim_Galore (v 0.6.4). The filtered high-quality reads were mapped to the Nicotiana tabacum genome [104] (ftp://ftp.solgenomics.net/genomes/Nicotiana_tabacum/edwards_et_al_2017/annotation/) using HISAT2 (v2.0.5) with default parameters. Reads with a mapping quality below 20 were excluded for further analyses using SAMtools (v1.9). StringTie (v1.3.4d) and FeatureCounts (v1.6.4) were used to obtain transcripts per million (TPM) and read counts for all annotated genes in the tobacco genome. Correlation between the biological triplicates was determined using Pearson correlation coefficient with TPM value.

A Principal Component Analysis (PCA) was performed using prcomp utilities in the R package. Differential gene expression levels were analyzed using DESeq2. Genes with expression levels of |log2(fold change)| > 1 and padj. < 0.05 were retained for downstream analyses. Gene Ontology (GO) enrichment analyses were conducted by using the online website AgriGO (http://systemsbiology.cau.edu.cn/agriGOv2/specises_analysis.php?&SpeciseID=15& latin=Nicotiana_tabacum), which is specifically designed for agricultural species.

qRT-PCR assay

For qRT-PCR assay, total RNA was extracted from the samples with similar treatment as we did for RNA-seq. The extracted RNA was treated with DNase for removal of genomic DNA contamination followed by reverse-transcription for synthesis of the first cDNA strand. qRT-PCR was conducted using the SYBR mix. NtActin7 was used as an internal control. Relative expression levels of genes examined were calculated from three independent biological replicates and expressed as 2-ΔΔCT [63]. Primers for qRT-PCR assay are listed in Additional file 2: Table S9.

Identification of NtCBF genes

For identification of cold-responsive NtCBF genes, We followed the NtCBF gene named by dehydration-responsive element-binding (DREB) protein in the previous literatures [41]. To extract the corresponding gene id, the DREB protein sequences were further aligned to Nicotiana tabacum protein genome (ftp://ftp.solgenomics.net/genomes/Nicotiana_tabacum/edwards_et_al_2017/annotation/) using Protein-Protein BLAST (v2.9.0+). Detailed information of 21 CBF genes is listed in Additional file 2: Table S6.

Generation of the K326 transgenic line

The Nitab4.5_0001083g0010.1(NtCBF2) gene identified in the GRN was chose for preparing a CRISPR-Cas9 vector (Additional file 2: Table S9). The vector was used for generation of transgenic plant using leaf disc transformation as the published procedures [105]. Homozygous transgenic plants were screening through sequencing and grown in the growth chamber for collecting more seeds for the cold treatment experiments.

Co-expression network analyses

Weighted Gene Co-Expression Network Analysis (WGCNA) package [50] was used in R for co-expression network analyses. The 9583 DEGs with TPM values greater than 1 at least at one time point in two cultivars were used for WGCNA (https://horvath.genetics.ucla.edu/html/CoexpressionNetwork/Rpackages/WGCNA/Tutorials/). A convenient 1-step network construction and module detection function were used for generation of co-expression gene network and identification of each module with distinct functions. The eigenvalues were used to indicate association of genes between each module and at each time point within each cultivar. All co-expressed DEGs were grouped into 17 cultivar- and time-point-specific modules. Each network was visualized using Cytoscape (v.3.7.2).

Genetic regulatory networks (GRNs) analyses

We followed the user manual of the R package GENIE3 (v1.12.0), (https://bioconductor.riken.jp/packages/release/bioc/vignettes/GENIE3/inst/doc/GENIE3.html) [58], to generate a regulatory network containing 9583 DEGs at each time point in the two varieties and 11 identified hub TFs as regulators. Subsequently, we set the threshold (weight > =0.2) to obtain a regulatory network diagram with a strong regulatory relationship between 11 TFs and 3308 DEGs (Additional file 1: Fig. S12) (See Additional file 2: Table S7 for detailed information of the regulation network diagram). The GRNs were visualized using Cytoscape (v.3.7.2).

Statement of plant material used in the study

Study protocol must comply with relevant institutional, national, and international guidelines and legislation.

Availability of data and materials

The raw RNA-seq data used in this study were deposited in the NCBI Gene Expression Omnibus (GEO) (http://www.ncbi.nlm.nih.gov/geo/) with accession number GSE173352. The materials of this study were provided by Hunan Tobacco Research Institute. A voucher specimen of this material has not been deposited in a publicly available herbarium. Correspondence and requests for materials should be addressed to Wenli Zhang (wzhang25@njau.edu.cn) or Risheng Hu (495298768@qq.com).

Abbreviations

ABI3/VP1:

ABSCISIC ACID INSENSITIVE 3/VIVIPAROUS 1, RAV

ARR:

Type A response regulator

CBF:

C-REPEAT BINDING FACTOR

COR:

Cold Regulated

CPL:

CTD phosphatase-like

DEGs:

Differentially expressed genes

DREB:

Dehydration-responsive element-binding

ERFs:

Ethylene responsive element binding factors

GO:

Gene Ontology

GRNs:

Genetic regulatory networks

PCA:

Principal Component Analysis

PKs:

Protein kinases

PP2A:

Ser/Thr phosphatase type 2A

PSP:

3-phosphoserine phosphatase

RT:

Room temperature

TPM:

Transcripts per million

TFs:

Transcription factors

TT:

Tai tobacco

WGCNA:

Weighted Gene Co-Expression Network Analysis

YT:

Yan tobacco

References

  1. Cramer GR, Urano K, Delrot S, Pezzotti M, Shinozaki K. Effects of abiotic stress on plants: a systems biology perspective. BMC Plant Biol. 2011;11:163.

    Article  PubMed  PubMed Central  Google Scholar 

  2. Yadav SK. Cold stress tolerance mechanisms in plants. A review Agronomy for Sustainable Development 2010;30:515–527.

  3. Guo X, Liu D, Chong K. Cold signaling in plants: insights into mechanisms and regulation. J Integr Plant Biol. 2018;60:745–56.

    Article  PubMed  Google Scholar 

  4. Sanghera GS, Wani SH, Hussain W, Singh NB. Engineering cold stress tolerance in crop plants. Curr Genomics. 2011;12:30–43.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  5. Zhu J, Dong CH, Zhu JK. Interplay between cold-responsive gene regulation, metabolism and RNA processing during plant cold acclimation. Curr Opin Plant Biol. 2007;10:290–5.

    Article  CAS  PubMed  Google Scholar 

  6. Theocharis A, Clement C, Barka EA. Physiological and molecular changes in plants grown at low temperatures. Planta. 2012;235:1091–105.

    Article  CAS  PubMed  Google Scholar 

  7. Thomashow MF. PLANT COLD ACCLIMATION: freezing tolerance genes and regulatory mechanisms. Annu Rev Plant Physiol Plant Mol Biol. 1999;50:571–99.

    Article  CAS  PubMed  Google Scholar 

  8. Hoermiller II, Naegele T, Augustin H, Stutz S, Weckwerth W, Heyer AG. Subcellular reprogramming of metabolism during cold acclimation in Arabidopsis thaliana. Plant Cell Environ. 2017;40:602–10.

    Article  CAS  PubMed  Google Scholar 

  9. Chinnusamy V, Zhu JK, Sunkar R. Gene regulation during cold stress acclimation in plants. Methods Mol Biol. 2010;639:39–55.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  10. Zhu JK. Abiotic stress signaling and responses in plants. Cell. 2016;167:313–24.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  11. Knight MR, Knight H. Low-temperature perception leading to gene expression and cold tolerance in higher plants. New Phytol. 2012;195:737–51.

    Article  CAS  PubMed  Google Scholar 

  12. Ritonga FN, Chen S. Physiological and molecular mechanism involved in cold stress tolerance in plants. Plants (Basel). 2020;9:560.

    Article  CAS  PubMed Central  Google Scholar 

  13. Stockinger EJ, Gilmour SJ, Thomashow MF. Arabidopsis thaliana CBF1 encodes an AP2 domain-containing transcriptional activator that binds to the C-repeat/DRE, a cis-acting DNA regulatory element that stimulates transcription in response to low temperature and water deficit. Proc Natl Acad Sci U S A. 1997;94:1035–40.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  14. Gilmour SJ, Zarka DG, Stockinger EJ, Salazar MP, Houghton JM, Thomashow MF. Low temperature regulation of the Arabidopsis CBF family of AP2 transcriptional activators as an early step in cold-induced COR gene expression. Plant J. 1998;16:433–42.

    Article  CAS  PubMed  Google Scholar 

  15. Chinnusamy V, Zhu J, Zhu JK. Cold stress regulation of gene expression in plants. Trends Plant Sci. 2007;12:444–51.

    Article  CAS  PubMed  Google Scholar 

  16. Calixto CPG, Guo W, James AB, Tzioutziou NA, Entizne JC, Panter PE, et al. Rapid and dynamic alternative splicing impacts the Arabidopsis cold response transcriptome. Plant Cell. 2018;30:1424–44.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  17. Penfield S. Temperature perception and signal transduction in plants. New Phytol. 2008;179:615–28.

    Article  CAS  PubMed  Google Scholar 

  18. Wang DZ, Jin YN, Ding XH, Wang WJ, Zhai SS, Bai LP, et al. Gene regulation and signal transduction in the ICE-CBF-COR signaling pathway during cold stress in plants. Biochemistry (Mosc). 2017;82:1103–17.

    Article  CAS  Google Scholar 

  19. Barah P, Jayavelu ND, Rasmussen S, Nielsen HB, Mundy J, Bones AM. Genome-scale cold stress response regulatory networks in ten Arabidopsis thaliana ecotypes. BMC Genomics. 2013;14:722.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  20. Li Q, Byrns B, Badawi MA, Diallo AB, Danyluk J, Sarhan F, et al. Transcriptomic insights into Phenological development and cold tolerance of wheat grown in the field. Plant Physiol. 2018;176:2376–94.

    Article  CAS  PubMed  Google Scholar 

  21. Svensson JT, Crosatti C, Campoli C, Bassi R, Stanca AM, Close TJ, et al. Transcriptome analysis of cold acclimation in barley Albina and xantha mutants. Plant Physiol. 2006;141:257–70.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  22. Guan S, Xu Q, Ma D, Zhang W, Xu Z, Zhao M, et al. Transcriptomics profiling in response to cold stress in cultivated rice and weedy rice. Gene. 2019;685:96–105.

    Article  CAS  PubMed  Google Scholar 

  23. Li M, Sui N, Lin L, Yang Z, Zhang Y. Transcriptomic profiling revealed genes involved in response to cold stress in maize. Funct Plant Biol. 2019;46:830–44.

    Article  CAS  PubMed  Google Scholar 

  24. Calzadilla PI, Maiale SJ, Ruiz OA, Escaray FJ. Transcriptome response mediated by cold stress in Lotus japonicus. Front Plant Sci. 2016;7:374.

    Article  PubMed  PubMed Central  Google Scholar 

  25. Zhang Y, Zhang Y, Lin Y, Luo Y, Wang X, Chen Q, et al. A transcriptomic analysis reveals diverse regulatory networks that respond to cold stress in strawberry (Fragariaxananassa). Int J Genomics. 2019;2019:7106092.

    PubMed  PubMed Central  Google Scholar 

  26. Kang WH, Sim YM, Koo N, Nam JY, Lee J, Kim N, et al. Transcriptome profiling of abiotic responses to heat, cold, salt, and osmotic stress of Capsicum annuum L. Sci Data. 2020;7:17.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  27. Yang QS, Gao J, He WD, Dou TX, Ding LJ, Wu JH, et al. Comparative transcriptomics analysis reveals difference of key gene expression between banana and plantain in response to cold stress. BMC Genomics. 2015;16:446.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  28. Chen H, Chen X, Chen D, Li J, Zhang Y, Wang A. A comparison of the low temperature transcriptomes of two tomato genotypes that differ in freezing tolerance: Solanum lycopersicum and Solanum habrochaites. BMC Plant Biol. 2015;15:132.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  29. Li Y, Wang X, Ban Q, Zhu X, Jiang C, Wei C, et al. Comparative transcriptomic analysis reveals gene expression associated with cold adaptation in the tea plant Camellia sinensis. BMC Genomics. 2019;20:624.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  30. Zhang X, Teixeira da Silva JA, Niu M, Li M, He C, Zhao J, et al. Physiological and transcriptomic analyses reveal a response mechanism to cold stress in Santalum album L. leaves. Sci Rep. 2017;7:42165.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  31. Sun LL, Wang YB, Wang RL, Wang RT, Zhang P, Ju Q, et al. Physiological, transcriptomic, and metabolomic analyses reveal zinc oxide nanoparticles modulate plant growth in tomato. Environmental Science-Nano. 2020;7:3587–604.

    Article  CAS  Google Scholar 

  32. Jin J, Zhang H, Zhang J, Liu P, Chen X, Li Z, et al. Integrated transcriptomics and metabolomics analysis to characterize cold stress responses in Nicotiana tabacum. BMC Genomics. 2017;18:496.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  33. Goossens A, Hakkinen ST, Laakso I, Seppanen-Laakso T, Biondi S, De Sutter V, et al. A functional genomics approach toward the understanding of secondary metabolism in plant cells. Proc Natl Acad Sci U S A. 2003;100:8595–600.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  34. Schweiger R, Schwenkert S. Protein-protein interactions visualized by bimolecular fluorescence complementation in tobacco protoplasts and leaves. J Vis Exp. 2014;9:51327.

    Google Scholar 

  35. Zhao L, Liu F, Xu W, Di C, Zhou S, Xue Y, et al. Increased expression of OsSPX1 enhances cold/subfreezing tolerance in tobacco and Arabidopsis thaliana. Plant Biotechnol J. 2009;7:550–61.

    Article  CAS  PubMed  Google Scholar 

  36. Khodakovskaya M, McAvoy R, Peters J, Wu H, Li Y. Enhanced cold tolerance in transgenic tobacco expressing a chloroplast omega-3 fatty acid desaturase gene under the control of a cold-inducible promoter. Planta. 2006;223:1090–100.

    Article  CAS  PubMed  Google Scholar 

  37. Kodama H, Hamada T, Horiguchi G, Nishimura M, Iba K. Genetic enhancement of cold tolerance by expression of a gene for chloroplast [omega]-3 fatty acid desaturase in transgenic tobacco. Plant Physiol. 1994;105:601–5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  38. Zhuo C, Wang T, Guo Z, Lu S. Overexpression of MfPIP2-7 from Medicago falcata promotes cold tolerance and growth under NO3 (−) deficiency in transgenic tobacco plants. BMC Plant Biol. 2016;16:138.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  39. Wei Y, Chen H, Wang L, Zhao Q, Wang D, Zhang T. Cold acclimation alleviates cold stress-induced PSII inhibition and oxidative damage in tobacco leaves. Plant Signal Behav. 2022;17:2013638.

    Article  PubMed  CAS  Google Scholar 

  40. Ya J, Zhang C, Yang H, Yang Y, Huang C, Tian Y. Lu X. proteomic analysis of cold stress responses in tobacco seedlings. Afr J Biotechnol. 2011;10:18991–9004.

    Google Scholar 

  41. Zhao Q, Xiang X, Liu D, Yang A, Wang Y. Tobacco transcription factor NtbHLH123 confers tolerance to cold stress by regulating the NtCBF pathway and reactive oxygen species homeostasis. Front Plant Sci. 2018;9:381.

    Article  PubMed  PubMed Central  Google Scholar 

  42. Wang C, Deng P, Chen L, Wang X, Ma H, Hu W, et al. A wheat WRKY transcription factor TaWRKY10 confers tolerance to multiple abiotic stresses in transgenic tobacco. PLoS One. 2013;8:e65120.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  43. Lin P, Shen C, Chen H, Yao XH, Lin J. Improving tobacco freezing tolerance by co-transfer of stress-inducible CbCBF and CbICE53 genes. Biol Plant. 2017;61:520–8.

    Article  CAS  Google Scholar 

  44. Khan SA, Li MZ, Wang SM, Yin HJ. Revisiting the role of plant transcription factors in the Battle against abiotic stress. Int J Mol Sci. 2018;19:1634.

    Article  PubMed Central  CAS  Google Scholar 

  45. Singh K, Foley RC, Onate-Sanchez L. Transcription factors in plant defense and stress responses. Curr Opin Plant Biol. 2002;5:430–6.

    Article  CAS  PubMed  Google Scholar 

  46. Kitsios G, Doonan JH. Cyclin dependent protein kinases and stress responses in plants. Plant Signal Behav. 2011;6:204–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  47. Jonak C, Kiegerl S, Ligterink W, Barker PJ, Huskisson NS, Hirt H. Stress signaling in plants: a mitogen-activated protein kinase pathway is activated by cold and drought. Proc Natl Acad Sci U S A. 1996;93:11274–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  48. Mizoguchi T, Ichimura K, Shinozaki K. Environmental stress response in plants: the role of mitogen-activated protein kinases. Trends Biotechnol. 1997;15:15–9.

    Article  CAS  PubMed  Google Scholar 

  49. Yang L, Wu K, Gao P, Liu X, Li G, Wu Z. GsLRPK, a novel cold-activated leucine-rich repeat receptor-like protein kinase from Glycine soja, is a positive regulator to cold stress tolerance. Plant Sci. 2014;215-216:19–28.

    Article  CAS  PubMed  Google Scholar 

  50. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 2008;9:559.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  51. Purty RS, Sachar M, Chatterjee S. Structural and expression analysis of salinity stress responsive phosphoserine phosphatase from Brassica juncea (L.). J Proteomics Bioinformatics. 2017;10:119–27.

    Article  Google Scholar 

  52. Zhao KX, Chu SS, Zhang XD, Wang L, Rono JK, Yang ZM. AtWRKY21 negatively regulates tolerance to osmotic stress in Arabidopsis. Environ Exp Bot. 2020;169:103920.

    Article  CAS  Google Scholar 

  53. Mathe C, Garda T, Freytag C, M MH. The role of serine-threonine protein phosphatase PP2A in plant oxidative stress signaling-facts and hypotheses. Int J Mol Sci 2019;20:3028.

  54. Javed T, Shabbir R, Ali A, Afzal I, Zaheer U, Gao SJ. Transcription factors in plant stress responses: challenges and potential for sugarcane improvement. Plants (Basel). 2020;9:491.

    Article  CAS  Google Scholar 

  55. Xie Z, Nolan TM, Jiang H, Yin Y. AP2/ERF transcription factor regulatory networks in hormone and abiotic stress responses in Arabidopsis. Front Plant Sci. 2019;10:228.

    Article  PubMed  PubMed Central  Google Scholar 

  56. Dong NQ, Sun Y, Guo T, Shi CL, Zhang YM, Kan Y, et al. UDP-glucosyltransferase regulates grain size and abiotic stress tolerance associated with metabolic flux redirection in rice. Nat Commun. 2020;11:2629.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  57. Bang W, Kim S, Ueda A, Vikram M, Yun D, Bressan RA, et al. Arabidopsis carboxyl-terminal domain phosphatase-like isoforms share common catalytic and interaction domains but have distinct in planta functions. Plant Physiol. 2006;142:586–94.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  58. Huynh-Thu VA, Irrthum A, Wehenkel L, Geurts P. Inferring regulatory networks from expression data using tree-based methods. PLoS One. 2010;5:e12776.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  59. Kim CY, Vo KTX, Nguyen CD, Jeong DH, Lee SK, Kumar M, et al. Functional analysis of a cold-responsive rice WRKY gene, OsWRKY71. Plant Biotechnol Rep. 2016;10:13–23.

    Article  CAS  Google Scholar 

  60. Zou C, Jiang W, Yu D. Male gametophyte-specific WRKY34 transcription factor mediates cold sensitivity of mature pollen in Arabidopsis. J Exp Bot. 2010;61:3901–14.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  61. Yang Y, Liu J, Zhou X, Liu S, Zhuang Y. Identification of WRKY gene family and characterization of cold stress-responsive WRKY genes in eggplant. Peer J. 2020;8:e8777.

    Article  PubMed  PubMed Central  Google Scholar 

  62. Yuan Y, Fang L, Karungo SK, Zhang L, Gao Y, Li S, et al. Overexpression of VaPAT1, a GRAS transcription factor from Vitis amurensis, confers abiotic stress tolerance in Arabidopsis. Plant Cell Rep. 2016;35:655–66.

    Article  CAS  PubMed  Google Scholar 

  63. Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2(−Delta Delta C(T)) method. Methods. 2001;25:402–8.

    Article  CAS  PubMed  Google Scholar 

  64. Lamalakshmi Devi E, Kumar S, Singh B, Sharma S, Beemrote A, Devi C, Chongtham S, Chongtham H, Yumlembam R, Athokpam H, et al. Adaptation Strategies and Defence Mechanisms of Plants During Environmental Stress. In: Medicinal Plants and Environmental Challenges. Edited by Ghorbanpour M, Varma A. Cham: Springer International Publishing; 2017:359–413.

  65. Chang YN, Zhu C, Jiang J, Zhang H, Zhu JK, Duan CG. Epigenetic regulation in plant abiotic stress responses. J Integr Plant Biol. 2020;62:563–80.

    Article  CAS  PubMed  Google Scholar 

  66. Krasensky J, Jonak C. Drought, salt, and temperature stress-induced metabolic rearrangements and regulatory networks. J Exp Bot. 2012;63:1593–608.

    Article  CAS  PubMed  Google Scholar 

  67. Haak DC, Fukao T, Grene R, Hua Z, Ivanov R, Perrella G, et al. Multilevel regulation of abiotic stress responses in plants. Front Plant Sci. 2017;8:1564.

    Article  PubMed  PubMed Central  Google Scholar 

  68. Bashir K, Matsui A, Rasheed S, Seki M. Recent advances in the characterization of plant transcriptomes in response to drought, salinity, heat, and cold stress. F1000Res. 2019;8:F1000.

    Article  PubMed  PubMed Central  Google Scholar 

  69. Winfield MO, Lu C, Wilson ID, Coghill JA, Edwards KJ. Plant responses to cold: transcriptome analysis of wheat. Plant Biotechnol J. 2010;8:749–71.

    Article  CAS  PubMed  Google Scholar 

  70. Li ZB, Zeng XY, Xu JW, Zhao RH, Wei YN. Transcriptomic profiling of cotton Gossypium hirsutum challenged with low-temperature gradients stress. Sci Data. 2019;6:197.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  71. Hu R, Zhu X, Xiang S, Zhan Y, Zhu M, Yin H, et al. Comparative transcriptome analysis revealed the genotype specific cold response mechanism in tobacco. Biochem Biophys Res Commun. 2016;469:535–41.

    Article  CAS  PubMed  Google Scholar 

  72. Bonora M, Patergnani S, Rimessi A, De Marchi E, Suski JM, Bononi A, et al. ATP synthesis and storage. Purinergic Signal. 2012;8:343–57.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  73. Linka N, Theodoulou FL, Haslam RP, Linka M, Napier JA, Neuhaus HE, et al. Peroxisomal ATP import is essential for seedling development in Arabidopsis thaliana. Plant Cell. 2008;20:3241–57.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  74. Rapacz M, Wolanin B, Hura K, Tyrka M. The effects of cold acclimation on photosynthetic apparatus and the expression of COR14b in four genotypes of barley (Hordeum vulgare) contrasting in their tolerance to freezing and high-light treatment in cold conditions. Ann Bot. 2008;101:689–99.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  75. Savitch LV, Barker-Astrom J, Ivanov AG, Hurry V, Oquist G, Huner NP, et al. Cold acclimation of Arabidopsis thaliana results in incomplete recovery of photosynthetic capacity, associated with an increased reduction of the chloroplast stroma. Planta. 2001;214:295–303.

    Article  CAS  PubMed  Google Scholar 

  76. Huner NP, Oquist G, Hurry VM, Krol M, Falk S, Griffith M. Photosynthesis, photoinhibition and low temperature acclimation in cold tolerant plants. Photosynth Res. 1993;37:19–39.

    Article  CAS  PubMed  Google Scholar 

  77. Stitt M, Hurry V. A plant for all seasons: alterations in photosynthetic carbon metabolism during cold acclimation in Arabidopsis. Curr Opin Plant Biol. 2002;5:199–206.

    Article  CAS  PubMed  Google Scholar 

  78. Savitch L, Leonardos E, Król M, Jansson S, Grodzinski B, Huner N, et al. Two different strategies for light utilization in photosynthesis in relation to growth and cold acclimation. Plant Cell Environment. 2002;25:761–71.

    Article  CAS  Google Scholar 

  79. Takahashi D, Li B, Nakayama T, Kawamura Y, Uemura M. Plant plasma membrane proteomics for improving cold tolerance. Front Plant Sci. 2013;4:90.

    PubMed  PubMed Central  Google Scholar 

  80. Zheng G, Tian B, Zhang F, Tao F, Li W. Plant adaptation to frequent alterations between high and low temperatures: remodelling of membrane lipids and maintenance of unsaturation levels. Plant Cell Environ. 2011;34:1431–42.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  81. Zhang H, Dong J, Zhao X, Zhang Y, Ren J, Xing L, et al. Research Progress in membrane lipid metabolism and molecular mechanism in Peanut cold tolerance. Front Plant Sci. 2019;10:838.

    Article  PubMed  PubMed Central  Google Scholar 

  82. Mittler R, Kim Y, Song L, Coutu J, Coutu A, Ciftci-Yilmaz S, et al. Gain- and loss-of-function mutations in Zat10 enhance the tolerance of plants to abiotic stress. FEBS Lett. 2006;580:6537–42.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  83. Davletova S, Schlauch K, Coutu J, Mittler R. The zinc-finger protein Zat12 plays a central role in reactive oxygen and abiotic stress signaling in Arabidopsis. Plant Physiol. 2005;139:847–56.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  84. Torres-Galea P, Hirtreiter B, Bolle C. Two GRAS proteins, SCARECROW-LIKE21 and PHYTOCHROME a SIGNAL TRANSDUCTION1, function cooperatively in phytochrome a signal TRANSDUCTION. Plant Physiol. 2013;161:291–304.

    Article  CAS  PubMed  Google Scholar 

  85. Torres-Galea P, Huang LF, Chua NH, Bolle C. The GRAS protein SCL13 is a positive regulator of phytochrome-dependent red light signaling, but can also modulate phytochrome a responses. Mol Gen Genomics. 2006;276:13–30.

    Article  CAS  Google Scholar 

  86. Jung JH, Domijan M, Klose C, Biswas S, Ezer D, Gao M, et al. Phytochromes function as thermosensors in Arabidopsis. Science. 2016;354:886–9.

    Article  CAS  PubMed  Google Scholar 

  87. Fujimoto SY, Ohta M, Usui A, Shinshi H, Ohme-Takagi M. Arabidopsis ethylene-responsive element binding factors act as transcriptional activators or repressors of GCC box-mediated gene expression. Plant Cell. 2000;12:393–404.

    CAS  PubMed  PubMed Central  Google Scholar 

  88. Sorek N, Sorek H, Kijac A, Szemenyei HJ, Bauer S, Hematy K, et al. The Arabidopsis COBRA protein facilitates cellulose crystallization at the plasma membrane. J Biol Chem. 2014;289:34911–20.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  89. Li Y, Qian Q, Zhou Y, Yan M, Sun L, Zhang M, et al. BRITTLE CULM1, which encodes a COBRA-like protein, affects the mechanical properties of rice plants. Plant Cell. 2003;15:2020–31.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  90. Song CP, Agarwal M, Ohta M, Guo Y, Halfter U, Wang P, et al. Role of an Arabidopsis AP2/EREBP-type transcriptional repressor in abscisic acid and drought stress responses. Plant Cell. 2005;17:2384–96.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  91. Sharoni AM, Nuruzzaman M, Satoh K, Shimizu T, Kondoh H, Sasaya T, et al. Gene structures, classification and expression models of the AP2/EREBP transcription factor family in rice. Plant Cell Physiol. 2011;52:344–60.

    Article  CAS  PubMed  Google Scholar 

  92. Kizis D, Pages M. Maize DRE-binding proteins DBF1 and DBF2 are involved in rab17 regulation through the drought-responsive element in an ABA-dependent pathway. Plant J. 2002;30:679–89.

    Article  CAS  PubMed  Google Scholar 

  93. Liu C, Zhang T. Expansion and stress responses of the AP2/EREBP superfamily in cotton. BMC Genomics. 2017;18:118.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  94. Chen L, Han J, Deng X, Tan S, Li L, Li L, et al. Expansion and stress responses of AP2/EREBP superfamily in Brachypodium distachyon. Sci Rep. 2016;6:21623.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  95. Cao S, Wang Y, Li X, Gao F, Feng J, Zhou Y. Characterization of the AP2/ERF transcription factor family and expression profiling of DREB subfamily under cold and osmotic stresses in Ammopiptanthus nanus. Plants (Basel). 2020;9:455.

    Article  CAS  PubMed Central  Google Scholar 

  96. Espinoza C, Liang Y, Stacey G. Chitin receptor CERK1 links salt stress and chitin-triggered innate immunity in Arabidopsis. Plant J. 2017;89:984–95.

    Article  CAS  PubMed  Google Scholar 

  97. Nguyen KH, Ha CV, Nishiyama R, Watanabe Y, Leyva-Gonzalez MA, Fujita Y, et al. Arabidopsis type B cytokinin response regulators ARR1, ARR10, and ARR12 negatively regulate plant responses to drought. Proc Natl Acad Sci U S A. 2016;113:3090–5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  98. Pais SM, Tellez-Inon MT, Capiati DA. Serine/threonine protein phosphatases type 2A and their roles in stress signaling. Plant Signal Behav. 2009;4:1013–5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  99. Wang S, Guo T, Wang Z, Kang J, Yang Q, Shen Y, et al. Expression of three related to ABI3/VP1 genes in Medicago truncatula caused increased stress resistance and branch increase in Arabidopsis thaliana. Front Plant Sci. 2020;11:611.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  100. Sinha AK, Jaggi M, Raghuram B, Tuteja N. Mitogen-activated protein kinase signaling in plants under abiotic stress. Plant Signal Behav. 2011;6:196–203.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  101. Lin L, Wu J, Jiang M, Wang Y. Plant mitogen-activated protein kinase cascades in environmental stresses. Int J Mol Sci. 2021;22:1543.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  102. Shou H, Bordallo P, Fan JB, Yeakley JM, Bibikova M, Sheen J, et al. Expression of an active tobacco mitogen-activated protein kinase kinase kinase enhances freezing tolerance in transgenic maize. Proc Natl Acad Sci U S A. 2004;101:3298–303.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  103. Dong H, Wu C, Luo C, Wei M, Qu S, Wang S. Overexpression of MdCPK1a gene, a calcium dependent protein kinase in apple, increase tobacco cold tolerance via scavenging ROS accumulation. PLoS One. 2020;15:e0242139.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  104. Edwards KD, Fernandez-Pozo N, Drake-Stowe K, Humphry M, Evans AD, Bombarely A, et al. A reference genome for Nicotiana tabacum enables map-based cloning of homeologous loci implicated in nitrogen utilization efficiency. BMC Genomics. 2017;18:448.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  105. Gao J, Wang G, Ma S, Xie X, Wu X, Zhang X, et al. CRISPR/Cas9-mediated targeted mutagenesis in Nicotiana tabacum. Plant Mol Biol. 2015;87:99–110.

    Article  CAS  PubMed  Google Scholar 

Download references

Acknowledgements

We thank the Bioinformatic Center in Nanjing Agricultural University for providing facilities to assist sequencing data analysis.

Funding

This work was supported by grants from the National Natural Science Foundation of China (32070561); and from the Key Science and Technology Program of Hunan Provincial Tobacco Corporation (17-18Aa05, HN2021KJ04 and HN2021KJ02).

Author information

Author notes

  1. Zhenyu Luo, Zhicheng Zhou and Yangyang Li contributed equally to this work.

Authors and Affiliations

  1. State Key Laboratory for Crop Genetics and Germplasm Enhancement, CIC-MCP, Nanjing Agricultural University, No.1 Weigang, Nanjing, 210095, Jiangsu, China

    Zhenyu Luo, Shentong Tao, Xuejiao Cheng & Wenli Zhang

  2. Hunan Tobacco Research Institute, Changsha, 410128, Hunan, China

    Zhicheng Zhou, Yangyang Li, Zheng-Rong Hu, Jia-Shuo Yang & Risheng Hu

Authors
  1. Zhenyu Luo
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Zhicheng Zhou
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Yangyang Li
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Shentong Tao
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Zheng-Rong Hu
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Jia-Shuo Yang
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Xuejiao Cheng
    View author publications

    You can also search for this author in PubMed Google Scholar

  8. Risheng Hu
    View author publications

    You can also search for this author in PubMed Google Scholar

  9. Wenli Zhang
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

WLZ and RSH conceived the study. WLZ and XJC instructed the study. ZYL analyzed the data. ZCZ, YYL, ZRH and JSY performed the experiments. STT assisted with data analysis. Z.Y.L., Z.C.Z. and Y.Y.L. WLZ wrote the manuscript with contribution from all authors. All authors read and approved the manuscript.

Corresponding authors

Correspondence to Xuejiao Cheng, Risheng Hu or Wenli Zhang.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Luo, Z., Zhou, Z., Li, Y. et al. Transcriptome-based gene regulatory network analyses of differential cold tolerance of two tobacco cultivars. BMC Plant Biol 22, 369 (2022). https://doi.org/10.1186/s12870-022-03767-7

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12870-022-03767-7

Keywords

BMC Plant Biology

ISSN: 1471-2229

Contact us