To gain the profiles of gene expression in Thellungiella under cold condition, eight cDNA libraries were constructed using roots and leaves under normal (control) and low temperature (cold) for 24 h, respectively. The cDNA libraries were sequenced by Illumina Hiseq2000 platform using the paired end method. After removing low quality, N-containing and adaptor-contaminated reads, a total of 96,305,447 clean reads were generated, with an average of ∼12 million reads per libraries (Table 1). Approximately 90% reads from leaves and 85% reads from roots were mapped to the Thellungiella reference genome, and about 7% reads from leaves and 3.5% reads from roots were mapped to multiple regions, respectively (Table 1). The RNA-seq raw sequencing data from this study have been submitted in the SRA database¹ under BioProjects: PRJNA377594 (SRA: SRP101369).

Among the 29,284 genes deposited in Thellungiella genome database, 22,414 (76.5%) genes were detected in the control and cold-treated libraries. A total of 2,782 and 1,430 differentially expressed genes (DEGs) were identified from leaves and roots, respectively (Figure 1 and Supplementary Tables S1, S2). Under cold treatment, 579 and 1,691 genes were up-regulated, 851 and 1,091 genes were down-regulated in roots and leaves, respectively (Figure 1). Our results demonstrated that the expression patterns of the majority of DEGs were different in roots and leaves. For example, among the 1,691 up-regulated genes in leaves, only 269 (15.9%) were also induced in roots (Figure 2). Interestingly, some genes showed the opposite expression trend in roots and leaves upon cold treatment. For example, fifteen DEGs were up-regulated in root, but down-regulated in leaves. There were 68 DEGs down-regulated in root, but up-regulated in leaves (Figure 2). Moreover, 463 DEGs were with the same expression trend in roots and leaves, including 269 up-regulated and 194 down-regulated DEGs (Figure 2 and Supplementary Table S1).

The DEGs were characterized by the assignment of gene ontology (GO) terms using Blast2GO (Conesa et al., 2005) program. GO enrichment analysis showed 241 and 374 GO terms were significant enriched in root and leaf, respectively (Supplementary Table S3). Then, we employed WEGO web-based tool to visualize the biological process, molecular function and cellular component main categories (Figure 3). In the cellular component category, the terms of “cell”, “cell part”, “organelle”, and “membrane” were enriched, implying the potential contribution of cell and cell structure in the process of Thellungiella response to cold condition. For the category of molecular function, “binding”, “catalytic activity” and “nucleic acid binding transcription factor activity” were the top terms. The most abundant terms of biological process were “cellular process”, “metabolic process”, “response to stimulus” and “single-organism process”, suggesting a high degree of metabolic activity changes upon cold treatment (Figure 3).

In order to obtain more biological information for understanding the molecular mechanism and regulatory network of Thellungiella cold tolerance, KEGG enrichment pathway analysis was performed. By applying a cut-off criterion of Q-value < 0.05 and P-value < 1E^-05, the result showed that ten and fifteen pathways were significantly enriched from roots and leaves, respectively (Table 2). Previous studies demonstrated that many metabolic changes for enhancing freeze tolerance in Arabidopsis, such as increasing accumulation of soluble sugars and other compatible osmolytes (Wanner and Junttila, 1999). We found six pathways, “biosynthesis of secondary metabolites”, “metabolic pathways”, “nitrogen metabolism”, “tryptophan metabolism”, “cysteine and methionine metabolism”, and “sulfur metabolism” were enriched both in roots and leaves in response to cold. Interestingly, all these pathways were involved in a particular metabolic process, suggesting that the metabolic process was activated via cold treatment (Table 2). Proline is one of the most effective organic osmolytes in plants, and there is a positive correlation between proline accumulation and plant stress tolerance. A number of studies showed that proline played beneficial roles in plants when exposed to cold condition (Hayat et al., 2012). Previous studies demonstrated that Thellungiella contained higher levels of proline than Arabidopsis under non-stressed condition (Kant et al., 2006). Twenty-six genes involved in “Arginine and proline metabolism” pathway were found to be regulated by cold in our study. Most of these DEGs (18 of 26) were up-regulated, suggesting that these cold induced genes might promote the accumulation of proline in Thellungiella to enhance cold tolerance.

