# A De Novo Transcriptome Analysis Identifies Cold-Responsive Genes in the Seeds ofTaxillus chinensis (DC.) Danser

Authors: Jine Fu; Lingyun Wan; Lisha Song; Lili He; Ni Jiang; Hairong Long; Juan Huo; Xiaowen Ji; Fengyun Hu; Shugen Wei; Limei Pan
Journal: BioMed Research International
 (2022)

---

## Abstract

Taxillus chinensis (DC.) Danser, a parasitic plant of the Loranthaceae family, grows by attacking other plants. It has a long history of being used in Chinese medicine to treat multiple chronic diseases. We previously observed that T. chinensis seeds are sensitive to cold. In this study, we performed transcriptome sequencing for T. chinensis seeds treated by cold (0°C) for 0 h, 12 h, 24 h, and 36 h. TRINITY assembled 257,870 transcripts from 223,512 genes. The GC content and N50 were calculated as 42.29% and 1,368, respectively. Then, we identified 42,183 CDSs and 35,268 likely proteins in the assembled transcriptome, which contained 1,622 signal peptides and 6,795 transmembrane domains. Next, we identified 17,217 genes (FPKM>5) and 2,333 differentially expressed genes (DEGs) in T. chinensis seeds under cold stress. The MAPK pathway, as an early cold response, was significantly enriched by the DEGs in the T. chinensis seeds after 24 h of cold treatment. Known cold-responsive genes encoding abscisic acid-associated, aquaporin, C-repeat binding factor (CBF), cold-regulated protein, heat shock protein, protein kinase, ribosomal protein, transcription factor (TF), zinc finger protein, and ubiquitin were deregulated in the T. chinensis seeds under cold stress. Notably, the upregulation of CBF gene might be the consequences of the downregulation of MYB and GATA TFs. Additionally, we identified that genes encoding CDC20, YLS9, EXORDIUM, and AUX1 and wound-responsive family protein might be related to novel mechanisms of T. chinensis seeds exposed to cold. This study is first to report the differential transcriptional induction in T. chinensis seeds under cold stress. It will improve our understanding of parasitic plants in response to cold and provide a valuable resource for future studies.


---

## Body

## 1. Introduction

Taxillus chinensis (DC.) Danser, also named “Sang Ji Sheng,” is a parasitic plant from the Loranthaceae family and grows by attacking other plants, such as Aceraceae, Anacardiaceae, Euphorbiaceae, Fabaceae, Fagaceae, Juglandaceae, Moraceae, Rosaceae, and Rutaceae [1]. Because the leaves and stems of T. chinensis have been used to treat angina pectoris, arrhythmia, hypertension, rheumatism, stroke, and threatened abortion [1], it has a long history of being used in the Chinese medicine. It is mainly distributed in the southern and southwestern areas of China, probably due to the warm and humid climate. Our knowledge of T. chinensis is very limited, especially its responses to abiotic stresses.Among the abiotic stresses, temperature is an important factor that affects the plant physiological processes [2]. Cold stress, including chilling (<20°C) and freezing (<0°C), exerts a substantial effect on the plant health and potentially limits the plant growth, development, yield, and geographical distribution [3]. As a mechanism to combat low environmental temperature, plants have developed a series of response mechanisms to adjust gene expression and further to enhance their cold tolerance [4]. In recent decades, many studies have been elucidated the molecular mechanisms involved in the plant cold acclimation. CBF/DREB- (C-repeat binding factors/dehydration factor/dehydration responsive element binding factor-) dependent signalling was characterized as a key and conserved regulatory mechanism of many plants in response to cold, together with the CBF activators (e.g., ICE1, CAMTA3, and BZR1/BES1) and repressors (e.g., MYB15, PIFs, and EIN3) [5, 6]. By binding to the cis-acting elements, CBFs induce the expression of numerous downstream cold-responsive (COR) genes [3]. Upon the exposure to cold, some protective genes are transcribed to enhance the cold tolerance of plants, such as heat shock proteins (HSPs) [7], protein kinases (PKs; e.g., SNF1-related protein kinases 2.6/open stomata 1) [8], cryoprotective proteins [4], zinc finger proteins [6], and various metabolites [9]. Additionally, many transcription factors (TFs) have been reported to regulating gene expression in plants under cold stress, such as ethylene response factors [10], NACs [11], MYBs [12], bHLHs [13], and WRKYs [14]. In addition, CBF-independent regulatory pathways have also been characterized in plants to enhance their cold tolerance, including the plant hormones auxin, abscisic acid, ethylene, gibberellins, and jasmonic acid [15]. However, little is known about the molecular mechanisms activated in T. chinensis seeds in response to cold.Transcriptome sequencing, a next-generation sequencing technology, is an efficient method to detect gene expression profiles and elucidate the breadth of molecular mechanisms involved in many physiological processes [16]. It has been widely used to identify key genes and factors involved in the responses to abiotic/biotic stresses in plants due to its advantages in the large-scale functional assignment of genes, thorough qualitative and quantitative analyses of gene expression, improved sensitivity, and accurate profiling of eukaryotic transcriptomes for both model and nonmodel organisms [17, 18]. Importantly, it facilitates analyses of gene expression in organisms whose genomes are not accessible. For example, using transcriptome sequencing and de novo assembly analysis, Liu et al. investigated the transcriptome of Rumex patientia during cold stress and identified 66 genes that are putatively involved in the response to cold stress, including members of the MYB, AP2/ERF, CBF, Znf, bZIP, NAC, and COR families [19]. Mohamed Sathik et al. employed transcriptome sequencing to study abiotic stress responsive genes in Hevea brasiliensis [20]. Fu et al. analyzed the transcriptome of Elymus nutans under cold stress and identified 26 hub genes playing a central role in the response to cold [21]. Similar to these studies, transcriptome sequencing technology will enable us to study changes in gene expression changes in T. chinensis seeds under cold stress.Previously, our lab identified genes that are expressed in response to water loss inT. chinensis seeds [22], such as RD22, HSP, and various TFs (MYB, WRKY, and ethylene-responsive transcription factors), and reported the regulatory miRNAs in T. chinensis seeds in response to cold [23], such as miR408, miR393b, miR946, ath-miR779.2, miR398, and miR9662. Interestingly, ICE3, IAA13, and multiple TFs (e.g., WRKY and CRF4 and TCP4) were shown to be targets of the dysregulated miRNAs identified in the T. chinensis seeds under cold stress [23]. In the present study, we used the same material as the miRNA study and investigated changes in gene expression changes in T. chinensis seeds in response to cold using transcriptome sequencing technology. This study is the first to analyze the changes in gene expression in T. chinensis seeds under cold stress, and our results will improve our understanding of molecular mechanisms of cold stress in T. chinensis seeds.


## 2. Materials and Methods

### 2.1. Sample Collection and Cold Treatment

The original seeds ofT. chinensis were obtained from three mulberry trees in the wild, and no permissions were required to collect these samples. The seeds were then confirmed by a senior botanist and deposited in the herbarium of Guangxi Botanical Garden of Medicinal Plants (acc. S0001794). For the cold treatment experiment, we selected 300 T. chinensis seeds with similar appearances, including sizes, weights, and health conditions. We observed that the T. chinensis seeds were sensitive to temperature and that 0°C was a suitable temperature to study the cold-responsive genes [23]. The seeds were divided into three groups—no treatment (A0) and cold treatment for 12 h (A1), 24 h (A2), and 36 h (A3). Then, we examined the viability of A0, A1, A2, and A3 by immersing the seeds in a solution of 1% (w/v) 2,3,5-triphenyl tetrazolium chloride (TTC, Sigma), as previously described [22].


### 2.2. Total RNA Isolation, cDNA Library Construction, and Deep Sequencing

Total RNA was extracted from theT. chinensis seeds (100 mg) exposed to cold for 0 h (A0), 12 h (A1), 24 h (A2), and 36 h (A3) using TRIzol reagent, as previously described [22]. Next, an Agilent 2100 Bioanalyzer (Agilent Technologies) was used to determine the quantity and quality of the total RNA. Equal amounts of total RNA (1 μg) were used to construct the cDNA libraries for transcriptome sequencing. Briefly, mRNAs were enriched with magnetic oligo (dT) beads and fragmented into ~200 bp fragments, followed by the double strand cDNA library construction using random hexamer N6 primers. Next, the double strand cDNA libraries were end repaired by adding a phosphate at the 5′-end and sticky “A” to the 3′-end. After the sequencing primers were ligated to each library and multiple libraries were pooled using the index technology, the pooled library was sequenced on the BGISEQ-500 RS platform with a paired-end 150 (PE150) strategy at BGI-Shenzhen.


### 2.3. Read Cleaning and De Novo Assembly of the Transcriptome

Raw reads were processed using SOAPnuke to remove the low-quality reads, reads with adaptors, and contaminating reads [24]. The obtained clean reads were quality controlled using FASTQC, as previously described [25]. Then, the clean reads of all samples assembled into the transcriptome using the TRINITY software with default parameters, according to a published protocol [26]. We next aligned all the clean reads to the assembled transcriptome using Bowtie2 and determined the global gene expression profile using the RSEM (RNA-Seq by Expectation-Maximization) method [27, 28]. The fragments per million reads per kilo base mapped (FPKM) method was used to normalize gene expression. We filtered genes expressed at low levels (FPKM<1) from the assembled transcriptome for quality control of the assembled transcriptome. CD-HIT was used to cluster the assembled T. chinensis genes [29]. BUSCO was used to evaluate the completeness of the assembled transcriptome [30].


### 2.4. Transcriptome Annotation

We annotated the assembledT. chinensis transcriptome by mapping it to public databases. In detail, the transcriptome was aligned to the NT database using BLASTn and aligned to the COG, KEGG, NR, and SwissProt databases using BLASTx. Then, open reading frames (ORFs) were predicted in the assembled transcriptome and aligned to the InterPro database using InterProScan. Gene Ontology annotations were retrieved from the mapping results of InterPro and NR. KEGG pathway annotation results were retrieved from the KEGG pathway mapping results. In combination with the BLAST results, we fetched the possible coding sequence (CDS) regions in the assembled transcriptome. Likely proteins encoded by the assembled transcriptome were predicted by TransDecoder and were used to identify signal peptides and transmembrane helices with SignalP (v5.0b) and TMHMM (v2.0c), respectively, according to the manufacturers’ protocols.


### 2.5. Gene Expression Profile and Differential Expression Analysis

We aligned the clean reads of each sample to the assembled transcriptome using Bowtie2 to profile gene expression in theT. chinensis seeds under cold stress [27]. Then, the RSEM method was used to count the reads aligned to each gene and normalize the gene expression using the FPKM method [28]. We used edgeR to calculate the p values and false discovery rate (FDR) to identify differentially expressed genes in the T. chinensis seeds under cold stress. Some cut-offs were used to identify the DEGs, including log2 foldchangelog2FC>1 or < -1, p value < 0.05, coefficient of variationCV<0.7, and false discovery rateFDR<0.1.


### 2.6. Functional Analysis

Using the GO and KEGG pathway annotations for the assembledT. chinensis seed transcriptome, we analyzed the enriched GO terms (biological processes, molecular functions, and cellular components) and KEGG pathways of the DEGs. Initially, p values were calculated using the Fisher’s exact test to show the significance, and q values were calculated using the R package “qvalue” for quality control of the p values.


### 2.7. qRT-PCR

We selected 11 genes for qRT-PCR validation, and 18S rRNA was used as an internal control. Forward and reverse primers were predicted using Primer3 and synthesized at BGI (Shenzhen) (TableS1). The procedure for the qRT-PCR experiment was the same as that in our previous study [22]. Then, we calculated the ΔCt value to determine the expression levels of target genes in a sample and ΔΔCt value to assess the difference in gene expression between two samples. We used relative normalized expression (RNE) to show the gene expression changes: RNE=2−∆∆Ct and calculated its log2 value (log2RNE) to compare the changes in gene expression identified using RNA-Seq and qRT-PCR.


## 2.1. Sample Collection and Cold Treatment

The original seeds ofT. chinensis were obtained from three mulberry trees in the wild, and no permissions were required to collect these samples. The seeds were then confirmed by a senior botanist and deposited in the herbarium of Guangxi Botanical Garden of Medicinal Plants (acc. S0001794). For the cold treatment experiment, we selected 300 T. chinensis seeds with similar appearances, including sizes, weights, and health conditions. We observed that the T. chinensis seeds were sensitive to temperature and that 0°C was a suitable temperature to study the cold-responsive genes [23]. The seeds were divided into three groups—no treatment (A0) and cold treatment for 12 h (A1), 24 h (A2), and 36 h (A3). Then, we examined the viability of A0, A1, A2, and A3 by immersing the seeds in a solution of 1% (w/v) 2,3,5-triphenyl tetrazolium chloride (TTC, Sigma), as previously described [22].


## 2.2. Total RNA Isolation, cDNA Library Construction, and Deep Sequencing

Total RNA was extracted from theT. chinensis seeds (100 mg) exposed to cold for 0 h (A0), 12 h (A1), 24 h (A2), and 36 h (A3) using TRIzol reagent, as previously described [22]. Next, an Agilent 2100 Bioanalyzer (Agilent Technologies) was used to determine the quantity and quality of the total RNA. Equal amounts of total RNA (1 μg) were used to construct the cDNA libraries for transcriptome sequencing. Briefly, mRNAs were enriched with magnetic oligo (dT) beads and fragmented into ~200 bp fragments, followed by the double strand cDNA library construction using random hexamer N6 primers. Next, the double strand cDNA libraries were end repaired by adding a phosphate at the 5′-end and sticky “A” to the 3′-end. After the sequencing primers were ligated to each library and multiple libraries were pooled using the index technology, the pooled library was sequenced on the BGISEQ-500 RS platform with a paired-end 150 (PE150) strategy at BGI-Shenzhen.


## 2.3. Read Cleaning and De Novo Assembly of the Transcriptome

Raw reads were processed using SOAPnuke to remove the low-quality reads, reads with adaptors, and contaminating reads [24]. The obtained clean reads were quality controlled using FASTQC, as previously described [25]. Then, the clean reads of all samples assembled into the transcriptome using the TRINITY software with default parameters, according to a published protocol [26]. We next aligned all the clean reads to the assembled transcriptome using Bowtie2 and determined the global gene expression profile using the RSEM (RNA-Seq by Expectation-Maximization) method [27, 28]. The fragments per million reads per kilo base mapped (FPKM) method was used to normalize gene expression. We filtered genes expressed at low levels (FPKM<1) from the assembled transcriptome for quality control of the assembled transcriptome. CD-HIT was used to cluster the assembled T. chinensis genes [29]. BUSCO was used to evaluate the completeness of the assembled transcriptome [30].


## 2.4. Transcriptome Annotation

We annotated the assembledT. chinensis transcriptome by mapping it to public databases. In detail, the transcriptome was aligned to the NT database using BLASTn and aligned to the COG, KEGG, NR, and SwissProt databases using BLASTx. Then, open reading frames (ORFs) were predicted in the assembled transcriptome and aligned to the InterPro database using InterProScan. Gene Ontology annotations were retrieved from the mapping results of InterPro and NR. KEGG pathway annotation results were retrieved from the KEGG pathway mapping results. In combination with the BLAST results, we fetched the possible coding sequence (CDS) regions in the assembled transcriptome. Likely proteins encoded by the assembled transcriptome were predicted by TransDecoder and were used to identify signal peptides and transmembrane helices with SignalP (v5.0b) and TMHMM (v2.0c), respectively, according to the manufacturers’ protocols.


## 2.5. Gene Expression Profile and Differential Expression Analysis

We aligned the clean reads of each sample to the assembled transcriptome using Bowtie2 to profile gene expression in theT. chinensis seeds under cold stress [27]. Then, the RSEM method was used to count the reads aligned to each gene and normalize the gene expression using the FPKM method [28]. We used edgeR to calculate the p values and false discovery rate (FDR) to identify differentially expressed genes in the T. chinensis seeds under cold stress. Some cut-offs were used to identify the DEGs, including log2 foldchangelog2FC>1 or < -1, p value < 0.05, coefficient of variationCV<0.7, and false discovery rateFDR<0.1.


## 2.6. Functional Analysis

Using the GO and KEGG pathway annotations for the assembledT. chinensis seed transcriptome, we analyzed the enriched GO terms (biological processes, molecular functions, and cellular components) and KEGG pathways of the DEGs. Initially, p values were calculated using the Fisher’s exact test to show the significance, and q values were calculated using the R package “qvalue” for quality control of the p values.


## 2.7. qRT-PCR

We selected 11 genes for qRT-PCR validation, and 18S rRNA was used as an internal control. Forward and reverse primers were predicted using Primer3 and synthesized at BGI (Shenzhen) (TableS1). The procedure for the qRT-PCR experiment was the same as that in our previous study [22]. Then, we calculated the ΔCt value to determine the expression levels of target genes in a sample and ΔΔCt value to assess the difference in gene expression between two samples. We used relative normalized expression (RNE) to show the gene expression changes: RNE=2−∆∆Ct and calculated its log2 value (log2RNE) to compare the changes in gene expression identified using RNA-Seq and qRT-PCR.


## 3. Results

### 3.1. Overview of Deep Sequencing Results and De Novo Assembly

We previously observed thatT. chinensis seeds were sensitive to temperature and that 0°C was a suitable temperature to study the cold responsive genes in T. chinensis seeds [23]. Thus, we performed paired-end transcriptome sequencing in triplicate for T. chinensis seeds exposed to 0°C for 0 h (A0), 12 h (A1), 24 h (A2), and 36 h (A3). Initially, we obtained 728.82 million raw reads and 725.59 million clean reads after data cleaning. Then, the clean reads of all samples were merged and used for the de novo assembly and analysis. TRINITY assembled 257,870 transcripts derived from 223,512 genes. After the expression levels of the assembled genes were estimated, we filtered genes at low levels (FPKM<1) from the assembled transcriptome and clustered the genes. Finally, we obtained 111,390 transcripts corresponding to 98,514 T. chinensis genes. The statistical values of the assembled transcriptome are shown in Table 1. The GC content and N50, a statistical measure of the average length of a sequence set, were calculated as 42.29% and 1,368, respectively, for the assembled T. chinensis transcriptome. Next, we evaluated the length distribution of the assembled T. chinensis transcriptome. As shown in Figure 1(a), 39,302 transcripts (39.89%) with lengths between 200 and 300 nt and 3,341 transcripts (1.25%) longer than 4,000 nt were obtained. The completeness of the assembled T. chinensis transcriptome was evaluated using BUSCO, which identified 249 (97.7%) complete (197 complete and single-copy and 52 complete and duplicated), 4 fragmented, and 2 missing BUSCOs.Table 1
Overview of transcriptome sequencing and de novo assembly of the T. chinensis seeds under cold stress.
A0A1A2A3A0_R1A0_R2A0_R3A1_R1A1_R2A1_R3A2_R1A2_R2A2_R3A3_R1A3_R2A3_R3Raw reads649190286561502065285442629602786286955064468970452237005316792653119914629519766489278863348464Clean reads648758986557470865243956629190446282569064424278444639425209393452094700629125906484780663310034Trinity assembly257870 transcripts and 223512 genesAfter filter (FPKM>1)147570 transcripts and 116545 genesAfter cluster111390 transcripts and 98514 genesGC (%)42.29N501368 for transcripts and 1008 for genesGenes (FPKM>5)164691403615750158531504113258161521445114654155961349815029Genes (averageFPKM>5)15414146581496314849Figure 1
De novo assembly transcriptome of T. chinensis seeds under cold stress and its annotation. (a) Length distribution of transcripts assembled for the T. chinensis seeds under cold stress. (b) Numbers of transcripts aligned to different known databases. (c) Number of transcripts aligned to genes from different species in NR. (d) COG function classification of the assembled transcriptome. (e) Gene Ontology annotation for the assembled transcriptome. (f) KEGG pathway annotation for the assembled transcriptome.
(a)(b)(c)(d)(e)(f)


### 3.2. Annotation of the AssembledT. chinensis Transcriptome