The “biosynthesis of unsaturated fatty acids” pathway was also enriched in Thellungiella leaves. A total of 22 genes involved in this pathway were affected by low temperature (Table 2). Previous study suggested that membrane lipid composition, especially the concentration of unsaturated fatty acid, is highly correlated with plant freezing tolerance (Thomas et al., 2012). The expression regulation of genes for unsaturated fatty acid synthesis might be a key factor contributing to cold tolerance in Thellungiella.

We observed that the pathways related to biosynthesis of phenylpropanoid and flavonoid were enriched in roots of Thellungiella. In addition, “Photosynthesis - antenna proteins” and “Photosynthesis” were enriched. Most of the DEGs in “Photosynthesis-antenna proteins” and “Photosynthesis” pathways were down-regulated in leaves, suggesting the adverse effect of low temperature on photosynthetic carbohydrate metabolism and photochemical reaction (Supplementary Table S2).

In Arabidopsis, environmental temperature affected the expression of clock component related genes, such as timing of cab expression 1 (TOC1), GIGANTEA (GI), circadian clock associated 1 (CCA1), and late elongated hypocotyl (LHY). Here, our results demonstrated that the pathway of “circadian rhythm - plant” was enriched in Thellungiella after cold treatment. The expression of many genes involved in circadian rhythm including LHY, CCA1, TOC1 and GI were also induced or inhibited in Thellungiella upon cold treatment (Table 2 and Supplementary Table S2).

To validate the RNA-seq data, quantitative real-time PCR (qRT-PCR) was performed for 25 genes with different expression levels and functional assignments (Supplementary Table S4). Among them, seven genes were significantly induced upon cold treatment both in roots and leaves. Two of these genes encoded salt and low temperature response protein, one encoded dehydrin. A putative phosphatidylethanolamine-binding protein gene, and three genes with unknown function were also selected for qRT-PCR validation. Other ten selected genes were down-regulated both in roots and leaves. These genes included one aquaporin TIP2-1, one NAC domain protein, one stress-induced sti1-like protein coding genes and several genes with unknown function. In addition, eight genes that down-regulated in roots but up-regulated in leaves were also selected for qRT-PCR analysis. These genes included O-methyltransferase family protein, peroxidase, cytochrome P450, leucine-rich repeat receptor-like protein kinase, MYB and zine finger AN1 domain-contain protein coding genes (Supplementary Table S4). The qRT-PCR results of these 25 genes in leaves were in agreement with the RNA-seq data. The relative expression level (log₂ cold/control) estimated by RNA-seq and qRT-PCR were strongly correlated (R²=0.9695) in leaves (Supplementary Table S4 and Figures 4, 5). In roots, the majority of genes (23 of 25) showed similar expression patterns except Thhalv10004977m and Thhalv10002333m, which were down-regulated in RNA-seq, but up-regulated in qRT-PCR. In roots, the Pearson’s coefficient was 0.7664 which was lower than that in leaves (Figure 4 and Supplementary Figure S2).

Studies showed that transcription factors played important roles in plant response to low temperature and other adverse stresses (Singh et al., 2002; Chen and Zhu, 2004; Agarwal P.K. et al., 2006). We observed that many transcript factors were response to cold in Thellungiella. Among the 1,430 DEGs in roots, 159 genes (11.1%) were annotated as different types of transcription factors (Table 3). About half of these transcription factor genes (74 of 159) were up-regulated, the rest 85 genes were down-regulated. According to functional annotation, these transcription factors were classified into 18 categories, such as abscisic acid responsive, NAC, Zinc finger domain, AP2, MYB, bHLH and WRKY etc. The transcription factors Zinc finger domains (27 genes, 17.0%), ethylene responsive (11 genes, 6.9%) and calcium ion binding (10 genes, 6.3%) were the three major families of the cold-regulated transcription factors in roots.

Among the 2,782 DEGs in leaves, 232 (8.3%) were transcription factors, including 131 up-regulated and 101 down-regulated genes. These transcription factors were classified into 18 categories (Table 3). Among these differential expressed transcription factors in leaves, Zinc finger domain, WRKY and MYB are enriched. Heat shock proteins (HSPs) were implicated in plant heat stress tolerance (Nover et al., 2002). Here, we found five HSPs were up-regulated and four HSPs were down-regulated upon cold treatment. WRKY transcription factors played important roles in plant responses to biotic and abiotic stress (Eulgem et al., 2000). Twelve differentially expressed WRKY transcription factor genes were identified in leaves, and eleven of them were up-regulated. NAC transcription factors were involved in many aspects of plant growth and development, and response to abiotic stress (Nuruzzaman et al., 2013). In our RNA-seq data, nine and three differentially expressed NAC transcription factor genes were identified in leaves and roots, respectively. These results indicated that HSF, WRKY and NAC transcription factors were involved in plant responses to various stresses, and suggested that cold stress might share common molecular mechanism with other abiotic stresses. In addition, we noticed that the expression of ten calcium ion binding transcription factors was all suppressed in roots upon cold treatment, while most of them were induced in leaves, implying that these genes might be functioning in different ways in roots and leaves (Table 3 and Supplementary Table S1). These results suggested the existence of differences in cold responsive regulatory networks between roots and leaves in Thellungiella.

Previous studies showed that some MYB transcription factors were involved in cold stress tolerance, such as AtMYB15 in Arabidopsis (Agarwal M. et al., 2006), MYB55 and OsMYBS3 in rice (Su et al., 2010; Elkereamy et al., 2012). In our data, the expression of six and ten MYB genes was altered after cold treatment in roots and leaves, respectively (Table 3 and Supplementary Tables S1, S2). Interestingly, Thhalv10008152m encoding a MYB transcription factor displayed opposite expression trend in roots and leaves. In roots, the expression of this gene was inhibited by cold treatment, while it was induced in leaves (Figure 5).

Basic helix-loop-helix (bHLH)-type transcription factors played important regulatory roles in diverse biological processes in plants (Toledoortiz et al., 2003). The latest evidences showed that bHLH transcription factor acted as positive regulators of CBF-pathway and conferred cold tolerance in plants (Feng et al., 2012; Peng et al., 2013; Xu et al., 2014). We found that the expression of four and four bHLH genes was altered after cold treatment in Thellungiella roots and leaves, respectively (Table 3 and Supplementary Tables S1, S2). However, only one of these bHLH gene (Thhalv10025656m) was induced by cold, and the others were down-regulated in response to cold treatment. Functional annotation showed that Thhalv10025656m was homologous of bHLH69 of Arabidopsis, and named as TsbHLH69. In Arabidopsis, bHLH69 contributed to the regulation of circadian periodicity by reducing the expression of LHY and TOC1 (Hanano et al., 2008; Harmer, 2009). These data suggested that the cold-induced TsbHLH69 might participate in the regulation of rhythm of Thellungiella under cold condition.

In this study, we found that the expression of many genes related to plant hormone biosynthesis and signaling were altered upon cold treatment (Supplementary Tables S1, S2). A total of 134 DEGs were related to plant hormone biosynthesis and signaling in leaves, 71 of which were up-regulated and 73 of which were down-regulated. DELLA, an key factor in gibberellins (GA) signaling, was also involved in the signal transduction of other hormones suggesting that DELLA functions as a modulator of plant development and response to stresses (Achard et al., 2003; Willige et al., 2007). Studies showed that DELLA contributed to CBF1-induced cold acclimation and was considered as components of CBF1-mediated cold stress response (Achard et al., 2008). We found the expression of a DELLA encoding gene (Thhalv10015535m) was slightly induced in both roots (1.40 fold) and leaves (1.39 fold) upon cold treatment. Meanwhile, gibberellin 2-oxidase (Thhalv10008179m), an enzyme inactivating the bioactive gibberellins (GAs), was significantly induced (2.6-fold) (Supplementary Table S2). These results implied that GA metabolism and signaling might contribute to cold stress tolerance in Thellungiella.