We next annotated the assembledT. chinensis transcriptome by mapping it to multiple databases. A total of 59,830 transcripts were annotated, of which 37,643, 52,000, 25,421, 26,725, 16,092, 21,442, and 26,866 were aligned to the NR, NT, SwissProt, KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway, COG (Clusters of Orthologous Groups), GO (Gene Ontology), and InterPro databases, respectively (Figure 1(b)). Notably, 10,509 transcripts were annotated by all these databases (Figure 1(b)). Among the NR mapping results, we found that the top four species hits to the assembled T. chinensis transcriptome were Vitis vinifera (10,313 transcripts), Theobroma cacao (1,980 transcripts), Nelumbo nucifera (1,828 transcripts), and Ziziphus jujuba (1,481 transcripts) (Figure 1(c)). COG annotation (Figure 1(d)) revealed 1,126, 730, 198, and 31 transcripts from “signal transduction mechanism,” “energy production and conversion,” “defense mechanism,” and “cell motility,” respectively. Among the 26,276 transcripts with GO annotations (Figure 1(e)), we found that the top four GO terms were “metabolic process” (13,530 transcripts), “cellular process” (12,603 transcripts), “catalytic activity” (10,471 transcripts), and “binding” (9,847 transcripts). Additionally, we identified 1,061 and 192 transcripts related to “signaling” and “signal transducer activity,” respectively (Figure 1(e)). By mapping the T. chinensis transcriptome to the KEGG pathway database, we identified transcripts from five categories (Figure 1(f)), including cellular processes, environmental information processing, genetic information processing, human diseases, metabolism, and organismal systems. Interestingly, metabolism was enriched in most of the transcripts, and 4,896 transcripts were annotated from the “global and overview maps” of metabolisms. We detected 833 and 855 transcripts involved in the pathways of “signal transduction” and “environmental adaption,” respectively (Figure 1(f)).Using the BLAST results, we identified 38,479 coding sequences (CDSs) from the assembledT. chinensis transcriptome, which consisted of 24.21 million bases (mean length: 629, N50: 927, and GC content: 46.95%). Next, ESTScan identified 3,704 CDSs from the T. chinensis transcriptome, whose size was 1.26 million bases (mean length: 339, N50: 330, and GC content: 50.74%) [31]. Thus, in total, we obtained 42,183 CDSs for the T. chinensis seeds under cold stress (25.46 million bases, mean length: 603, N50: 891, and GC content: 47.14). Furthermore, we predicted 35,268 likely proteins derived from 30,728 of the T. chinensis transcripts using TransDecoder (https://github.com/TransDecoder/TransDecoder). SignalP and TMHMM identified 1,622 signal peptides and 6,795 transmembrane domains among the likely proteins of T. chinensis seeds. Annotations obtained from different perspectives improved our understanding of the assembled transcriptome and were useful for the identification of cold-responsive genes in the T. chinensis seeds. However, further experiments are required to explore some transcripts that were annotated without encoding capacity.


### 3.3. Gene Expression Profiles of theT. chinensis Seeds under Cold Stress

The viability ofT. chinensis seeds decreased quickly under cold stress, and we investigated the cold-responsive genes of T. chinensis seeds. Read mapping showed that 81.52%-84.95% of the clean reads were aligned to the assembled transcriptome. Then, we obtained 17,236 genes with FPKM>5 in the T. chinensis seeds under cold stress, of which 15,414, 14,658, 14,963, and 14,849 were distributed in A0, A1, A2, and A3, respectively (Figure 2(a), Table 1, and Table S2). Pearson’s correlation analysis revealed that the linear correlation between these samples was greater than 0.97, indicating a limited difference between them. The distribution of gene expression revealed that 58.30%-58.68% of detected genes in the T. chinensis seeds were expressed with FPKM values ranging from 10 to 99 (Figure 2(b)). Notably, 64, 69, 74, and 83 T. chinensis genes were expressed at more than 1000 FPKM in A0, A1, A2, and A3, respectively.Figure 2
Gene expression profiles and differential expression analysis for theT. chinensis seeds under cold stress. (a) Venn diagram of genes identified in T. chinensis seeds in responsive to cold. (b) Distribution of T. chinensis genes from different expression levels. (c) Numbers of differentially expressed genes identified in T. chinensis seeds under cold stress. Biological processes involved by the DEGs identified in A1 (d), A2 (e), and A3 (f) compared to A0. KEGG pathway enrichment analysis of DEGs identified in A1 (g), A2 (h), and A3 (i) compared to A0.
(a)(b)(c)(d)(e)(f)(g)(h)(i)We next compared the highly expressed genes in all four samples (TableS2), and 9 of the top 10 highly expressed genes were shared by all four samples. Among them, TR50621|c2_g12 was specific to A0, and TR36565|c1_g1 (lipid transfer protein) was specific to A1, A2, and A3. Then, we counted the numbers of some known cold responsive genes. As shown in Table 2, in the T. chinensis seed transcriptome, we identified 18, 22, 3, 9, 36, 529, 327, 402, 82, and 297 genes related to abscisic acid, aquaporin, C-repeat binding factor (CBF), cold stress (e.g., cold-regulated (COR) protein), heat shock protein (HSP), protein kinase (PK), ribosomal protein (RP), transcription factor (TF), zinc finger protein (ZFP), and ubiquitin, respectively.Table 2
Cold-responsive genes of T. chinensis seeds.
CategoryTotalA1 vs. A0∗A2 vs. A0A3 vs. A0SharedAbscisic acid (ABA)190/1Aquaporin226/02/06/01/0C-repeat binding factor (CBF)31/01/01/01/0Cold stress related90/1Cold-regulated protein (COR)30/1Heat shock protein (HSP)363/0Protein kinase (PK)52926/216/040/1313/0Ribosomal protein (RP)3275/16/1Transcription factor (TF)40215/2512/731/327/5AP2, AP2-like, AP2/ERF172/01/03/01/0Ethylene-responsive (ER)545/106/512/143/4GATA180/30/4MYB290/31/4TCP110/1WRKY355/23/04/42/0bHLH320/22/2bZIP81/01/0Heat stress101/02/02/01/0Zinc finger protein (ZFP)827/75/28/83/2ZAT174/35/15/43/1Ubiquitin2973/32/06/111/0∗Numerator and denominator represent the upregulated and downregulated genes in the comparison.


### 3.4. Differentially Expressed Genes in theT. chinensis Seeds under Cold Stress

We next identified differentially expressed genes (DEGs) in theT. chinensis seeds under cold stress. Compared to A0, we detected 983 (538 upregulated and 445 downregulated), 451 (315 upregulated and 136 downregulated), and 1,857 (1,101 upregulated and 756 downregulated) DEGs in A1, A2, and A3, respectively (Figure 2(c) and Table S3). As shown in Figure 2(d), DEGs in A1 and A0 were involved in biological processes, including “GO:0009909~regulation of flower development,” “GO:0071840~cellular component organization or biogenesis,” “GO:0006807~nitrogen compound metabolic process,” “GO:0006096 ~ glycolytic process,” “GO:0046037~GMP metabolic process,” and “GO:0009697~salicylic acid biosynthetic process.” However, after 24 h of cold treatment, the DEGs between A2 and A0 were enriched in the biological processes “GO:0006810 ~ transport,” “GO:0008643 ~ carbohydrate transport,” “GO:0071555~cell wall organization,” “GO:0034757~negative regulation of iron ion transport,” and “GO:0034220~ion transmembrane transport” (Figure 2(e)). After 36 h of cold treatment, the DEGs were found to participate in the biological processes such as “GO:0006970~response to osmotic stress,” “GO:0009733~response to auxin,” and “GO:0006281 ~ DNA repair” (Figure 2(f)). Interestingly, cold-responsive genes of T. chinensis seeds were found to be enriched in different pathways during the cold treatment. For example, DEGs in A1 and A0 were enriched in 5 KEGG pathways (e.g., “ubiquitin mediated proteolysis,” “plant-pathogen interaction,” “oocyte meiosis,” “glutathione metabolism,” and “cell cycle”) (Figure 2(g)); “MAPK signaling pathway-plant” was triggered by the DEGs in T. chinensis seeds after 12 h of cold treatment (Figure 2(h)), and pathways such as “plant-hormone signal transduction,” “oxidative phosphorylation,” and “flavone and flavonol biosynthesis” were initiated in the T. chinensis seeds after 36 h of cold treatment (Figure 2(i)). The number of DEGs from different gene families, such as aquaporin, CBF, COR, PK, RP, TF, ZFP, and ubiquitin proteins, is shown in Table 2.


### 3.5. Cold-Responsive Genes inT. chinensis Seeds

Prior to identifying early and late cold-responsive genes inT. chinensis seeds, we showed the expression levels of all DEGs in these samples (A0, A1, A2, and A3) by constructing a heat map (Figure 3(a)). Venn diagrams (Figure 3(b)) revealed that 231 upregulated and 77 downregulated genes were shared in seeds during the cold treatment for 12 h, 24 h, and 36 h. Among them, 1 aquaporin, 1 ABA-associated gene, 13 HSPs, 12 TFs, 5 ZFPs, and 1 ubiquitin gene were commonly deregulated (Table 2). Notably, among the commonly upregulated genes detected in T. chinensis seeds under cold stress, we identified 5, 2, 3, and 3 shared upregulated genes encoding CDC20-1, CDC20-2, YLS9, and EXORDIUM, respectively (Table S3). In addition, 332 (151 upregulated and 181 downregulated) and 1,225 (707 upregulated and 518 downregulated) DEGs were specifically deregulated in T. chinensis seeds at the early (A1) and late (A3) stages of cold treatment, respectively (Figure 3(b)). Among the genes that were specifically deregulated in A1, we identified TR43917|c0_g1 encoding a wound-responsive family protein, TR37450|c1_g1 encoding auxin transporter protein 1 (AUX1), and three genes encoding aquaporins (Figure 3(c) and Table S3). These genes might be related to the early response to cold in T. chinensis seeds. Among the genes specifically deregulated in A3, we found 11 genes encoding disease resistance proteins, 4 genes encoding ER TFs, TR15040|c0_g1 encoding auxin-responsive protein IAA1-like, TR37606|c0_g4 encoding ethylene receptor, and TR39623|c0_g1 encoding ethylene response factor 10 (ERF10) (Figure 3(c) and Table S3), which might function in the late cold responsive stage in T. chinensis seeds.Figure 3
Differentially expressed genes inT. chinensis seeds exposed to cold. (a) Heat map of all DEGs identified in this study. (b) Venn diagram of DEGs identified in A1, A2, and A3 compared to A0. (c) Specifically upregulated genes T. chinensis seeds in response to cold at early and late stages.
(a)(b)(c)


### 3.6. Transcription Factors

We next further analyzed the TF expression patterns in theT. chinensis seeds under cold stress. Four hundred two TF genes were annotated in the T. chinensis seed transcriptome (Table 2), and 76 were differentially expressed in response to cold (Table 3). ER TF was the largest (27 genes) class that was differentially expressed in T. chinensis seeds under cold stress (Table 3 and Figure 4(a)). Notably, the expression levels of some ER TF genes were increased in the T. chinensis seeds at the late stage of cold treatment but not in A1 or A2. The downregulation of ER TF genes might be related to the seed development and other seed activities at normal temperature. Similar to ER TF genes, the functions of some other TF classes, such as WRKY and bHLH, are difficult to determine, as we observed both up- and downregulation of these genes in the T. chinensis seeds under cold stress (Figure 4(b)). However, the upregulation of AP2 domain class TF genes and the downregulation of GATA and MYB TF genes might participate in the cold response of the T. chinensis seeds, as they showed consistent regulation at all three time points of the T. chinensis seeds under cold stress (Table 3 and Figure 4(b)).Table 3
Differentially expressed TF genes inT. chinensis seeds under cold stress.
ClassGeneIDA1 vs. A0A2 vs. A0A3 vs. A0DescriptionLog2FCFDRRegulationLog2FCFDRRegulationLog2FCFDRRegulationAP2TR41417|c0_g11.651.05E-03Up1.431.43E-02Up1.647.33E-04UpAP2 domain-containing transcription factor family proteinTR45521|c5_g11.201.87E-02Up0.991.34E-01NC1.221.09E-02UpAP2 domain class transcription factorTR42026|c0_g10.752.68E-01NC0.607.29E-01NC1.122.53E-02UpAP2 domain class transcription factorERTR58370|c0_g10.268.85E-01NC0.666.05E-01NC1.633.74E-04UpEthylene-responsive transcription factor, putativeTR9006|c0_g1-1.241.44E-02Down-0.882.36E-01NC-1.122.33E-02DownEthylene-responsive transcription factor homolog, partialTR5872|c0_g10.209.29E-01NC-0.041.00E+00NC-1.133.50E-02DownEthylene-responsive transcription factor ERF118-likeTR47899|c0_g31.191.96E-02Up0.931.79E-01NC1.344.16E-03UpEthylene-responsive transcription factor ERF118TR8021|c0_g11.961.80E-04NC1.983.72E-04Up3.451.12E-12UpEthylene-responsive transcription factor ERF109-likeTR34720|c0_g22.481.22E-04NC3.131.10E-06NC5.292.62E-21NCEthylene-responsive transcription factor ERF109TR43918|c0_g11.084.21E-02Up1.011.19E-01NC1.815.26E-05UpEthylene-responsive transcription factor ERF105-likeTR44435|c0_g1-0.189.35E-01NC-0.391.00E+00NC-1.353.95E-03DownEthylene-responsive transcription factor ERF105TR15305|c0_g2-1.365.82E-03Down-1.213.53E-02Down-2.181.19E-06DownEthylene-responsive transcription factor ERF096TR52068|c0_g1-0.742.68E-01NC-0.587.41E-01NC-1.815.43E-05DownEthylene-responsive transcription factor ERF073TR78189|c0_g1-0.764.70E-01NC-0.411.00E+00NC-1.334.60E-02DownEthylene-responsive transcription factor ERF053TR29163|c0_g10.901.86E-01NC1.101.13E-01NC1.851.09E-04UpEthylene-responsive transcription factor ERF023TR101055|c0_g11.795.32E-04NC1.373.22E-02Up1.513.51E-03UpEthylene-responsive transcription factor ERF021TR15121|c0_g1-2.775.79E-08Down-1.761.33E-03Down-3.156.13E-10DownEthylene-responsive transcription factor ERF020TR63055|c0_g1-2.874.87E-10Down-1.504.15E-03Down-2.995.23E-11DownEthylene-responsive transcription factor ERF020TR43748|c0_g22.333.45E-07Up1.955.74E-05Up2.739.01E-10UpEthylene-responsive transcription factor ERF017-likeTR44877|c0_g11.047.20E-02Up1.156.95E-02Up2.171.94E-06UpEthylene-responsive transcription factor ERF014-likeTR40856|c0_g1-1.836.38E-05Down-1.135.60E-02Down-1.221.09E-02DownEthylene-responsive transcription factor ERF010-likeTR34995|c0_g1-1.133.80E-02Down-0.685.93E-01NC-1.132.98E-02DownEthylene-responsive transcription factor ERF003-likeTR28226|c0_g1-1.252.08E-02Down-0.745.11E-01NC-1.349.29E-03DownEthylene-responsive transcription factor ERF003-likeTR48171|c0_g1-0.496.08E-01NC-0.508.86E-01NC-1.151.86E-02DownEthylene-responsive transcription factor CRF5-likeTR18144|c0_g10.496.13E-01NC0.764.22E-01NC1.672.76E-04UpEthylene-responsive transcription factor CRF4-likeTR58185|c0_g10.941.04E-01NC1.184.33E-02Up2.302.19E-07UpEthylene-responsive transcription factor 4-likeTR47594|c1_g91.569.17E-03UP1.532.48E-02Up1.992.01E-04UpEthylene-responsive transcription factor 3-likeTR47594|c1_g120.978.25E-02NC0.902.10E-01NC1.451.65E-03UpEthylene-responsive transcription factor 3TR32067|c0_g1-1.152.64E-02Down-0.862.67E-01NC-1.334.88E-03DownEthylene-responsive transcription factor 2-likeTR147950|c0_g1-2.266.20E-05Down-1.591.68E-02Down-3.051.90E-07NCEthylene-responsive transcription factor 1BTR33516|c0_g1-1.543.72E-03Down-1.111.17E-01NC-2.491.17E-06DownEthylene-responsive transcription factor 1BGATATR117820|c0_g1-1.2345246830.037374327Down-0.873.95E-01NC-1.854.42E-04DownGATA transcription factor 4TR36890|c0_g1-1.2438573550.015695749Down-0.951.80E-01NC-1.625.01E-04DownGATA transcription factor 4TR44658|c0_g1-1.2029034890.020192595Down-0.734.71E-01NC-1.005.33E-02DownGATA transcription factor 21 isoform X1TR35784|c0_g1-0.747316940.257655438NC-0.734.62E-01NC-1.422.26E-03DownGATA transcription factor 1 isoform X2WRKYTR40062|c0_g3-0.2156526880.914437172NC-0.311.00E+00NC-1.201.39E-02DownWRKY8 transcription factorTR49705|c0_g1-1.416022260.003732422Down-0.852.81E-01NC-0.652.91E-01NCWRKY3 transcription factorTR63021|c0_g11.9643202745.63E-05Up1.216.65E-02Up1.182.83E-02NCWRKY-type transcription factorTR44228|c1_g11.7545162630.000398175Up1.912.14E-04Up3.015.62E-11UpWRKY transcription factor 70TR49520|c0_g100.716372310.288662641NC0.616.76E-01NC1.063.21E-02UpWRKY transcription factor 7TR46552|c0_g12.1214382354.29E-06Up1.679.78E-04Up2.322.21E-07UpWRKY transcription factor 51TR40322|c0_g1-0.8852883850.131822094NC-0.763.98E-01NC-1.392.76E-03DownWRKY transcription factor 44TR45024|c0_g11.1786456080.02156575Up0.783.77E-01NC1.004.83E-02UpWRKY transcription factor 41TR36422|c0_g31.3400141990.008307407Up0.892.61E-01NC0.938.74E-02NCWRKY transcription factor 33TR37159|c0_g1-0.5711322630.48525436NC-0.538.32E-01NC-1.014.85E-02DownWRKY transcription factor 31TR39273|c0_g1-1.3307667260.007657745Down-0.803.58E-01NC-1.845.22E-05DownWRKY transcription factor 21bHLHTR46437|c2_g1-0.3987195730.722357889NC-0.459.81E-01NC-1.258.73E-03DownTranscription factor bHLH96TR46303|c0_g10.6338416980.398045132NC0.616.85E-01NC1.604.63E-04UpTranscription factor bHLH47TR50792|c4_g150.769701540.240486872NC0.587.54E-01NC1.142.13E-02UpTranscription factor bHLH35-likeTR45007|c0_g1-1.1628829990.023954374Down-0.773.83E-01NC-1.201.23E-02DownTranscription factor bHLH147TR33541|c0_g1-1.0923967570.041192648Down-0.675.72E-01NC-0.861.13E-01NCTranscription factor bHLH117MYBTR68102|c0_g20.5099894820.679216055NC0.271.00E+00NC1.134.42E-02UpTranscription factor MYB44-likeTR22911|c0_g1-1.0360551970.057505766Down-0.714.93E-01NC-1.151.93E-02DownTranscription factor MYB44-likeTR39216|c0_g1-0.2978652590.852027342NC-0.361.00E+00NC-1.063.98E-02DownTranscription factor MYB12-likeTR33910|c0_g1-0.7165442130.294088684NC-0.636.36E-01NC-1.277.70E-03DownR3 MYB transcription factorTR44548|c0_g1-1.3103757750.009591119Down-0.764.23E-01NC-0.712.33E-01NCR2R3-MYB transcription factorTR36770|c0_g1-1.2540127290.01222635Down-0.842.87E-01NC-0.918.14E-02NCR2R3 Myb30 transcription factorTR41284|c1_g1-0.899449110.141069891NC-0.695.58E-01NC-1.327.06E-03DownMYB family transcription factor APL isoform X2OtherTR140849|c0_g10.2109812640.947640239NC-0.041.00E+00NC-1.099.18E-02DownTranscription factor TCP20-likeTR42062|c0_g12.310687851.32E-05NC3.073.14E-09NC4.454.54E-19UpTranscription factor perianthia-like isoform X2TR33703|c0_g2-0.4060945960.731737523NC-0.528.98E-01NC-1.172.43E-02DownTranscription factor IBH1TR23952|c0_g1-1.607281510.004248943Down-0.943.10E-01NC-0.653.70E-01NCTranscription factor HsfA5aTR9240|c0_g2-0.5498364040.519567814NC-1.213.78E-02Down-0.801.56E-01NCTranscription factor HEC2-likeTR92214|c0_g1-2.1393019514.02E-05Down-1.275.06E-02Down-2.181.61E-05DownTranscription factor HEC2-likeTR48713|c0_g10.6045663090.459558903NC0.607.42E-01NC1.201.60E-02UpTranscription factor GTE7-likeTR49896|c0_g11.0155929880.06218293Up0.882.33E-01NC1.353.79E-03UpTranscription factor GTE7TR78444|c0_g10.296476660.876314583NC0.441.00E+00NC1.085.33E-02UpTranscription factor GT-3b-likeTR49572|c1_g5-1.1183099090.03538728Down-0.587.54E-01NC-0.188.67E-01NCTranscription factor family proteinTR152976|c0_g1-0.5930079910.467699069NC-0.685.69E-01NC-1.932.77E-05DownScarecrow-like transcription factor PAT1TR38765|c0_g10.9807573980.098741136NC0.971.93E-01NC1.541.19E-03UpNucleolar transcription factor 1 (LOC104416354)TR26896|c0_g1-1.2671525820.013342976Down-0.411.00E+00NC0.583.76E-01NCNucleolar transcription factor 1TR52018|c3_g10.5110798940.60834591NC1.099.48E-02NC2.474.46E-08UpNAC transcription factors 36TR45153|c0_g12.0027340242.01E-05Up1.612.13E-03Up1.607.49E-04UpHeat stress transcription factor A-4bTR44136|c1_g10.4968208180.613845686NC1.438.21E-03Up3.076.23E-12UpHeat stress transcription factor A-2d isoform X1TR22732|c1_g11.1719580760.023399719Up0.842.91E-01NC1.082.99E-02UpbZIP transcription factor bZIP16TR48040|c0_g1-1.3507235030.008908786Down-0.744.87E-01NC-0.732.40E-01NCTranscription factor KAN4 isoform X2Figure 4
Expression patterns of some important gene families inT. chinensis seeds under cold stress. (a) Heat map of ER TF genes in T. chinensis seeds under cold stress. (b) Heat map of AP2, bHLH, GATA, MYB, and WRKY TF genes in T. chinensis seeds under cold stress. (c) Heat maps of protein kinase (left), ribosomal protein (middle), and zinc finger protein (right) genes in T. chinensis seeds under cold stress.
(a)(b)(c)


### 3.7. Protein Kinases, Ribosomal Proteins, and Zinc Finger Proteins

As shown in Table2, 60, 11, and 20 genes encoding PK, RP, and ZFP were differentially expressed in the T. chinensis seeds under cold stress, and their expression patterns are presented in Figure 4(c). Although some genes encoding calcium dependent, serine-threonine (e.g., LRR receptor-like and PTI1-like tyrosine) PKs were upregulated in A1 compared to A0; more PK genes from these subtypes were upregulated in A3 (Figure 4(c), left). Interestingly, we identified 402 RP genes in the T. chinensis seed transcriptome, but only 11 of them were deregulated in response to cold stress, including four 40S (e.g., RPS3 and RPS4) and seven 60S (e.g., RPL6, RPL7, and RPL10) RP genes. The overall expression of RP genes was upregulated in T. chinensis seeds to defend against cold stress (Figure 4(c), middle). We also detected both upregulated and downregulated ZFP genes in T. chinensis seeds after cold treatment (Figure 4(c), right). Two major classes of ZFP genes were identified and presented divergent expression patterns: dof ZFP and ZAT. The dof ZFP genes were downregulated, while the ZAT genes were upregulated in T. chinensis seeds under cold stress.


### 3.8. qRT-PCR

We used qRT-PCR to validate the expression patterns of genes in theT. chinensis seeds under cold stress. Eleven genes (e.g., CBF1, GATA4, MYB44, ERF4, ERF010, and ZAT10) and an internal control (18S rRNA) were selected to perform the qRT-PCR validation experiment. The primers for these candidate genes used in qRT-PCR are listed in Table S1. After ∆∆Ct values were calculated, we used log2 relative normalized expression (log2RNE) to show the changes in the expression of target genes between two samples, similar to the transcriptome sequencing. Among the 33 events showing the expression patterns of the target genes in T. chinensis seeds under cold stress, we found that 31 (93.9%) exhibited consistent changes in both RNA-Seq and qRT-PCR experiments (Figure 5). Notably, the expression patterns of CBF1, SOBIR1, and ERF10 were confirmed by qRT-PCR. The high agreement of gene expression patterns obtained using RNA-Seq and qRT-PCR indicates that the genes identified in this study might function in the cold response of T. chinensis seeds.Figure 5
qRT-PCR validation.


## 3.1. Overview of Deep Sequencing Results and De Novo Assembly

We previously observed thatT. chinensis seeds were sensitive to temperature and that 0°C was a suitable temperature to study the cold responsive genes in T. chinensis seeds [23]. Thus, we performed paired-end transcriptome sequencing in triplicate for T. chinensis seeds exposed to 0°C for 0 h (A0), 12 h (A1), 24 h (A2), and 36 h (A3). Initially, we obtained 728.82 million raw reads and 725.59 million clean reads after data cleaning. Then, the clean reads of all samples were merged and used for the de novo assembly and analysis. TRINITY assembled 257,870 transcripts derived from 223,512 genes. After the expression levels of the assembled genes were estimated, we filtered genes at low levels (FPKM<1) from the assembled transcriptome and clustered the genes. Finally, we obtained 111,390 transcripts corresponding to 98,514 T. chinensis genes. The statistical values of the assembled transcriptome are shown in Table 1. The GC content and N50, a statistical measure of the average length of a sequence set, were calculated as 42.29% and 1,368, respectively, for the assembled T. chinensis transcriptome. Next, we evaluated the length distribution of the assembled T. chinensis transcriptome. As shown in Figure 1(a), 39,302 transcripts (39.89%) with lengths between 200 and 300 nt and 3,341 transcripts (1.25%) longer than 4,000 nt were obtained. The completeness of the assembled T. chinensis transcriptome was evaluated using BUSCO, which identified 249 (97.7%) complete (197 complete and single-copy and 52 complete and duplicated), 4 fragmented, and 2 missing BUSCOs.Table 1
Overview of transcriptome sequencing and de novo assembly of the T. chinensis seeds under cold stress.
A0A1A2A3A0_R1A0_R2A0_R3A1_R1A1_R2A1_R3A2_R1A2_R2A2_R3A3_R1A3_R2A3_R3Raw reads649190286561502065285442629602786286955064468970452237005316792653119914629519766489278863348464Clean reads648758986557470865243956629190446282569064424278444639425209393452094700629125906484780663310034Trinity assembly257870 transcripts and 223512 genesAfter filter (FPKM>1)147570 transcripts and 116545 genesAfter cluster111390 transcripts and 98514 genesGC (%)42.29N501368 for transcripts and 1008 for genesGenes (FPKM>5)164691403615750158531504113258161521445114654155961349815029Genes (averageFPKM>5)15414146581496314849Figure 1
De novo assembly transcriptome of T. chinensis seeds under cold stress and its annotation. (a) Length distribution of transcripts assembled for the T. chinensis seeds under cold stress. (b) Numbers of transcripts aligned to different known databases. (c) Number of transcripts aligned to genes from different species in NR. (d) COG function classification of the assembled transcriptome. (e) Gene Ontology annotation for the assembled transcriptome. (f) KEGG pathway annotation for the assembled transcriptome.
(a)(b)(c)(d)(e)(f)


## 3.2. Annotation of the AssembledT. chinensis Transcriptome

We next annotated the assembledT. chinensis transcriptome by mapping it to multiple databases. A total of 59,830 transcripts were annotated, of which 37,643, 52,000, 25,421, 26,725, 16,092, 21,442, and 26,866 were aligned to the NR, NT, SwissProt, KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway, COG (Clusters of Orthologous Groups), GO (Gene Ontology), and InterPro databases, respectively (Figure 1(b)). Notably, 10,509 transcripts were annotated by all these databases (Figure 1(b)). Among the NR mapping results, we found that the top four species hits to the assembled T. chinensis transcriptome were Vitis vinifera (10,313 transcripts), Theobroma cacao (1,980 transcripts), Nelumbo nucifera (1,828 transcripts), and Ziziphus jujuba (1,481 transcripts) (Figure 1(c)). COG annotation (Figure 1(d)) revealed 1,126, 730, 198, and 31 transcripts from “signal transduction mechanism,” “energy production and conversion,” “defense mechanism,” and “cell motility,” respectively. Among the 26,276 transcripts with GO annotations (Figure 1(e)), we found that the top four GO terms were “metabolic process” (13,530 transcripts), “cellular process” (12,603 transcripts), “catalytic activity” (10,471 transcripts), and “binding” (9,847 transcripts). Additionally, we identified 1,061 and 192 transcripts related to “signaling” and “signal transducer activity,” respectively (Figure 1(e)). By mapping the T. chinensis transcriptome to the KEGG pathway database, we identified transcripts from five categories (Figure 1(f)), including cellular processes, environmental information processing, genetic information processing, human diseases, metabolism, and organismal systems. Interestingly, metabolism was enriched in most of the transcripts, and 4,896 transcripts were annotated from the “global and overview maps” of metabolisms. We detected 833 and 855 transcripts involved in the pathways of “signal transduction” and “environmental adaption,” respectively (Figure 1(f)).Using the BLAST results, we identified 38,479 coding sequences (CDSs) from the assembledT. chinensis transcriptome, which consisted of 24.21 million bases (mean length: 629, N50: 927, and GC content: 46.95%). Next, ESTScan identified 3,704 CDSs from the T. chinensis transcriptome, whose size was 1.26 million bases (mean length: 339, N50: 330, and GC content: 50.74%) [31]. Thus, in total, we obtained 42,183 CDSs for the T. chinensis seeds under cold stress (25.46 million bases, mean length: 603, N50: 891, and GC content: 47.14). Furthermore, we predicted 35,268 likely proteins derived from 30,728 of the T. chinensis transcripts using TransDecoder (https://github.com/TransDecoder/TransDecoder). SignalP and TMHMM identified 1,622 signal peptides and 6,795 transmembrane domains among the likely proteins of T. chinensis seeds. Annotations obtained from different perspectives improved our understanding of the assembled transcriptome and were useful for the identification of cold-responsive genes in the T. chinensis seeds. However, further experiments are required to explore some transcripts that were annotated without encoding capacity.


## 3.3. Gene Expression Profiles of theT. chinensis Seeds under Cold Stress

The viability ofT. chinensis seeds decreased quickly under cold stress, and we investigated the cold-responsive genes of T. chinensis seeds. Read mapping showed that 81.52%-84.95% of the clean reads were aligned to the assembled transcriptome. Then, we obtained 17,236 genes with FPKM>5 in the T. chinensis seeds under cold stress, of which 15,414, 14,658, 14,963, and 14,849 were distributed in A0, A1, A2, and A3, respectively (Figure 2(a), Table 1, and Table S2). Pearson’s correlation analysis revealed that the linear correlation between these samples was greater than 0.97, indicating a limited difference between them. The distribution of gene expression revealed that 58.30%-58.68% of detected genes in the T. chinensis seeds were expressed with FPKM values ranging from 10 to 99 (Figure 2(b)). Notably, 64, 69, 74, and 83 T. chinensis genes were expressed at more than 1000 FPKM in A0, A1, A2, and A3, respectively.Figure 2
Gene expression profiles and differential expression analysis for theT. chinensis seeds under cold stress. (a) Venn diagram of genes identified in T. chinensis seeds in responsive to cold. (b) Distribution of T. chinensis genes from different expression levels. (c) Numbers of differentially expressed genes identified in T. chinensis seeds under cold stress. Biological processes involved by the DEGs identified in A1 (d), A2 (e), and A3 (f) compared to A0. KEGG pathway enrichment analysis of DEGs identified in A1 (g), A2 (h), and A3 (i) compared to A0.
(a)(b)(c)(d)(e)(f)(g)(h)(i)We next compared the highly expressed genes in all four samples (TableS2), and 9 of the top 10 highly expressed genes were shared by all four samples. Among them, TR50621|c2_g12 was specific to A0, and TR36565|c1_g1 (lipid transfer protein) was specific to A1, A2, and A3. Then, we counted the numbers of some known cold responsive genes. As shown in Table 2, in the T. chinensis seed transcriptome, we identified 18, 22, 3, 9, 36, 529, 327, 402, 82, and 297 genes related to abscisic acid, aquaporin, C-repeat binding factor (CBF), cold stress (e.g., cold-regulated (COR) protein), heat shock protein (HSP), protein kinase (PK), ribosomal protein (RP), transcription factor (TF), zinc finger protein (ZFP), and ubiquitin, respectively.Table 2
Cold-responsive genes of T. chinensis seeds.
CategoryTotalA1 vs. A0∗A2 vs. A0A3 vs. A0SharedAbscisic acid (ABA)190/1Aquaporin226/02/06/01/0C-repeat binding factor (CBF)31/01/01/01/0Cold stress related90/1Cold-regulated protein (COR)30/1Heat shock protein (HSP)363/0Protein kinase (PK)52926/216/040/1313/0Ribosomal protein (RP)3275/16/1Transcription factor (TF)40215/2512/731/327/5AP2, AP2-like, AP2/ERF172/01/03/01/0Ethylene-responsive (ER)545/106/512/143/4GATA180/30/4MYB290/31/4TCP110/1WRKY355/23/04/42/0bHLH320/22/2bZIP81/01/0Heat stress101/02/02/01/0Zinc finger protein (ZFP)827/75/28/83/2ZAT174/35/15/43/1Ubiquitin2973/32/06/111/0∗Numerator and denominator represent the upregulated and downregulated genes in the comparison.


## 3.4. Differentially Expressed Genes in theT. chinensis Seeds under Cold Stress

We next identified differentially expressed genes (DEGs) in theT. chinensis seeds under cold stress. Compared to A0, we detected 983 (538 upregulated and 445 downregulated), 451 (315 upregulated and 136 downregulated), and 1,857 (1,101 upregulated and 756 downregulated) DEGs in A1, A2, and A3, respectively (Figure 2(c) and Table S3). As shown in Figure 2(d), DEGs in A1 and A0 were involved in biological processes, including “GO:0009909~regulation of flower development,” “GO:0071840~cellular component organization or biogenesis,” “GO:0006807~nitrogen compound metabolic process,” “GO:0006096 ~ glycolytic process,” “GO:0046037~GMP metabolic process,” and “GO:0009697~salicylic acid biosynthetic process.” However, after 24 h of cold treatment, the DEGs between A2 and A0 were enriched in the biological processes “GO:0006810 ~ transport,” “GO:0008643 ~ carbohydrate transport,” “GO:0071555~cell wall organization,” “GO:0034757~negative regulation of iron ion transport,” and “GO:0034220~ion transmembrane transport” (Figure 2(e)). After 36 h of cold treatment, the DEGs were found to participate in the biological processes such as “GO:0006970~response to osmotic stress,” “GO:0009733~response to auxin,” and “GO:0006281 ~ DNA repair” (Figure 2(f)). Interestingly, cold-responsive genes of T. chinensis seeds were found to be enriched in different pathways during the cold treatment. For example, DEGs in A1 and A0 were enriched in 5 KEGG pathways (e.g., “ubiquitin mediated proteolysis,” “plant-pathogen interaction,” “oocyte meiosis,” “glutathione metabolism,” and “cell cycle”) (Figure 2(g)); “MAPK signaling pathway-plant” was triggered by the DEGs in T. chinensis seeds after 12 h of cold treatment (Figure 2(h)), and pathways such as “plant-hormone signal transduction,” “oxidative phosphorylation,” and “flavone and flavonol biosynthesis” were initiated in the T. chinensis seeds after 36 h of cold treatment (Figure 2(i)). The number of DEGs from different gene families, such as aquaporin, CBF, COR, PK, RP, TF, ZFP, and ubiquitin proteins, is shown in Table 2.


## 3.5. Cold-Responsive Genes inT. chinensis Seeds

Prior to identifying early and late cold-responsive genes inT. chinensis seeds, we showed the expression levels of all DEGs in these samples (A0, A1, A2, and A3) by constructing a heat map (Figure 3(a)). Venn diagrams (Figure 3(b)) revealed that 231 upregulated and 77 downregulated genes were shared in seeds during the cold treatment for 12 h, 24 h, and 36 h. Among them, 1 aquaporin, 1 ABA-associated gene, 13 HSPs, 12 TFs, 5 ZFPs, and 1 ubiquitin gene were commonly deregulated (Table 2). Notably, among the commonly upregulated genes detected in T. chinensis seeds under cold stress, we identified 5, 2, 3, and 3 shared upregulated genes encoding CDC20-1, CDC20-2, YLS9, and EXORDIUM, respectively (Table S3). In addition, 332 (151 upregulated and 181 downregulated) and 1,225 (707 upregulated and 518 downregulated) DEGs were specifically deregulated in T. chinensis seeds at the early (A1) and late (A3) stages of cold treatment, respectively (Figure 3(b)). Among the genes that were specifically deregulated in A1, we identified TR43917|c0_g1 encoding a wound-responsive family protein, TR37450|c1_g1 encoding auxin transporter protein 1 (AUX1), and three genes encoding aquaporins (Figure 3(c) and Table S3). These genes might be related to the early response to cold in T. chinensis seeds. Among the genes specifically deregulated in A3, we found 11 genes encoding disease resistance proteins, 4 genes encoding ER TFs, TR15040|c0_g1 encoding auxin-responsive protein IAA1-like, TR37606|c0_g4 encoding ethylene receptor, and TR39623|c0_g1 encoding ethylene response factor 10 (ERF10) (Figure 3(c) and Table S3), which might function in the late cold responsive stage in T. chinensis seeds.Figure 3
Differentially expressed genes inT. chinensis seeds exposed to cold. (a) Heat map of all DEGs identified in this study. (b) Venn diagram of DEGs identified in A1, A2, and A3 compared to A0. (c) Specifically upregulated genes T. chinensis seeds in response to cold at early and late stages.
(a)(b)(c)


## 3.6. Transcription Factors

We next further analyzed the TF expression patterns in theT. chinensis seeds under cold stress. Four hundred two TF genes were annotated in the T. chinensis seed transcriptome (Table 2), and 76 were differentially expressed in response to cold (Table 3). ER TF was the largest (27 genes) class that was differentially expressed in T. chinensis seeds under cold stress (Table 3 and Figure 4(a)). Notably, the expression levels of some ER TF genes were increased in the T. chinensis seeds at the late stage of cold treatment but not in A1 or A2. The downregulation of ER TF genes might be related to the seed development and other seed activities at normal temperature. Similar to ER TF genes, the functions of some other TF classes, such as WRKY and bHLH, are difficult to determine, as we observed both up- and downregulation of these genes in the T. chinensis seeds under cold stress (Figure 4(b)). However, the upregulation of AP2 domain class TF genes and the downregulation of GATA and MYB TF genes might participate in the cold response of the T. chinensis seeds, as they showed consistent regulation at all three time points of the T. chinensis seeds under cold stress (Table 3 and Figure 4(b)).Table 3
Differentially expressed TF genes inT. chinensis seeds under cold stress.
ClassGeneIDA1 vs. A0A2 vs. A0A3 vs. A0DescriptionLog2FCFDRRegulationLog2FCFDRRegulationLog2FCFDRRegulationAP2TR41417|c0_g11.651.05E-03Up1.431.43E-02Up1.647.33E-04UpAP2 domain-containing transcription factor family proteinTR45521|c5_g11.201.87E-02Up0.991.34E-01NC1.221.09E-02UpAP2 domain class transcription factorTR42026|c0_g10.752.68E-01NC0.607.29E-01NC1.122.53E-02UpAP2 domain class transcription factorERTR58370|c0_g10.268.85E-01NC0.666.05E-01NC1.633.74E-04UpEthylene-responsive transcription factor, putativeTR9006|c0_g1-1.241.44E-02Down-0.882.36E-01NC-1.122.33E-02DownEthylene-responsive transcription factor homolog, partialTR5872|c0_g10.209.29E-01NC-0.041.00E+00NC-1.133.50E-02DownEthylene-responsive transcription factor ERF118-likeTR47899|c0_g31.191.96E-02Up0.931.79E-01NC1.344.16E-03UpEthylene-responsive transcription factor ERF118TR8021|c0_g11.961.80E-04NC1.983.72E-04Up3.451.12E-12UpEthylene-responsive transcription factor ERF109-likeTR34720|c0_g22.481.22E-04NC3.131.10E-06NC5.292.62E-21NCEthylene-responsive transcription factor ERF109TR43918|c0_g11.084.21E-02Up1.011.19E-01NC1.815.26E-05UpEthylene-responsive transcription factor ERF105-likeTR44435|c0_g1-0.189.35E-01NC-0.391.00E+00NC-1.353.95E-03DownEthylene-responsive transcription factor ERF105TR15305|c0_g2-1.365.82E-03Down-1.213.53E-02Down-2.181.19E-06DownEthylene-responsive transcription factor ERF096TR52068|c0_g1-0.742.68E-01NC-0.587.41E-01NC-1.815.43E-05DownEthylene-responsive transcription factor ERF073TR78189|c0_g1-0.764.70E-01NC-0.411.00E+00NC-1.334.60E-02DownEthylene-responsive transcription factor ERF053TR29163|c0_g10.901.86E-01NC1.101.13E-01NC1.851.09E-04UpEthylene-responsive transcription factor ERF023TR101055|c0_g11.795.32E-04NC1.373.22E-02Up1.513.51E-03UpEthylene-responsive transcription factor ERF021TR15121|c0_g1-2.775.79E-08Down-1.761.33E-03Down-3.156.13E-10DownEthylene-responsive transcription factor ERF020TR63055|c0_g1-2.874.87E-10Down-1.504.15E-03Down-2.995.23E-11DownEthylene-responsive transcription factor ERF020TR43748|c0_g22.333.45E-07Up1.955.74E-05Up2.739.01E-10UpEthylene-responsive transcription factor ERF017-likeTR44877|c0_g11.047.20E-02Up1.156.95E-02Up2.171.94E-06UpEthylene-responsive transcription factor ERF014-likeTR40856|c0_g1-1.836.38E-05Down-1.135.60E-02Down-1.221.09E-02DownEthylene-responsive transcription factor ERF010-likeTR34995|c0_g1-1.133.80E-02Down-0.685.93E-01NC-1.132.98E-02DownEthylene-responsive transcription factor ERF003-likeTR28226|c0_g1-1.252.08E-02Down-0.745.11E-01NC-1.349.29E-03DownEthylene-responsive transcription factor ERF003-likeTR48171|c0_g1-0.496.08E-01NC-0.508.86E-01NC-1.151.86E-02DownEthylene-responsive transcription factor CRF5-likeTR18144|c0_g10.496.13E-01NC0.764.22E-01NC1.672.76E-04UpEthylene-responsive transcription factor CRF4-likeTR58185|c0_g10.941.04E-01NC1.184.33E-02Up2.302.19E-07UpEthylene-responsive transcription factor 4-likeTR47594|c1_g91.569.17E-03UP1.532.48E-02Up1.992.01E-04UpEthylene-responsive transcription factor 3-likeTR47594|c1_g120.978.25E-02NC0.902.10E-01NC1.451.65E-03UpEthylene-responsive transcription factor 3TR32067|c0_g1-1.152.64E-02Down-0.862.67E-01NC-1.334.88E-03DownEthylene-responsive transcription factor 2-likeTR147950|c0_g1-2.266.20E-05Down-1.591.68E-02Down-3.051.90E-07NCEthylene-responsive transcription factor 1BTR33516|c0_g1-1.543.72E-03Down-1.111.17E-01NC-2.491.17E-06DownEthylene-responsive transcription factor 1BGATATR117820|c0_g1-1.2345246830.037374327Down-0.873.95E-01NC-1.854.42E-04DownGATA transcription factor 4TR36890|c0_g1-1.2438573550.015695749Down-0.951.80E-01NC-1.625.01E-04DownGATA transcription factor 4TR44658|c0_g1-1.2029034890.020192595Down-0.734.71E-01NC-1.005.33E-02DownGATA transcription factor 21 isoform X1TR35784|c0_g1-0.747316940.257655438NC-0.734.62E-01NC-1.422.26E-03DownGATA transcription factor 1 isoform X2WRKYTR40062|c0_g3-0.2156526880.914437172NC-0.311.00E+00NC-1.201.39E-02DownWRKY8 transcription factorTR49705|c0_g1-1.416022260.003732422Down-0.852.81E-01NC-0.652.91E-01NCWRKY3 transcription factorTR63021|c0_g11.9643202745.63E-05Up1.216.65E-02Up1.182.83E-02NCWRKY-type transcription factorTR44228|c1_g11.7545162630.000398175Up1.912.14E-04Up3.015.62E-11UpWRKY transcription factor 70TR49520|c0_g100.716372310.288662641NC0.616.76E-01NC1.063.21E-02UpWRKY transcription factor 7TR46552|c0_g12.1214382354.29E-06Up1.679.78E-04Up2.322.21E-07UpWRKY transcription factor 51TR40322|c0_g1-0.8852883850.131822094NC-0.763.98E-01NC-1.392.76E-03DownWRKY transcription factor 44TR45024|c0_g11.1786456080.02156575Up0.783.77E-01NC1.004.83E-02UpWRKY transcription factor 41TR36422|c0_g31.3400141990.008307407Up0.892.61E-01NC0.938.74E-02NCWRKY transcription factor 33TR37159|c0_g1-0.5711322630.48525436NC-0.538.32E-01NC-1.014.85E-02DownWRKY transcription factor 31TR39273|c0_g1-1.3307667260.007657745Down-0.803.58E-01NC-1.845.22E-05DownWRKY transcription factor 21bHLHTR46437|c2_g1-0.3987195730.722357889NC-0.459.81E-01NC-1.258.73E-03DownTranscription factor bHLH96TR46303|c0_g10.6338416980.398045132NC0.616.85E-01NC1.604.63E-04UpTranscription factor bHLH47TR50792|c4_g150.769701540.240486872NC0.587.54E-01NC1.142.13E-02UpTranscription factor bHLH35-likeTR45007|c0_g1-1.1628829990.023954374Down-0.773.83E-01NC-1.201.23E-02DownTranscription factor bHLH147TR33541|c0_g1-1.0923967570.041192648Down-0.675.72E-01NC-0.861.13E-01NCTranscription factor bHLH117MYBTR68102|c0_g20.5099894820.679216055NC0.271.00E+00NC1.134.42E-02UpTranscription factor MYB44-likeTR22911|c0_g1-1.0360551970.057505766Down-0.714.93E-01NC-1.151.93E-02DownTranscription factor MYB44-likeTR39216|c0_g1-0.2978652590.852027342NC-0.361.00E+00NC-1.063.98E-02DownTranscription factor MYB12-likeTR33910|c0_g1-0.7165442130.294088684NC-0.636.36E-01NC-1.277.70E-03DownR3 MYB transcription factorTR44548|c0_g1-1.3103757750.009591119Down-0.764.23E-01NC-0.712.33E-01NCR2R3-MYB transcription factorTR36770|c0_g1-1.2540127290.01222635Down-0.842.87E-01NC-0.918.14E-02NCR2R3 Myb30 transcription factorTR41284|c1_g1-0.899449110.141069891NC-0.695.58E-01NC-1.327.06E-03DownMYB family transcription factor APL isoform X2OtherTR140849|c0_g10.2109812640.947640239NC-0.041.00E+00NC-1.099.18E-02DownTranscription factor TCP20-likeTR42062|c0_g12.310687851.32E-05NC3.073.14E-09NC4.454.54E-19UpTranscription factor perianthia-like isoform X2TR33703|c0_g2-0.4060945960.731737523NC-0.528.98E-01NC-1.172.43E-02DownTranscription factor IBH1TR23952|c0_g1-1.607281510.004248943Down-0.943.10E-01NC-0.653.70E-01NCTranscription factor HsfA5aTR9240|c0_g2-0.5498364040.519567814NC-1.213.78E-02Down-0.801.56E-01NCTranscription factor HEC2-likeTR92214|c0_g1-2.1393019514.02E-05Down-1.275.06E-02Down-2.181.61E-05DownTranscription factor HEC2-likeTR48713|c0_g10.6045663090.459558903NC0.607.42E-01NC1.201.60E-02UpTranscription factor GTE7-likeTR49896|c0_g11.0155929880.06218293Up0.882.33E-01NC1.353.79E-03UpTranscription factor GTE7TR78444|c0_g10.296476660.876314583NC0.441.00E+00NC1.085.33E-02UpTranscription factor GT-3b-likeTR49572|c1_g5-1.1183099090.03538728Down-0.587.54E-01NC-0.188.67E-01NCTranscription factor family proteinTR152976|c0_g1-0.5930079910.467699069NC-0.685.69E-01NC-1.932.77E-05DownScarecrow-like transcription factor PAT1TR38765|c0_g10.9807573980.098741136NC0.971.93E-01NC1.541.19E-03UpNucleolar transcription factor 1 (LOC104416354)TR26896|c0_g1-1.2671525820.013342976Down-0.411.00E+00NC0.583.76E-01NCNucleolar transcription factor 1TR52018|c3_g10.5110798940.60834591NC1.099.48E-02NC2.474.46E-08UpNAC transcription factors 36TR45153|c0_g12.0027340242.01E-05Up1.612.13E-03Up1.607.49E-04UpHeat stress transcription factor A-4bTR44136|c1_g10.4968208180.613845686NC1.438.21E-03Up3.076.23E-12UpHeat stress transcription factor A-2d isoform X1TR22732|c1_g11.1719580760.023399719Up0.842.91E-01NC1.082.99E-02UpbZIP transcription factor bZIP16TR48040|c0_g1-1.3507235030.008908786Down-0.744.87E-01NC-0.732.40E-01NCTranscription factor KAN4 isoform X2Figure 4
Expression patterns of some important gene families inT. chinensis seeds under cold stress. (a) Heat map of ER TF genes in T. chinensis seeds under cold stress. (b) Heat map of AP2, bHLH, GATA, MYB, and WRKY TF genes in T. chinensis seeds under cold stress. (c) Heat maps of protein kinase (left), ribosomal protein (middle), and zinc finger protein (right) genes in T. chinensis seeds under cold stress.
(a)(b)(c)


## 3.7. Protein Kinases, Ribosomal Proteins, and Zinc Finger Proteins

As shown in Table2, 60, 11, and 20 genes encoding PK, RP, and ZFP were differentially expressed in the T. chinensis seeds under cold stress, and their expression patterns are presented in Figure 4(c). Although some genes encoding calcium dependent, serine-threonine (e.g., LRR receptor-like and PTI1-like tyrosine) PKs were upregulated in A1 compared to A0; more PK genes from these subtypes were upregulated in A3 (Figure 4(c), left). Interestingly, we identified 402 RP genes in the T. chinensis seed transcriptome, but only 11 of them were deregulated in response to cold stress, including four 40S (e.g., RPS3 and RPS4) and seven 60S (e.g., RPL6, RPL7, and RPL10) RP genes. The overall expression of RP genes was upregulated in T. chinensis seeds to defend against cold stress (Figure 4(c), middle). We also detected both upregulated and downregulated ZFP genes in T. chinensis seeds after cold treatment (Figure 4(c), right). Two major classes of ZFP genes were identified and presented divergent expression patterns: dof ZFP and ZAT. The dof ZFP genes were downregulated, while the ZAT genes were upregulated in T. chinensis seeds under cold stress.


## 3.8. qRT-PCR

We used qRT-PCR to validate the expression patterns of genes in theT. chinensis seeds under cold stress. Eleven genes (e.g., CBF1, GATA4, MYB44, ERF4, ERF010, and ZAT10) and an internal control (18S rRNA) were selected to perform the qRT-PCR validation experiment. The primers for these candidate genes used in qRT-PCR are listed in Table S1. After ∆∆Ct values were calculated, we used log2 relative normalized expression (log2RNE) to show the changes in the expression of target genes between two samples, similar to the transcriptome sequencing. Among the 33 events showing the expression patterns of the target genes in T. chinensis seeds under cold stress, we found that 31 (93.9%) exhibited consistent changes in both RNA-Seq and qRT-PCR experiments (Figure 5). Notably, the expression patterns of CBF1, SOBIR1, and ERF10 were confirmed by qRT-PCR. The high agreement of gene expression patterns obtained using RNA-Seq and qRT-PCR indicates that the genes identified in this study might function in the cold response of T. chinensis seeds.Figure 5
qRT-PCR validation.


## 4. Discussion

T. chinensis is a plant that is sensitive to cold, and this study was designed to investigate the changes in gene expression in T. chinensis seeds under cold stress. Due to the lack of genome sequences, we assembled the transcriptome for T. chinensis seeds using RNA-Seq, similar to the studies of other plants, such as Magnolia wufengensis [32], Passiflora edulis Sims [33], Hevea brasiliensis [20], and Rumex patientia [19]. We assembled 98,514 T. chinensis genes and 35,268 likely proteins expressed in the T. chinensis seeds in response to cold. Interestingly, many of the T. chinensis genes were predicted to have similarities with grape genes (Figure 1(c)). A number of cold-responsive genes are involved in multiple biological processes and pathways, such as plant hormone signal transduction [32, 33] and metabolic processes [19, 20, 32, 33]. The DEGs identified in the T. chinensis seeds under cold stress were determined to be enriched in these pathways (Figures 2(d)–2(i)). Some known cold-responsive genes, including ABA, aquaporin, CBF, COR, HSP, PK, RP, ZFP ubiquitin, and TFs, were identified (Table 2). In addition, we found that genes encoding CDC20-1, CDC20-2, YLS9, and EXORDIUM were upregulated in the T. chinensis seeds in response to cold (Table S3).Cold acclimation temperatures have the potential to induce profound changes in the plant transcriptome. Approximately 4% to 20% of the genes in theArabidopsis genome are affected by cold [34, 35]. Interestingly, different pathways were activated in the T. chinensis seeds during cold treatment (Figures 2(g)–2(i)). The “MAPK signalling pathway–plant” was significantly enriched in the DEGs in A2, but not A1, compared to A0 (Figure 2(h)). This pathway has been reported to modulate plant tolerance to multiple abiotic stresses, such as drought, salt, cold, and heat [36]. The MAPK cascades have been reported to convert the environmental stimuli into cellular responses and to negatively regulate freezing tolerance via the phosphorylation of ICE1, a basic-helix-loop-helix transcription factor that regulates the expression of CBF genes [37], which mediate cold-inducible transcription and play a key role in freezing tolerance and cold acclimation by binding to the C-repeat/DRE element [38]. MAPK pathway activation is a rapid response to cold, as MPK3, MPK4, and MPK6 are rapidly activated after cold treatment [37]. Compared to A0, we identified 10, 14, and 21 genes from the MAPK pathway that were differentially expressed in the T. chinensis seeds after exposure to cold for 12 h, 24 h, and 36 h, respectively (Figures 2(h) and 2(i)). Among them, MPK9 was commonly upregulated at all three time points while two MAPKK kinases (NPK1-like and YODA) were specifically upregulated in A3 (Table S3). NPK1-like and its clustered genes, namely, other NPK-like genes (e.g., OsNPKL2, 3, and 4), are induced by cold in rice [39]. YODA was shown to be upstream of MKK4/MKK5, a negative regulator of freezing tolerance [37], and downstream of the ER receptor in regulating coordinated local cell proliferation, which shapes the morphology of plant organs [40].Interestingly, the ICE1 TF gene (TR144797|c0_g1) was not dysregulated inT. chinensis seeds under cold stress; however, the CBF1 gene (TR16191|c0_g1) was upregulated during this process (Table S2). The Arabidopsis CBF genes are transcribed within a short period of exposure to cold stress and subsequently induce the COR gene expression [41, 42]. We detected two COR genes that were deregulated in the T. chinensis seeds under cold stress: TR39356|c0_g1 (cold-regulated 413 plasma membrane protein 2-like) and TR38199|c0_g1 (cold-regulated gibberellin-regulated protein 1 LTCOR12) (Table S3). In Arabidopsis, CBFs were shown to be negatively regulated by a cold-induced C2H2 zinc finger transcription factor gene, ZAT12 [43], which is also regulated by ICE1 [44], and to be inducers of the C2H2 transcription factor gene ZAT10 [45]. In the present study, we observed the upregulation of ZAT10 and the downregulation of ZAT12 (Table S3) in T. chinensis seeds under cold stress. CBF2, ZAT12, and ZAT10 were shown to regulate 172, 67, and 54 COR genes in Arabidopsis [46]. Thus, the upregulation of CFB and COR genes might be regulated by other modulators, such as MYB15, a negative regulator of CBFs in Arabidopsis [47], and GATA TFs, which were observed to be downregulated in T. chinensis seeds under cold stress (Figure 4(b)).LRR-RLK (LRR receptor-like serine/threonine-protein kinase) might be another group of PKs participating in the early response to cold stimuli inT. chinensis seeds, as we identified 8 LRR-RLK genes (5, 1, and 5 in A1, A2, and A3, respectively) that were differentially expressed in the T. chinensis seeds under cold stress (Figure 4(c) and Table S3). LRR-RLKs are a well-known class of RLKs that play important roles in plant growth, development, hormone perception, and responses to biotic/abiotic stresses [48]. They have been reported to be positive regulators of cold tolerance in Glycine soja [48], Glycine max [48], and Oryza sativa [49]. We also observed 10 aquaporin genes that were upregulated in the T. chinensis seeds under cold stress (Table 2 and Table S3). Low environmental temperature has been proven to inhibit water uptake by roots [50]. In chilling-sensitive rice, low temperature treatment resulted in a gradual increase in the expression of aquaporin genes [51]. Li et al. found that the aquaporin gene GhTIP1 was substantially upregulated in cotyledons but downregulated in roots within a few hours after cotton seedlings were exposed to cold [52]. These studies, including the present study, support the hypothesis that aquaporin genes might be involved in the response to cold stress in plants. RP genes might be another class of cold-induced genes in plants that potentially enhances the cold tolerance, as both large and small RP subunits have been shown to be induced in soybean in response to cold [53]. In tobacco, ribosomal proteins Rps2, Rps4, and Rpl20 are essential for cell survival, and Rpl33 is required for sustaining a sufficient plastid translation capacity in cold temperatures [54]. The upregulation of RP genes was also observed in Hippophae rhamnoides under cold and freeze stress [55]. In T. chinensis seeds, we observed the upregulation of RP genes (Table 2 and Table S3), including RPL36, RPS3, RPS4, and RPS8. However, additional experiments are needed to elucidate the regulatory network of RP genes and their functions in T. chinensis seeds under cold stress, along with some other genes, such as CDC20, YLS9, EXORDIUM, AUX1, and TR43917|c0_g1 (wound-responsive family protein).


## 5. Conclusions

In conclusion, we investigated the transcriptome ofT. chinensis seeds in response to cold stress. The MAPK pathway might participate in the early response to cold, and the upregulation of CBF genes might be mediated by other regulators, such as MYB and GATA TFs, rather than ICE1. The deregulation of the RP genes, CDC20, YLS9, EXORDIUM, and AUX1, might be a novel mechanisms activated in T. chinensis seeds in response to cold. This transcriptome study is the first to analyze T. chinensis seeds under cold stress. The results will improve our understanding of mechanisms regulating gene expression in plants under cold stress. More importantly, the results from this study will provide a valuable resource for future studies of T. chinensis and benefit the breeding program of T. chinensis.



---

*Source: 9247169-2022-07-06.xml*