Studies demonstrated that the expression of COR genes was strongly induced after plants were shifted to cold temperature (Hajela et al., 1990; Griffith et al., 2007). RNA-seq results revealed five COR genes including COR27, COR47, and three COR15 were dramatically induced in both roots and leaves under cold treatment (Table 4). Interestingly, we found that the expression levels of all three COR15 increased by 1,000-fold in leaves after cold treatment. Among all these COR genes, COR15 represented the most induced genes, whose expression level increased more than 8,000-fold to compare with the control. Sequence alignment results showed that the DNA sequences of these three COR15 genes were highly similar, implying these genes might have the same functions in cold tolerance (Supplementary Figure S1). The COR47 gene encoding a member of dehydrin has been isolated in several plant species. Previous studies showed that the promoter region of COR47 contained C-repeat, dehydration-responsive element, low temperature-responsive element (CRT/DRE/LTREs) and ABA regulatory element (ABRE). COR47 was the downstream gene of CBF/DREB1, and CBF/DREB1 binds to the promoter of COR47 to induce its expression (Jaglo-Ottosen et al., 1998). Functional annotation and sequence alignment showed that Thhalv10027405m, a TsCBF/DREB1 gene, was a homolog of Arabidopsis CBF1/2/3. RNA-seq data revealed that Thhalv10027405m was up-regulated in both roots and leaves under cold treatment (Table 4 and Supplementary Table S2). The expression profiles of TsCBF/DREB1 and COR47 in Thellungiella was consistent with their counterparts in Arabidopsis, suggesting that CBF/DREB1 and COR47 genes might play important roles in DREB1/CBF-dependent cold signaling pathway in Thellungiella.

In addition, our RNA-seq data indicated that there might be more genes involved in DREB1/CBF pathway in Thellungiella. For example, we found DREB2B gene (Thhalv10021161m) was up-regulated in leaves under cold treatment (Supplementary Table S2). In addition to COR47 gene, we observed other two dehydrin genes, Thhalv10010821m and Thhalv10008313m were also induced by cold treatment (Figure 5 and Supplementary Table S2).

In Arabidopsis, INDUCER OF CBF EXPRESSION 1 (ICE1) could induce CBF3 expression by binding to its promoter. It was an important upstream regulatory factor in DREB1/CBF signaling pathway (Chinnusamy et al., 2003). In Thellungiella, we observed that the expression of ICE1 (Thhalv10004059m) were up-regulated in roots upon cold treatment (1.57-fold). However, the expression level of ICE1 was similar in cold treated and control samples in leaves. These results indicated that the function of ICE1 might be different in Arabidopsis and Thellungiella.

Aquaporins (or water channel proteins) played a crucial role in plant water relations. According to their distinct sub-cellular localization, aquaporins could be divided into four subgroups, the tonoplast intrinsic proteins (TIPs), plasma membrane intrinsic proteins (PIPs), Nod26-like intrinsic membrane proteins (NIPs) and small basic intrinsic proteins (SIPs). TIPs were abundantly expressed in vacuole (Maurel, 2007) and played an important role in maintenance of the intracellular space by controlling the water influx in vacuole (Leitão et al., 2014). A member of tonoplast intrinsic proteins GhTIP1;1 was responsive to cold stress and contributed to freezing-tolerance in cotton (Li D.D. et al., 2009). In the current study, five cold regulated aquaporin genes including three TIPs, one PIPs and one NIPs were identified. RNA-seq and qRT-PCR results demonstrated that TIP2-1 (Thhalv10021303m) was significantly down-regulated in both roots and leaves when the plants were under cold treatment. Another two TIP genes, Thhalv10011730m (also annotated as TIP2-1) and Thhalv10002084m (TIP4-1), were also significantly down-regulated upon cold treatment. These results suggested that the decreased expression of TIP genes might be beneficial to reduce water in/out of vacuole which is important for maintaining the stability of the cells in cold condition (Figure 5 and Supplementary Tables S1, S2). Moreover, we observed that the expression of PIP2-7 (Thhalv10025940m) was significantly decreased, suggesting that the water in/out of the cell might be reduced.

Seeds of T. salsuginea (Shandong ecotype) were surface-sterilized and plated on 1/2 MS-agar plates for synchronize germination at 4°C for a week. The plates were moved to growth chamber with 16 h light at 26°C with light intensity of 3000 lx and 8 h dark at 22°C. Seven-day-old seedlings were transferred to soil for 5 weeks. For cold treatment, seedlings were exposed to a growth chamber under a 16/8 h light/dark regime at 8/4°C for 24 h. For the control, the seedlings were grown at 22/16°C. Leaves and roots of control and cold treated seedlings were collected, frozen immediately in liquid nitrogen, and stored at -80°C. To minimize the plant to plant variation, nine individuals were used as one biological replicate, the tissue from nine individuals were pooled into one independent biological replicate. For both leaves and roots, three biological replicates were prepared.

Total RNA was isolated from leaves and roots using RNAiso Reagent (Takara, China), and treated with DNase I (Takara, China) to remove the contaminated genomic DNA according to the manufacturer’s protocols. RNA quality was detected by electrophoresis on 1.0% agarose gels and NanoDrop. The mRNA was enriched and cleaved into short fragments (about 200 nt). The mRNA fragments were used as templates to synthesize the first strand cDNA using random hexamer-primer. The first strand cDNA was further incubated with DNA polymerase I, buffer, dNTPs and RNase H to synthesize the second strand. Following end repair, a single nucleotide (adenine) was added, and then sequencing adaptors were ligated to the fragments. Finally, the fragments were purified and enriched with PCR amplification to construct cDNA library. Two biological replicates of each sample were used for RNA-Seq via Illumina HiSeq^TM2000 platform by Beijing Genomics Institute (BGI).

To acquire clean reads, the low-quality reads, adaptor sequences, and empty reads were removed. All clean reads were mapped with the genome sequences of T. salsuginea² using SOAP2 program under the criterion of no more than two mismatches in the alignment (Li R. et al., 2009). The gene expression level was calculated using RPKM (Reads per Kb per Million reads) method according to previous studies (Morrissy et al., 2009; Qu et al., 2016). The relative gene expression level between different samples was calculated by log₂ ratio. Differentially expressed genes (DEGs) were identified using NOIseq under the criteria of probability ≥0.8 and the absolute value of log²Ratio ≥ 1. The probability (P-value) was calculated according to the manufacturer’s protocol with the default parameters (Tarazona et al., 2011).

Gene Ontology annotation was conducted using Blast2GO (Conesa et al., 2005) by comparing all DEGs with GO terms in the database, which covered three domains: cellular component, molecular function and biological process³. The significantly enriched GO terms in DEGs were identified using hypergeometric test comparing to the genome background under the standard of p-value ≤ 0.05. Then GO picture was generated using WEGO⁴ (Web Gene Ontology Annotation Plot) (Ye et al., 2006). KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway analysis was performed by mapping the DEGs to specific biochemical pathways in KEGG database⁵. Significantly enriched metabolic pathways or signal transduction pathways were identified using enrichment analysis by comparing DEGs with the whole genome background.

The primers used for qRT-PCR validation were designed using primer premier 5.0 software⁶ and were listed in Supplementary Table S5. The primers were designed according to the transcript sequences of T. salsuginea downloaded from database ⁷. qRT-PCR was performed in ABI7500 Real Time System (Applied Biosystems) using SYBR Green I (Roche) with the following reaction: 94°C 10 min; 94°C 15 s, 60°C 10 s and 72°C 25 s for 40 cycles. All reactions were performed in biological triplicates, ubiquitin, and actin were used as internal reference genes. The relative expression of genes was calculated by the software of ABI7500 Real-Time PCR System using the 2^{-△△ Ct} method.