Online ISSN 2092-9293 Print ISSN 1976-9571

Genes & Genomics https://doi.org/10.1007/s13258-018-0681-5

RESEARCH ARTICLE

Transcriptome profiling of rubber tree (Hevea brasiliensis) discovers candidate regulators of the cold stress response Xiao‑Xiao Gong1 · Bing‑Yu Yan1 · Jin Hu1 · Cui‑Ping Yang1 · Yi‑Jian Li2 · Jin‑Ping Liu1 · Wen‑Bin Liao3 Received: 27 October 2017 / Accepted: 28 February 2018 © The Genetics Society of Korea and Springer Science+Business Media B.V., part of Springer Nature 2018

Abstract Tropical plant rubber tree (Hevea brasiliensis) is the sole source of commercial natural rubber and low-temperature stress is the most important limiting factor for its cultivation. To characterize the gene expression profiles of H. brasiliensis under the cold stress and discover the key cold stress-induced genes. Three cDNA libraries, CT (control), LT2 (cold treatment at 4 °C for 2 h) and LT24 (cold treatment at 4 °C for 24 h) were constructed for RNA sequencing (RNA-Seq) and gene expression profiling. Quantitative real time PCR (qRT-PCR) was conducted to validate the RNA-Seq and gene differentially expression results. A total of 1457 and 2328 differentially expressed genes (DEGs) in LT2 and LT24 compared with CT were respectively detected. Most significantly enriched KEGG pathways included flavonoid biosynthesis, phenylpropanoid biosynthesis, plant hormone signal transduction, cutin, suberine and wax biosynthesis, Pentose and glucuronate interconversions, phenylalanine metabolism and starch and sucrose metabolism. A total of 239 transcription factors (TFs) were differentially expressed following 2 h or/and 24 h of cold treatment. Cold-response transcription factor families included ARR-B, B3, BES1, bHLH, C2H, CO-like, Dof, ERF, FAR1, G2-like, GRAS, GRF, HD-ZIP, HSF, LBD, MIKC-MADS, M-type MADS, MYB, MYB-related, NAC, RAV, SRS, TALE, TCP, Trihelix, WOX, WRKY, YABBY and ZF-HD. The genome-wide transcriptional response of rubber tree to the cold treatments were determined and a large number of DEGs were characterized including 239 transcription factors, providing important clues for further elucidation of the mechanisms of cold stress responses in rubber tree. Keywords  Rubber tree · Hevea brasiliensis · Cold stress response · RNA sequencing · Transcriptome · Transcription factor

Introduction Xiao-Xiao Gong, Bing-Yu Yan and Jin Hu have equally contributed to this work. Electronic supplementary material  The online version of this article (https​://doi.org/10.1007/s1325​8-018-0681-5) contains supplementary material, which is available to authorized users. * Jin‑Ping Liu [email protected] * Wen‑Bin Liao [email protected] 1



Hainan Key Laboratory for Sustainable Utilization of Tropical Bioresources, Tropical Agriculture and Forestry Institute, Hainan University, Haikou 570228, Hainan, China

2



Service Center of Science and Technology, Rubber Research Institute, Chinese Academy of Tropical Agricultural Sciences, Danzhou 571737, Hainan, China

3

Institute of Tropical Bioscience and Biotechnology, Chinese Academy of Tropical Agricultural Sciences, Haikou 571101, Hainan, China



Rubber tree (Hevea brasiliensis) is the almost exclusive source of commercial natural rubber. Due to its outstanding properties such as high elasticity, resilience and toughness which make it superior to synthetic rubber, natural rubber is required for the production of more than 50,000 products including tires, footwear and medical devices (Priyadarshan 2011; Gronover et al. 2011). Rubber tree is indigenous to the tropical regions of South America (Amazon Basin), and recently domesticated and transplanted to South–East Asia. Presently, the main rubber-producing countries include Thailand, Indonesia, Malaysia, India and China. The increased global demand for natural rubber results in expansion of rubber tree cultivation into the non-traditional rubbergrowing zones with suboptimal climate conditions (Pushparajah 1983; Priyadarshan 2011). Of all stress situations, low-temperature stress in winter is the most important limiting factor for rubber tree cultivation (Huang and Pan 1992;

13

Vol.:(0123456789)



Priyadarshan 2011; Gronover et al. 2011). Under chilling stress, Hevea plants experience a series of physiological and biochemical dysfunctions such as changes in membrane permeability and impairments of membrane systems, damages of photosynthetic apparatus and reduction of photosynthetic capacity, increase in cell osmotic pressure and cell dehydration, disintegration of membranous organelles, etc. (He et al. 1982; Lin and Yang 1994; Mai et al. 2009; Tian et al. 2016; Wang et al. 1978, 2012). Chilling injuries displayed by fieldgrown Hevea trees include defoliation, dieback, bark splitting and bleeding, and death of trunk and branches (Chen et al. 2012; Kan et al. 2009; Luo et al. 2012). For most temperate plants, low-temperature stress triggers a highly complex reprogramming of the transcriptome and subsequently leads to physiochemical modifications and structures remodeling which enhance their cold tolerance (i.e. cold acclimation) (Chinnusamy et al. 2007; Theocharis et  al. 2012; Thomashow 1999). For example, 4–20% of the Arabidopsis genome showed changes in gene expression patterns in response to low temperature (Chinnusamy et al. 2010; Hannah et al. 2005; Lee et al. 2005; Provart et al. 2003). A wide variety of cold responsive (COR) genes have been identified from cold-acclimated plants and these COR genes encode proteins linked with the detoxification of reactive oxygen species (ROS), the biosynthesis of carbohydrates and lipids and secondary metabolites, and the accumulation of osmolytes, cryoprotectants, transporters, chaperones, transporters, heat-shock proteins (HSPs), late embryogenesis abundant (LEA) proteins (Van Buskirk and Thomashow 2006; Chinnusamy et al. 2007; Fowler and Thomashow 2002; Knight and Knight 2012; Theocharis et al. 2012; Thomashow 1999; Zhao et al. 2015). Presently, the best studied signal transduction and regulatory pathways for cold response is dehydration-responsive element-binding protein 1 s (DREB1s)/C-repeat binding factors (CBFs) pathway. After cold induction, DREB1/CBF, a class of APETALA2/ETHYLENE RESPONSE FACTOR family transcription factors, bind to cis-elements in the promoters of COR genes and activate their expression (Stockinger et al. 1997; Liu et al. 1998; Gilmour et al. 2000, 2004). Inducer of CBF EXPRESSION1 (ICE1), a MYC-type transcription factor, acts as a master upstream regulator and controls DREB1/ CBF pathway (Thomashow 1999; Chinnusamy et al. 2007). Moreover, hormones play an important signaling role in cold stress and control the CBF-dependent and -independent pathway (Eremina et al. 2016; Sharma and Laxmi 2016). The above-mentioned mechanisms of cold acclimation are largely derived from the studies of the freeze-tolerant, model plant Arabidopsis thaliana and other temperate plants. For the tropical plants such as rubber tree which has little tolerance to low temperatures, the molecular mechanisms underlying the cold stress response have been scarcely investigated. Specifically, a fundamental question is whether

13

Genes & Genomics

the key transcription factors ICE-DREB1/CBF members which contribute significantly to the cold acclimation and freezing tolerance in many temperate plants also play critical roles in the regulation of cold stress in the nonacclimated, tropical plant rubber tree. RNA-sequencing (RNA-Seq), a powerful molecular technology and ‘omics’ tool, enables the gene expression study at whole transcriptome level and the identification of genes involved in specific biological processes in both model and non-model plants (Strickler et al. 2012; Varshney et al. 2009). With the introduction of next-generation RNA-seq, it has been extensively used to characterize gene expression patterns in diverse plant processes (Ding et al. 2017; Hua et al. 2017; Li et al. 2017; Wang et al. 2017a, b). In this study, RNA-seq and transcriptome analysis were used to characterize the gene expression profiles of H. brasiliensis under the cold stress to increase the understanding of transcriptional regulation of cold stress response and discover the novel cold stress-induced genes which might be important for rubber tree genetic improvement of cold tolerance through genetic engineering. Transcription factors (TFs) generally act as master regulators of expressions of many cold stress-responsive genes, so they are good candidates for genetic engineering to develop cold-tolerant rubber tree. Therefore, in current study, more emphasis was put on the identification of cold-responsive TFs.

Materials and methods Plant materials and low temperature treatment One-year old bud-grafting plants of rubber tree clone Reyan 7-33-97, a high-yielding and cold-resistant rubber tree clone, were used for low temperature treatment. The plants were incubated in a growth chamber under weak light (cool-white fluorescent light at 35 µmol m−2 s−1) and a 16/8-h photoperiod at 4 °C and 75% relative humidity for cold treatments. More than three were harvested after for 0 h (CT), 2 h (LT2) and 24 h (LT24) of treatment, respectively. The samples were cut into 1–2 cm and immediately frozen in liquid nitrogen and stored at − 80 °C prior to transporting to Beijing CapitalBio Corporation (Beijing, China) for RNASeq analysis.

RNA‑Seq library construction and sequencing Total RNAs were extracted from a pool of more than three biological replicates of the stems using RNeasy Plant Mini Kit (Qiagen) according to the manufacturer’s instructions. RNA was checked for quality with Agilent 2100 BioAnalyzer (Agilent) and quantified using QUBIT RNA ASSAY KIT (Invitrogen). For mRNA library construction, at least

Genes & Genomics

20 mg of total RNA samples were prepared using the NEBNext Ultra RNA Library Prep Kit for Illumina (NEB) and Multiplex Oligos (NEB). Subsequently, RNA sequencing was performed on an Illumina HiSeq 2500 platform.

Data processing and DEG analysis Raw reads were cleaned by removing reads containing an adaptor, low-quality reads with quality scores less than Q30 and reads in which unknown nucleotides accounted for more than 5%. The Cufflinks pipeline (Roberts et al. 2011) was used to assemble transcripts individually for each sample group and to estimate their abundance. Afterwards we employed Cuffmerge (Trapnell et al. 2012a, b) to produce a union of transcripts from all three sample groups. The resultant transcripts were annotated with the reference genome sequence information (Tang et al. 2016) and filtered to discard the transcripts in which the total length of exons is lower than 100 bp. All clean reads were deposited in the NCBI Sequence Read Archive (SRA, http://www.ncbi. nlm.nih.gov/Trace​s/sra) database with accession number SRP119284. Cuffdiff (Trapnell et al. 2012a, b) was used to estimate gene expression levels for each group and analyze the DEGs between samples with the ‘blind’ dispersion method (for the samples with no biological replicates). Transcript level was calculated with the relative abundance of the transcripts and was given by FPKM (Fragments Per Kilo bases per Million fragments). The following formula was used: FPKM = Total exon fragments/mapped reads (Millions) × exon length (Kb). The default criteria to identify DEGs was |log2 (fold change)| ≧ 1 and p-value ≦ 0.05. Gene differential expression analysis was based on the independent statistical hypothesis test of expression levels of thousands of genes. The resultant p-value was corrected with the false discovery rate (FDR). For gene expression profiling, Gene Ontology (GO) enrichment analysis was conducted using agriGO (http:// bioin​fo.cau.edu.cn/agriG​O/) (Du et  al. 2010). Pathway enrichment analysis was performed using KOBAS 2.0 (http://kobas.​ cbi.pku.edu.cn/) (Xie et al. 2011), which incorporated 5 pathway databases, i.e. Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa and Goto 2000), Pathway Interaction Database (PID) (Schaefer et al. 2009),

BioCyc (Karp et al. 2005), Reactome (Matthews et al. 2009; Croft et al. 2011) and PANTHER (Thomas et al. 2003).

Quantitative real‑time PCR analysis Quantitative real time PCR (qRT-PCR) was performed to validate the RNA-Seq and DGE results as described previously (Liu et al. 2015, 2016). Twelve DEGs were randomly selected and Eukaryotic Initiation Factor 2 (eIF2) was used as the internal control. All of the samples were measured in triplicate. The relative expressions of the genes were calculated using the 2­ −△△Ct method. All of the primers that were used in this study are listed in Additional file 1.

Results and discussion RNA sequencing and transcriptome assembly The Illumina HiSeq 2500 platform was used to sequence the three cDNA libraries, i.e. CT (4 °C treatment for 0 h), LT2 (4 °C treatment for 2 h) and LT24 (4 °C treatment for 24 h). We obtained 90,662,140, 84,883,838 and 81,138,988 raw reads from CT, LT2 and LT24, respectively (Table 1). After filtering, 88,006,896, 82,452,710 and 78,805,348 clean reads were respectively produced from CT, LT2 and LT24, with the average Q30 quality scores of the three transcriptomes above 96% and the percentage of mapped reads over 88%. These high-quality reads were then assembled to yield 39,133, 39,095 and 42,684 non-redundant genes in CT, LT2 and LT24, respectively (Table 2).

Identification of differentially expressed genes (DEGs) in response to cold stress A total of 539 up-regulated genes (Additional file 2) and 918 down-regulated genes (Additional file 3) were detected between LT2 and CT libraries, and 1436 up-regulated genes (Additional file 4) and 892 down-regulated genes (Additional file 5) were detected between LT24 and CT libraries (Fig. 1a). It is noted that the number of up-regulated genes increased with longer cold treatment but the number of down-regulated genes largely remained constant between two time points. As shown in venn diagram (Fig. 1b), 362

Table 1  Overview of transcriptome sequencing and mapping Sample Raw reads number Clean reads number (%) CT LT2 LT24

90,662,140 84,883,838 81,138,988

88,006,896 (97.07) 82,452,710 (97.14) 78,805,348 (97.12)

Clean bases

Q30 (%) Mapped reads (%)

13,160,716,238 96.01 12,333,970,211 96.04 11,776,107,225 96.07

Uniquely mapped reads (%)

78,146,104 (88.80) 72,497,490 (92.77) 72,812,177 (88.31) 67,158,452 (92.24) 69,576,584 (88.29) 64,229,483 (92.31)

Multiple mapped reads (%) 5,648,614 (7.23) 5,653,725 (7.76) 5,347,101 (7.69)

13



Genes & Genomics

Table 2  Statistics of transcriptome assembly results Sample

Gene number

Transcript number

Exon total length (bp)

Average transcript length (bp)

Max transcript Min transcript N50 length (bp) length (bp) length (bp)

CT LT2 LT24

39,133 39,095 42,684

58,155 58,072 63,927

99,333,262 99,663,338 127,589,761

1708 1716 1996

16,630 17,307 23,735

Fig. 1  The statistics of DEGs obtained in this study. a Numbers of DEGs between LT2 and CT (a), between LT24 and CT, and between LT24 and LT2. b Venn diagram of the numbers of up-regulated DEGs between LT2 and CT, between LT24 and CT, and between LT24 and LT2. c Venn diagram of the numbers of down-regulated DEGs between LT2 and CT, between LT24 and CT, and between LT24 and LT2

and 1259 up-regulated DEGs were specifically detected in LT2 and LT24, and 177 DEGs were up-regulated both in LT2 and LT24 compared with CT. In addition, 529 and 503 down-regulated DEGs were uniquely detected in LT2 and LT24, and 389 DEGs were down-regulated both in LT2 and LT24 compared with CT (Fig. 1c).

Gene ontology classification and KEGG pathway analysis of cold‑responsive genes The Gene ontology (GO) was used to classify the functions of all DEGs according to their biological process, cellular

13

150 150 150

2109 2129 2573

component and molecular function (Fig. 2). In the GO terms of “biological process”, the categories of “metabolic process” and “cell process” were prominent. The categories of “membrane part” and “cell” were dominant in the GO terms of “cellular component”. In the GO terms of “molecular function”, a high percentage of DEGs were associated with “catalytic activity” and “binding”. To understand the biological function of DEGs, we mapped the cold-responsive DEGs to terms in the Kyoto Encyclopedia of Genes and Genomes (KEGG) database, and 183 and 215 significantly enriched metabolic pathways or signal transduction pathways were identified in DEGs in LT2 and LT24 compared with CT, respectively. Top 20 significantly enriched pathways after 2- and 24-h cold treatments were shown in Fig. 3. The results showed that, of the top 20 significantly enriched pathways, seven pathways including Flavonoid biosynthesis, Phenylpropanoid biosynthesis, Plant hormone signal transduction, Cutin, suberine and wax biosynthesis, Pentose and glucuronate interconversions, Phenylalanine metabolism and Starch and sucrose metabolism, were significantly enriched in both LT2 and LT24 samples (Additional file 6). Terpenoid backbone biosynthesis, Galactose metabolism, Brassinosteroid biosynthesis, Diterpenoid biosynthesis, Protein digestion and absorption, Calcium signaling pathway, Stilbenoid, diarylheptanoid and gingerol biosynthesis, Circadian rhythm—plant, Steroid biosynthesis, Carotenoid biosynthesis, Butanoate metabolism, Cysteine and methionine metabolism, and Steroid hormone biosynthesis were specifically enriched in LT2 sample, and Plant–pathogen interaction, Flavone and flavonol biosynthesis, Fatty acid elongation, Leishmaniasis, Pertussis, Apoptosis, NF-kappa B signaling pathway, Toxoplasmosis, Tolllike receptor signaling pathway, Measles and Chagas disease (American trypanosomiasis) specifically enriched in LT24 sample, suggesting that different pathways are involved in the cold stress response of rubber tree at different time points. These KEGG annotations provide a valuable resource for investigating the pathways and processes involved in the early-stage and late-stage cold stress response. Flavonoids are a group of plant polyphenolic secondary metabolites including anthocyanins, flavonols, flavanols and proanthocyanidins (PAs) or condensed tannins (Petrussa et al. 2013). Cold stress has been shown to induce

Genes & Genomics

Fig. 2  GO classification of the annotated DEGs between LT2 and CT (a) and between LT24 and CT (b). The biological process, the cellular component and the molecular function are indicated by blue, green and red, respectively

anthocyanin synthesis in diverse species (Chalker-Scott 1999; Christie et al. 1994) and flavonoid biosynthesis-related genes in grape skin are up-regulated by low temperature combined with the presence of light (Azuma et al. 2012).

Ahmed et al. (2015) and Schulz et al. (2016) provided evidence that flavonoids are determinants of freezing tolerance and cold acclimation in Arabidopsis thaliana and Brassica rapa. It is well known that cold stress increases plant

13



Genes & Genomics

Fig. 3  Top 20 significantly enriched pathways between LT2 and CT (a) and between LT24 and CT (b). Five pathway databases (KEGG, PID, BioCyc, Reactome and PANTHER) are indicated by five different colors

phenylpropanoid phenolic production (Solecka 1997). The majority of the genes in the phenylpropanoid pathway were up-regulated by cold stress in Lotus japonicus (Calzadilla et al. 2016). Furthermore, “Phenylpropanoid biosynthesis” was the most highly DEG-enriched KEGG pathways under cold stress (Kaplan et al. 2007; Khaledian et al. 2015; Zhan et al. 2016). However, majority of genes in these two processes were down-regulated, contrary to the induction of these two processes in the cold acclimation of freezingtolerant plants. This might partially explain why the rubber tree is freezing-intolerant. In addition, groups of DEGs involved in different plant hormone signal transduction pathways, such as abscisic acid (ABA), JA (Jasmonic acid), auxin, ethylene (ET), cytokinin

13

(CK), gibberellic acid (GA), brassinosteroid (BR) and salicylic acid (SA) pathways, were identified in rubber tree subjected to cold stress (Additional file 6). For example, in ABA signaling, Protein Phosphatase 2C (PP2C) was shown to be down-regulated but ABA receptor PYRABACT IN RESISTANCE (PYR1)/PYR1-LIKE (PYL) up-regulated both in LT2 and LT24 compared with CT. The genes encoding the components of ET signaling such as ethylene receptor, EIN3binding F-box protein and ethylene-responsive transcription factor were all down-regulated in LT2 and LT24. It is not surprising since cold stress responses in plants are highly sophisticated events and complex phytohormone signaling cascades are essentially utilized to induce changes in coldresponsive gene expression (Eremina et al. 2016).

Genes & Genomics

Cold‑responsive TFs in rubber tree A total of 239 TFs detected by our RNA sequencing and transcriptome analysis were differentially expressed following 2 h or/and 24 h of cold treatment (Table 3, Additional file 7). TF families showing especially strong responses to cold treatment included ARR-B, B3, BES1, bHLH, C2H, CO-like, Dof, ERF, FAR1, G2-like, GRAS, GRF, HD-ZIP, HSF, LBD, MIKC-MADS, M-type MADS, MYB, MYBrelated, NAC, RAV, SRS, TALE, TCP, Trihelix, WOX, WRKY, YABBY and ZF-HD. Many of these families of TFs were reported to play central roles in plant response to abiotic stress including cold stress (Lindemose et al. 2013) and to be utilized to improve plant abiotic stress tolerance by gene transfer technology (Wang et al. 2016). In particular, of differentially expressed bHLH TFs, MYC4, bHLH13, bHLH14 and bHLH111 exhibited a coldinducible response and other bHLH DNA-binding family protein genes showed a cold-repressed expression. MYC4 is one of the important JAZ-interacting transcription factors that regulate jasmonate (JA) responses (Fernández-Calvo et al. 2011; Niu et al. 2011). Song et al. (2013) identified the bHLH subgroup IIId transcription factors such as bHLH3, bHLH13, bHLH14 and bHLH17 as novel targets of JAZs which act as transcription repressors to negatively regulate JA responses. JA is a lipid-derived phytohormone that essential for plant defense and recent studies showed evidence supporting the involvement of JA in tolerance to cold stress (Hu et al. 2017). Of C2H2-type zinc fingers displayed differential expression, ZAT12 was notably up-regulated both in LT2 and LT24, and ZAT10, ZAT11, Zinc finger CCCH domain-containing protein 29, 39 and 48 were up-regulated in LT24, compared with C. In Arabidopsis, ZAT12 and ZAT10 were rapidly induced in response to low temperature in parallel with CBF genes, and they function together with other regulons to form a complex low temperature regulatory network (Vogel et al. 2005; Park et al. 2015). In addition, ZAT10 and ZAT11 were also implicated in Dendrobium officinale cold acclimation (Wu et al. 2016). Furthermore, several zinc-finger transcription factors were discovered to be responsive to cold treatment and have been identified as novel cold-stress regulators in A. thaliana and rice (Brown 2005; Figueiredo et al. 2012; Van Buskirk and Thomashow 2006). The APETALA 2/ethylene-responsive element binding factor (AP2/ERF) family includes four major subfamilies: the AP2, RAV, ERF and dehydration-responsive elementbinding protein (DREB) subfamilies, and the AP2/ERF superfamily plays an important regulatory role in the plant responses to various stresses (Mizoi et al. 2012). In particular, CBF transcription factors such as DREB1A/CBF3, DREB1B/CBF1 and DREB1C/CBF2 were all induced by cold and the DREB1/CBF regulon appears to be conserved

for cold acclimation in many temperate plants including Arabidopsis. In this study, ERF109, ERF061, ERF098 and an unidentified ERF showed high expression under cold treatment at the time points of both 2 and 24 h, RAP2-11, ERF114, ERF012, DREB1D, DREB2F. DREB2A and DREB2C demonstrated up-regulation at the time points of 24 h and no DREB1As was detected to be up-regulated. ERF109 has been reported to retard PCD and improve salt tolerance in Arabidopsis (Bahieldin et al. 2016) and also showed strong response to cold treatment in Nicotiana tabacum (Jin et  al. 2017). Dehydration-responsive element-binding protein 1D was identified to be up-regulated in both banana and plantain after 3 h of cold stress (Yang et al. 2015). However, the members of DREB1/CBF regulon hadn’t been detected as DEGs in the tropical plant cassava which belongs to the same family (Euphorbiaceae) (An et al. 2012; Zeng et al. 2014). Although the DREB1/CBF regulon plays a pivotal role in gene regulation and appears to be one of the main regulatory pathways during cold acclimation, other regulatory pathways still exist (Fowler and Thomashow 2002). For example, in Arabidopsis, co-regulation of COR genes is not limited to members of the DREB1/ CBF regulon and 11 other TFs were rapidly cold-induced in parallel with the “first-wave” DREB1/CBF genes (Park et al. 2015). The DREB1/CBFs regulon only accounts for a small percentage of the COR genes and about two-thirds of the COR genes are not regulated by it (Fowler and Thomashow 2002; Zhao et al. 2015). Probably the DREB1/CBF is not the main regulatory pathway in the cold-stress response in rubber tree and cassava. On the other hand, possibly because it is well-established that DREB1/CBF regulon plays a central role in cold acclimation, tropical plants such as rubber tree and cassava are therefore freezing-intolerant. MYB transcription factors has been proved to be active players in abiotic stress signaling (Roy 2016). In our study, Members of the MYB transcription factor family were also differentially regulated by cold stress. Specially, Myb-related protein Myb4, WER, ODORANT1, GAMYB, Myb-related protein 308, MYB44 were induced by cold stress. MYB4 regulation of cold stress response in plants has been shown by the enhancement of freezing tolerance of Arabidopsis through the overexpression of the rice MYB4 (Vannini et al. 2004). MYB44 which has possible implication in redox balancing under stress conditions, has been characterized as phosphorylation-dependent positive regulator of salt stress signaling (Chinnusamy et al. 2010; Persak and Pitzschke 2014; Seo et al. 2012) and plays multiple roles in ABA signaling, stress responses and leaf senescence (Jaradat et al. 2013). Recently, MYB44 was found to be rapidly up-regulated by low temperature (Park et al. 2015) and it is also possibly a cold stress transcriptional regulator. A majority of WRKY gene family members displayed up-regulation at the time point of 24 h, including WRKY6,

13



Genes & Genomics

Table 3  DEGs of TF in LT2 and L24 compared with CT Tracking ID

AGI locus

Family

Description

log2 (fold change) LT2 vs. CT LT24 vs. CT

XLOC_043671 XLOC_005214 XLOC_039496 XLOC_009019 XLOC_000485 XLOC_007831 XLOC_007911 XLOC_008902 XLOC_009977 XLOC_010379 XLOC_012038 XLOC_013133 XLOC_013993 XLOC_018133 XLOC_019830 XLOC_022935 XLOC_024385 XLOC_025634 XLOC_026602 XLOC_027543 XLOC_032816 XLOC_032963 XLOC_034890 XLOC_043311 XLOC_043613 XLOC_044760 XLOC_045471 XLOC_005673 XLOC_009840 XLOC_011924 XLOC_015858 XLOC_022283 XLOC_037288 XLOC_045875 XLOC_046168 XLOC_047062 XLOC_001178 XLOC_002891 XLOC_005543 XLOC_011278 XLOC_017092 XLOC_026769 XLOC_027140 XLOC_028823 XLOC_031590 XLOC_032275 XLOC_033311 XLOC_038632

13

AT2G25180 AT3G19184 AT3G19184 AT1G78700 AT5G65640 AT4G17880 AT5G41315 AT1G72210 AT1G22490 AT5G41315 AT5G41315 AT1G25330 AT1G73830 AT3G24140 AT1G22490 AT1G01260 AT5G50915 AT4G00870 AT4G20970 AT4G20970 AT5G50915 AT1G73830 AT4G17880 AT4G17880 AT1G31050 AT2G42280 AT5G65640 AT5G11260 AT3G30530 AT1G08320 AT3G49760 AT3G56850 AT3G49760 AT3G58120 AT4G35900 AT2G40620 AT4G35610 AT1G80730 AT4G06634 AT1G26610 AT4G35280 AT5G10970 AT1G10480 AT5G25160 AT5G59820 AT1G66140 AT1G27730 AT1G80730

ARR-B B3 B3 BES1 bHLH bHLH bHLH bHLH bHLH bHLH bHLH bHLH bHLH bHLH bHLH bHLH bHLH bHLH bHLH bHLH bHLH bHLH bHLH bHLH bHLH bHLH bHLH bZIP bZIP bZIP bZIP bZIP bZIP bZIP bZIP bZIP C2H2 C2H2 C2H2 C2H2 C2H2 C2H2 C2H2 C2H2 C2H2 C2H2 C2H2 C2H2

Two-component response regulator ARR18 B3 domain-containing protein At3g19184 B3 domain-containing protein Os01g0234100 BES1/BZR1 homolog protein 4 Transcription factor bHLH93 Transcription factor MYC4 Transcription factor GLABRA 3 Transcription factor bHLH96 Transcription factor bHLH94 Transcription factor GLABRA 3 Transcription factor GLABRA 3 Transcription factor bHLH75 Transcription factor BEE 3 Transcription factor FAMA Transcription factor bHLH94 Transcription factor bHLH13 Transcription factor bHLH137 Transcription factor bHLH14 Transcription factor bHLH36 – Transcription factor bHLH137 Transcription factor BEE 3 Transcription factor MYC4 Transcription factor MYC4 Transcription factor bHLH111 Transcription factor bHLH130 Transcription factor bHLH93 Transcription factor HY5 Ocs element-binding factor 1 TGACG-sequence-specific DNA-binding protein TGA-2.1 bZIP transcription factor TRAB1 ABSCISIC ACID-INSENSITIVE 5-like protein 2 bZIP transcription factor TRAB1 Transcription factor VIP1 Protein FD Transcription factor RF2b Zinc finger protein ZAT3 Zinc finger protein 1 – – – Zinc finger protein 1 – – Zinc finger protein ZAT12 Zinc finger protein 4 Zinc finger protein ZAT10 Zinc finger protein 1

− 1.41164 − 1.59891 2.47607 − 1.25528 − 1.71893 − 2.13039

− 1.94686 − 1.78755 − 1.61862 2.28982 − 1.45409 5.89944 − 1.62784 − 3.23034 2.88505 1.98701 1.00578 − 1.17887 − 1.85783 1.36057 − 2.57835 − 2.459 − 8.24254 − 1.72374

10.14625 6.38309 − 2.22394

1.62964 − 1.37639 − 1.39057

2.58574 − 9.44487 11.14090 − 2.41315 2.27344 − 2.20449 − 2.74305 − 1.597 − 1.5722 − 1.86036 − 5.71194 11.54219 − 2.64301 2.44175 2.20542 12.63442 − 4.38027 − 3.18728

− 1.83685 − 2.20383

− 2.3809

9.98031 2.90793 10.64626 3.44065 9.90357

10.08131 7.24239 3.03866 1.90288

Genes & Genomics Table 3  (continued) Tracking ID

AGI locus

Family

Description

log2 (fold change) LT2 vs. CT LT24 vs. CT

XLOC_040708 XLOC_042748 XLOC_043490 XLOC_047070 XLOC_048731 XLOC_008322 XLOC_009094 XLOC_024021

AT2G28710 AT2G37430 AT3G57670 AT1G27730 AT5G10970 AT1G03790 AT2G40140 AT4G25440

C2H2 C2H2 C2H2 C2H2 C2H2 C3H C3H C3H

XLOC_030957 XLOC_009193 XLOC_016809 XLOC_016552 XLOC_017972 XLOC_019692 XLOC_025714 XLOC_027188 XLOC_031230 XLOC_043888 XLOC_044590 XLOC_001499 XLOC_002618 XLOC_006330 XLOC_006451 XLOC_006452 XLOC_007298 XLOC_009702 XLOC_011234 XLOC_013335 XLOC_013353 XLOC_013428 XLOC_019369 XLOC_021334 XLOC_022039 XLOC_022103 XLOC_026166 XLOC_026489 XLOC_026588 XLOC_028201 XLOC_028547 XLOC_028633 XLOC_028859 XLOC_029572 XLOC_029573 XLOC_031218 XLOC_031823 XLOC_031824 XLOC_032132

AT3G19360 AT5G24930 AT3G02380 AT3G50410 AT1G21340 AT5G65590 AT3G47500 AT3G61850 AT5G66940 AT5G39660 AT2G28510 AT5G19790 AT1G75490 AT5G64750 AT5G52020 AT5G51990 AT3G23230 AT2G44940 AT3G16770 AT3G23230 AT3G23240 AT4G17500 AT4G34410 AT3G57600 AT4G17500 AT5G51190 AT3G23240 AT3G16770 AT1G28360 AT3G16770 AT5G13330 AT2G40340 AT2G40340 AT4G17490 AT4G17490 AT5G11590 AT4G17490 AT4G17490 AT1G64380

C3H CO-like CO-like Dof Dof Dof Dof Dof Dof Dof Dof ERF ERF ERF ERF ERF ERF ERF ERF ERF ERF ERF ERF ERF ERF ERF ERF ERF ERF ERF ERF ERF ERF ERF ERF ERF ERF ERF ERF

Zinc finger protein ZAT12 Zinc finger protein ZAT11 Protein TRANSPARENT TESTA 1 Zinc finger protein ZAT10 – − 1.92728 Zinc finger CCCH domain-containing protein 2 − 1.84158 Zinc finger CCCH domain-containing protein 29 1.45205 Zinc finger CCCH domain-containing protein 48|Zinc finger CCCH domain-containing protein 48 Zinc finger CCCH domain-containing protein 39 Zinc finger protein CONSTANS-LIKE 4 Zinc finger protein CONSTANS-LIKE 2 1.07243 Dof zinc finger protein DOF3.4 − 12.8364 Dof zinc finger protein DOF1.2 Dof zinc finger protein DOF5.7 Cyclic dof factor 2 Dof zinc finger protein DOF3.7 Dof zinc finger protein DOF5.8 − 2.28363 Cyclic dof factor 2 Dof zinc finger protein DOF2.1 Ethylene-responsive transcription factor RAP2-11 Dehydration-responsive element-binding protein 2D Ethylene-responsive transcription factor ERF114 Ethylene-responsive transcription factor ERF026 − 8.96750 Dehydration-responsive element-binding protein 1D Ethylene-responsive transcription factor ERF098 Ethylene-responsive transcription factor ERF034 − 1.28101 Ethylene-responsive transcription factor RAP2-3 − 2.1996 Ethylene-responsive transcription factor ERF098 11.14930 Ethylene-responsive transcription factor 1B − 1.49041 Ethylene-responsive transcription factor 1A − 1.45253 Ethylene-responsive transcription factor ERF109 2.11478 Dehydration-responsive element-binding protein 2F Ethylene-responsive transcription factor 2 Ethylene-responsive transcription factor ERF105 − 1.32326 Pathogenesis-related genes transcriptional activator PTI5 − 2.29317 Ethylene-responsive transcription factor RAP2-3 − 1.24667 Ethylene-responsive transcription factor 12 − 2.80585 Ethylene-responsive transcription factor RAP2-3 − 2.72552 Ethylene-responsive transcription factor ERF113 − 1.30671 Dehydration-responsive element-binding protein 2C|– Dehydration-responsive element-binding protein 2A Ethylene-responsive transcription factor 6 − 1.73809 Ethylene-responsive transcription factor 6 − 1.36377 Dehydration-responsive element-binding protein 3 − 2.204 Ethylene-responsive transcription factor 6 − 1.52189 Ethylene-responsive transcription factor 6 − 1.80645 Ethylene-responsive transcription factor ERF061 2.55105

6.59286 10.12581 − 5.00379 2.59781

4.0375 9.64243 − 2.25334 1.9411 2.86549 2.28813 2.15965 − 1.84848 2.18564 2.79894 11.08593 − 9.59053 11.19614 − 9.10501 4.79649 9.51601 − 4.66595 10.90287 − 9.23315 4.45676 9.17298 4.39914

− 2.49386 − 2.8679 2.19574 4.58957 − 1.58798 − 2.51516 − 1.72344 3.57873

13



Genes & Genomics

Table 3  (continued) Tracking ID

AGI locus

Family

Description

log2 (fold change) LT2 vs. CT LT24 vs. CT

XLOC_034060 XLOC_036374 XLOC_039717 XLOC_039780 XLOC_040748 XLOC_042119 XLOC_043161 XLOC_045559 XLOC_048452 XLOC_003558 XLOC_034482 XLOC_014843 XLOC_040804 XLOC_011170 XLOC_020846 XLOC_024317 XLOC_041743 XLOC_047674 XLOC_018934 XLOC_031013 XLOC_038343 XLOC_040315 XLOC_040972 XLOC_040985 XLOC_007211 XLOC_007212 XLOC_007085 XLOC_014960 XLOC_021747 XLOC_022418 XLOC_039450 XLOC_047103 XLOC_048734 XLOC_003689 XLOC_027541 XLOC_033678 XLOC_016148 XLOC_020701 XLOC_025619 XLOC_027103 XLOC_032959 XLOC_033587 XLOC_036491 XLOC_037804 XLOC_042997 XLOC_015802 XLOC_019939 XLOC_034738

13

AT1G21910 AT4G11140 AT2G40340 AT1G19210 AT2G40340 AT4G34410 AT3G23240 AT1G21910 AT4G34410 AT4G38180 AT2G27110 AT2G38300 AT4G37180 AT5G49300 AT5G66320 AT3G54810 AT5G66320 AT3G54810 AT4G17230 AT2G37650 AT5G48150 AT2G29060 AT5G48150 AT5G48150 AT3G13960 AT3G13960 AT5G53980 AT4G40060 AT5G06710 AT4G36740 AT5G06710 AT3G01470 AT5G53980 AT2G26150 AT2G26150 AT5G03720 AT3G47870 AT3G11090 AT3G27650 AT2G42430 AT3G02550 AT1G31320 AT3G27650 AT5G63090 AT2G45420 AT1G77980 AT3G02310 AT3G54340

ERF ERF ERF ERF ERF ERF ERF ERF ERF FAR1 FAR1 G2-like G2-like GATA​ GATA​ GATA​ GATA​ GATA​ GRAS GRAS GRAS GRAS GRAS GRAS GRF GRF HD-ZIP HD-ZIP HD-ZIP HD-ZIP HD-ZIP HD-ZIP HD-ZIP HSF HSF HSF LBD LBD LBD LBD LBD LBD LBD LBD LBD MIKC_MADS MIKC_MADS MIKC_MADS

Ethylene-responsive transcription factor ERF012 – Dehydration-responsive element-binding protein 2A Ethylene-responsive transcription factor ERF017 Dehydration-responsive element-binding protein 2C Ethylene-responsive transcription factor ERF109 Ethylene-responsive transcription factor 1B Ethylene-responsive transcription factor ERF012 Ethylene-responsive transcription factor ERF109 – – –|Putative Myb family transcription factor At1g14600 Transcription factor BOA GATA transcription factor 16 GATA transcription factor 5 GATA transcription factor 8 GATA transcription factor 5 GATA transcription factor 8 Scarecrow-like protein 13 Scarecrow-like protein 9 Chitin-inducible gibberellin-responsive protein 1 Scarecrow-like protein 33|Scarecrow-like protein 9 Chitin-inducible gibberellin-responsive protein 1 Chitin-inducible gibberellin-responsive protein 1 Growth-regulating factor 3 Growth-regulating factor 3 Homeobox-leucine zipper protein ATHB-52 Homeobox-leucine zipper protein ATHB-16 Homeobox-leucine zipper protein HAT14 Homeobox-leucine zipper protein ATHB-40 Homeobox-leucine zipper protein HAT14 Homeobox-leucine zipper protein HAT5 Homeobox-leucine zipper protein ATHB-52 Heat shock factor protein HSF30 Heat shock factor protein HSF30 Heat stress transcription factor A-3 LOB domain-containing protein 27 LOB domain-containing protein 21 LOB domain-containing protein 25 LOB domain-containing protein 16 LOB domain-containing protein 41 LOB domain-containing protein 4 LOB domain-containing protein 25 Protein LATERAL ORGAN BOUNDARIES LOB domain-containing protein 18 MADS-box transcription factor 13 MADS-box protein CMB1 Floral homeotic protein DEFICIENS

− 1.55259 2.68175

− 2.89744 − 4.49248 3.23031 1.69714 1.68571 − 1.60357 9.54480

2.48018

2.19774

− 1.89306 4.05085 2.21009

− 1.63332 − 3.7082

5.50717 6.25185 3.48957 9.98031 − 4.77681 5.00914 9.53798 − 1.85012 10.49706 − 2.12809 2.23632 4.29774 2.26741 1.6703 2.5355 4.37019 3.34205 3.56011 3.89221 10.61593 10.01019 − 3.53708 − 2.91975 − 2.76363 − 2.13162 − 11.97498 − 1.40491 − 2.90543 2.77192 8.99401 − 6.33878 − 9.26488 − 2.3217 − 2.50597

− 1.16013 − 1.72012 − 2.69763 − 3.21509 − 11.59676 − 8.60493 1.35559 − 9.08639 − 2.33521 − 9.74642

Genes & Genomics Table 3  (continued) Tracking ID

AGI locus

Family

Description

log2 (fold change) LT2 vs. CT LT24 vs. CT

XLOC_010605 XLOC_036605 XLOC_001521 XLOC_001572 XLOC_002968 XLOC_003366 XLOC_004805 XLOC_007287 XLOC_007536 XLOC_015924 XLOC_019156 XLOC_021958 XLOC_022522 XLOC_022724 XLOC_026886 XLOC_029789 XLOC_033634 XLOC_034297 XLOC_036451 XLOC_037973 XLOC_038614 XLOC_040954 XLOC_042804 XLOC_045005 XLOC_017221 XLOC_044502 XLOC_001505 XLOC_002326 XLOC_003030 XLOC_011394 XLOC_012090 XLOC_012833 XLOC_016192 XLOC_017472 XLOC_024147 XLOC_031523 XLOC_032863 XLOC_035723 XLOC_037351 XLOC_038970 XLOC_041934 XLOC_043846 XLOC_044096 XLOC_006163

AT5G60440 AT5G60440 AT1G22640 AT5G14750 AT2G31180 AT5G62470 AT4G21440 AT1G69560 AT1G69560 AT3G13540 AT1G22640 AT5G14750 AT5G02320 AT3G28910 AT4G12350 AT5G06100 AT5G49620 AT4G38620 AT1G69560 AT4G37260 AT5G14750 AT4G37260 AT4G37260 AT1G68320 AT1G75250 AT1G70000 AT5G61430 AT4G35580 AT5G24590 AT5G13180 AT1G61110 AT1G33060 AT4G17980 AT1G65910 AT5G64530 AT5G39610 AT5G64530 AT5G08790 AT4G17980 AT3G15510 AT1G28470 AT5G14490 AT5G18270 AT1G51120

M-type_MADS M-type_MADS MYB MYB MYB MYB MYB MYB MYB MYB MYB MYB MYB MYB MYB MYB MYB MYB MYB MYB MYB MYB MYB MYB MYB_related MYB_related NAC NAC NAC NAC NAC NAC NAC NAC NAC NAC NAC NAC NAC NAC NAC NAC NAC RAV

XLOC_043106 AT1G51120 RAV XLOC_015130 AT1G02065 SBP XLOC_002721 AT1G75520 SRS

Agamous-like MADS-box protein AGL62 Agamous-like MADS-box protein AGL62 Anthocyanin regulatory C1 protein Transcription factor WER Myb-related protein Myb4 Myb-related protein 306 Transcription factor MYB39 Transcription factor MYB44 Transcription factor MYB44 Myb-related protein 308 Anthocyanin regulatory C1 protein Transcription factor WER Myb-related protein 3R-1 Myb-related protein 306 Protein ODORANT1 Transcription factor GAMYB Transcription factor MYB108 Myb-related protein 308 Transcription factor MYB44 – Transcription factor WER Transcription factor MYB44 Transcription factor MYB44 Transcription factor MYB108 Transcription factor RADIALIS –|Transcription factor MYB1R1 NAC domain-containing protein 100 Protein FEZ NAC domain-containing protein 78 NAC transcription factor 25 NAC transcription factor 25 –|NAC domain-containing protein 89 NAC domain-containing protein 78 NAC domain-containing protein 7 NAC transcription factor ONAC010 NAC transcription factor ONAC010 NAC transcription factor ONAC010 NAC domain-containing protein 2 NAC domain-containing protein 78 NAC transcription factor NAM-B2 NAC domain-containing protein 8 NAC domain-containing protein 8 NAC domain-containing protein 100 AP2/ERF and B3 domain-containing transcription factor At1g51120 AP2/ERF and B3 domain-containing transcription factor At1g51120 Squamosa promoter-binding-like protein 8 Protein SHI RELATED SEQUENCE 5

6.03757 6.47764 − 3.7006 − 1.83973

9.44656 9.70988 − 3.17361 4.00503 − 2.02362

− 1.90541 − 10.38741 − 10.52491 − 1.13735 − 1.98597 − 2.24839 − 4.47875 − 3.32976 − 3.22121 9.93219 1.7587 − 1.84964 1.61938 4.529 − 2.18203 − 2.5997 1.60543 − 3.72907 4.3922 − 3.32801 1.56022 2.60462 − 3.56744 − 10.12741 − 11.00701 − 11.13921 − 1.82957 1.93911 1.89319 3.62125 3.49181 − 2.05717 6.81750 9.31361 5.94478 − 1.57553 − 1.74713

− 2.23896 − 1.22614 − 2.33805 1.49839

− 2.69549 − 2.49908 2.80721 4.38309 2.77918 9.23875

5.00453 − 1.21469 − 2.99098

− 2.26021

13



Genes & Genomics

Table 3  (continued) Tracking ID

AGI locus

Family

Description

log2 (fold change) LT2 vs. CT LT24 vs. CT

XLOC_020869 XLOC_034719 XLOC_008511 XLOC_021077 XLOC_026996 XLOC_028829 XLOC_006460 XLOC_039807 XLOC_046509 XLOC_048562 XLOC_001978 XLOC_015312 XLOC_021823 XLOC_028245 XLOC_029006 XLOC_033665 XLOC_035380 XLOC_041389 XLOC_042216 XLOC_005996 XLOC_036149 XLOC_000967 XLOC_002774 XLOC_004423 XLOC_007755 XLOC_009539 XLOC_011418 XLOC_014526 XLOC_020024 XLOC_021866 XLOC_022921

AT3G51060 AT3G51060 AT1G62360 AT1G62360 AT1G62360 AT4G36870 AT2G45680 AT1G68800 AT4G18390 AT1G68800 AT3G10040 AT3G24860 AT1G21200 AT3G54390 AT5G03680 AT5G03680 AT5G01380 AT1G76880 AT5G03680 AT1G20700 AT3G11260 AT4G31550 AT2G37260 AT1G29280 AT1G62300 AT3G56400 AT1G62300 AT5G64810 AT3G56400 AT2G47260 AT5G24110

SRS SRS TALE TALE TALE TALE TCP TCP TCP TCP Trihelix Trihelix Trihelix Trihelix Trihelix Trihelix Trihelix Trihelix Trihelix WOX WOX WRKY WRKY WRKY WRKY WRKY WRKY WRKY WRKY WRKY WRKY

XLOC_023998 AT1G29860 WRKY XLOC_026615 AT2G38470 WRKY XLOC_027736 AT5G64810 WRKY

Protein SHI RELATED SEQUENCE 1 − 2.83017 Protein SHI RELATED SEQUENCE 1 − 2.52588 Homeobox protein SBH1 1.59295 Homeobox protein SBH1 1.17173 Homeobox protein SBH1 1.0014 BEL1-like homeodomain protein 2 1.54988 Transcription factor TCP9 1.84253 Transcription factor TCP12 − 9.50904 Transcription factor TCP2 1.01143 Transcription factor TCP12 − 13.94380 − 8.98038 – − 1.61719 − 3.35225 – 2.07374 – − 1.46772 − 2.42717 – − 1.12654 Trihelix transcription factor PTL − 3.93523 − 11.22499 Trihelix transcription factor PTL − 2.06298 Trihelix transcription factor GT-3a 3.24021 – − 2.10905 Trihelix transcription factor PTL − 4.16242 − 11.06445 WUSCHEL-related homeobox 8 3.21184 WUSCHEL-related homeobox 5 − 9.79807 Probable WRKY transcription factor 11 WRKY transcription factor 44 − 2.06571 Probable WRKY transcription factor 65 2.67975 WRKY transcription factor 6 1.48518 Probable WRKY transcription factor 70 4.0632 WRKY transcription factor 6 1.93507 Probable WRKY transcription factor 51 6.68728 10.00602 Probable WRKY transcription factor 70 2.12805 Probable WRKY transcription factor 48 2.43677 Probable WRKY transcription factor 46|Probable WRKY transcrip5.41009 tion factor 30 Probable WRKY transcription factor 28 − 2.06225 Probable WRKY transcription factor 33 3.91559 Probable WRKY transcription factor 51 4.02861

and probable WRKY transcription factor 30, 33, 41, 46, 48, 51, 53, 65 and 70. The WRKY gene family has been suggested to participate in cold stress responses in Arabidopsis (Lee et al. 2005), rice (Ramamoorthy et al. 2008), Glycine max (Zhou et al. 2008), Vitis vinifera (Wang et al. 2014) and cassava (Wei et al. 2016). In particular, cucumber WRKY46 has conferred cold tolerance to transgenic plants and positively regulated the cold signaling pathway in an ABA-dependent manner (Zhang et al. 2016). Cassava WRKY65 showed up-regulation at all time points of cold treatment (Wei et al. 2016).

13

The expression of 12 genes was analyzed by Real time quantitative PCR (qRT-PCR) during cold stress to confirm results obtained from the RNA sequencing. These included MYC4, bHLH14, ZAT12, ZAT10, ERF061, ERF109, GATA5, GRAS1, MYB44, NAC100, NAC89, WRKY30 and DREB1D and their qRT-PCR expression patterns were well in agreement with the RNA sequencing results (Fig. 4). In addition, many TFs such as ARR18 (ARR-B family), Zinc finger protein CONSTANS-LIKE 2 (CO-like family), DOF1.2, DOF5.7, DOF3.7, Cyclic dof factor 2 and Dof zinc finger protein DOF2.1 (Dof family), Putative Myb family

Genes & Genomics Table 3  (continued) Tracking ID

AGI locus

Family

Description

log2 (fold change) LT2 vs. CT LT24 vs. CT

XLOC_028564 XLOC_028870 XLOC_029051 XLOC_033150 XLOC_033162

AT1G30650 AT2G38470 AT3G56400 AT4G01250 AT4G11070

WRKY WRKY WRKY WRKY WRKY

XLOC_033412 XLOC_035667 XLOC_043365 XLOC_044857 XLOC_048780 XLOC_004948 XLOC_026622 XLOC_011510 XLOC_033876 XLOC_039254

AT1G30650 AT2G38470 AT5G13080 AT1G80840 AT1G80840 AT1G08465 AT1G08465 AT5G15210 AT1G69600 AT4G24660

WRKY WRKY WRKY WRKY WRKY YABBY YABBY ZF-HD ZF-HD ZF-HD

Probable WRKY transcription factor 14 Probable WRKY transcription factor 33 Probable WRKY transcription factor 70 WRKY transcription factor 22 Probable WRKY transcription factor 53|Probable WRKY transcription factor 41 Probable WRKY transcription factor 14 Probable WRKY transcription factor 33 Probable WRKY transcription factor 75 – – Putative axial regulator YABBY 2 Putative axial regulator YABBY 2 Zinc-finger homeodomain protein 8 Zinc-finger homeodomain protein 11 –

− 2.46396

− 1.43419 − 9.73092

− 2.2637 − 1.85843

− 1.88767

− 5.86493 4.19123 4.04332 1.98289 3.22583

4.00902 − 9.86842 5.4644 4.79753 − 4.71246 − 1.88649 − 1.85038

This table corresponds to the Additional file 7 described with the reference genome sequence information (Tang et al. 2016). Tracking ID is gene or transcript ID generated and assigned by Cuffmerge (Trapnell et al. 2012a, b) in this transcriptome analysis. The blank means non-differential expression (|log2 (fold change)| < 1)

transcription factor At1g14600 (G2-like family), GATA 5 and 8 (GATA family), Chitin-inducible gibberellin-responsive protein 1, Scarecrow-like protein 9, 13 and 33 (GRAS family), Growth-regulating factor 3 (GRF family), Homeobox-leucine zipper protein HAT14 (HD-ZIP family), LOB domain-containing protein 27 (LBD family), MADS-box transcription factor 13 (MIKC-MADS family), Agamouslike MADS-box protein AGL62 (M-type MADS family), AP2/ERF and B3 domain-containing transcription factor At1g51120 (RAV family), Homeobox protein SBH1 and BEL1-like homeodomain protein 2 (TALE family), Transcription factor TCP2 (TCP family), Trihelix transcription factor GT-3a (Trihelix family), WUSCHEL-related homeobox 8 (WOX family), and NAC transcription factor 25, NAC transcription factor NAM-B2, Protein FEZ, NAC domain-containing protein 8, 78, 89 and 100 (NAC family)

were demonstrated to be up-regulated in present study. However, to our knowledge, there is no report of these transcription factors involved in cold stress response, suggesting that they might be specific to rubber tree.

Conclusions This study describes the genome-wide transcriptional response of rubber tree (H. brasiliensis) to the low-temperature treatments. A large number of DEGs were characterized including 239 transcription factors. The results will help to develop a better understanding the molecular basis of cold stress response in rubber tree, thus creating a start point for improving the cold stress tolerance of rubber tree with gene manipulation.

13



Genes & Genomics

Fig. 4  Expression patterns of 12 DEGs confirmed by qRT-PCR

Acknowledgements  This work was funded by the National Natural Science Foundation of China (No. 31560573) and the Hainan Province Major Science and Technology Project (ZDZX2013023). Author contributions  XXG, BYY and JH conducted bioinformatics studies and carried out qRT-PCR validation. CPY and YJL participated in collection of plant materials. JPL supervised the experiments and wrote the manuscript. WBL participated in the bioinformatics studies. All the authors read and approved the manuscript for publication.

Compliance with ethical standards  Conflict of interest  Xiao-Xiao Gong declares that she does not have conflict of interest. Bing-Yu Yan declares that she does not have conflict of interest. Jin Hu declares that he does not have conflict of interest. Cui-Ping Yang declares that she does not have conflict of interest. Yi-Jian Li declares that he does not have conflict of interest. Jin-Ping Liu declares that he does not have conflict of interest. Wen-Bin Liao declares that he does not have conflict of interest.

13

Ethical approval  This article does not contain any studies with human subjects or animals performed by any of the authors.

References Ahmed NU, Park JI, Jung HJ, Hur Y, Nou IS (2015) Anthocyanin biosynthesis for cold and freezing stress tolerance and desirable color in Brassica rapa. Funct Integr Genom 15:383–394 An D, Yang J, Zhang P (2012) Transcriptome profiling of low temperature-treated Cassava apical shoots showed dynamic responses of tropical plant to cold stress. BMC Genom 13:64 Azuma A, Yakushiji H, Koshita Y, Kobayashi S (2012) Flavonoid biosynthesis-related genes in grape skin are differentially regulated by temperature and light conditions. Planta 236:1067–1080 Bahieldin A, Atef A, Edris S, Gadalla NO, Ali HM, Hassan SM, AlKordy MA, Ramadan AM, Makki RM, Al-Hajar AS et al (2016) Ethylene responsive transcription factor ERF109 retards PCD and improves salt tolerance in plant. BMC Plant Biol 16:216

Genes & Genomics Brown RS (2005) Zinc finger proteins: getting a grip on RNA. Curr Opin Struct Biol 15:94–98 Calzadilla PI, Maiale SJ, Ruiz OA, Escaray FJ (2016) Transcriptome response mediated by cold stress in Lotus japonicus. Front Plant Sci 7:374 Chalker-Scott L (1999) Environmental significance of anthocyanins in plant stress responses. Photochem Photobiol 70:1–9 Chen B-Q, Cao J-H, Wang J-K, Wu Z-X, Tao Z-L, Chen J-M, Yang C-A, Xie G-S (2012) Estimation of rubber stand age in typhoon and chilling injury afflicted area with Landsat TM data: a case study in Hainan Island, China. For Ecol Manag 274:222–230 Chinnusamy V, Zhu J, Zhu JK (2007) Cold stress regulation of gene expression in plants. Trends Plant Sci 12:444–451 Chinnusamy V, Zhu JK, Sunkar R (2010) Gene regulation during cold stress acclimation in plants. Methods Mol Biol 639:39–55 Christie PJ, Alfenito MR, Walbot V (1994) Impact of low-temperature stress on general phenylpropanoid and anthocyanin pathways: enhancement of transcript abundance and anthocyanin pigmentation in maize seedlings. Planta 194:541–549 Croft D, O’Kelly G, Wu G, Haw R, Gillespie M, Matthews L, Caudy M, Garapati P, Gopinath G, Jassal B et al (2011) Reactome: a database of reactions, pathways and biological processes. Nucleic Acids Res 39:D691-D697 Ding J, Wang L, Ruan C (2017) Comparative transcriptome analysis of lipid biosynthesis in seeds and non-seed tissues of sea buckthorn. Genes Genom 39:1021–1033 Du Z, Zhou X, Ling Y, Zhang Z, Su Z (2010) agriGO: a GO analysis toolkit for the agricultural community. Nucleic Acids Res 38:W64–W70 Eremina M, Rozhon W, Poppenberger B (2016) Hormonal control of cold stress responses in plants. Cell Mol Life Sci 73:797–810 Fernández-Calvo P, Chini A, Fernández-Barbero G, Chico JM, Gimenez-Ibanez S, Geerinck J, Eeckhout D, Schweizer F, Godoy M, Franco-Zorrilla JM et al (2011) The Arabidopsis bHLH transcription factors MYC3 and MYC4 are targets of JAZ repressors and act additively with MYC2 in the activation of jasmonate responses. Plant Cell 23:701–715 Figueiredo DD, Barros PM, Cordeiro AM, Serra TS, Lourenço T, Chander S, Oliveira MM, Saibo NJM (2012) Seven zinc-finger transcription factors are novel regulators of the stress responsive gene OsDREB1B. J Exp Bot 63:3643–3656 Fowler S, Thomashow MF (2002) Arabidopsis transcriptome profiling indicates that multiple regulatory pathways are activated during cold acclimation in addition to the CBF cold response pathway. Plant Cell 14:1675–1690 Gilmour SJ, Sebolt AM, Salazar MP, Everard JD, Thomashow MF (2000) Overexpression of the Arabidopsis CBF3 transcriptional activator mimics multiple biochemical changes associated with cold acclimation. Plant Physiol 24:1854–1865 Gilmour SJ, Fowler SG, Thomashow MF (2004) Arabidopsis transcriptional activators CBF1, CBF2, and CBF3 have matching functional activities. Plant Mol Biol 54:767–781 Gronover CS, Wahler D, Prüfer D (2011) Natural rubber biosynthesis and physic-chemical studies on plant derived latex. In: Elnashar M (ed) Biotechnology of biopolymers. InTech, Croatia, pp 75–88 Hannah MA, Heyer AG, Hincha DK (2005) A global survey of gene regulation during cold acclimation in Arabidopsis thaliana. PLoS Genet 1:e26 He J, Yang H-J, Lin M-X (1982) Studies on cellular injuries in rubber trees under chilling stress. Acta Sci Nat Univ Amoi 21:84–91 Hu Y, Jiang Y, Han X, Wang H, Pan J, Yu D (2017) Jasmonate regulates leaf senescence and tolerance to cold stress: crosstalk with other phytohormones. J Exp Bot 68:1361–1369 Hua W, Kong W, Cao X-Y, Chen C, Liu Q, Li X, Wang Z (2017) Transcriptome analysis of Dioscorea zingiberensis identifies genes involved in diosgenin biosynthesis. Genes Genom 39:509–520

Huang ZD, Pan YQ (1992) Rubber cultivation under climatic stresses in China. In: Sethuraj MR, Mathew NM (eds) Natural rubber: biology, cultivation and technology. Elsevier, Dordrecht, pp 220–238 Jaradat MR, Feurtado JA, Huang D, Lu Y, Cutler AJ (2013) Multiple roles of the transcription factor AtMYBR1/AtMYB44 in ABA signaling, stress responses, and leaf senescence. BMC Plant Biol 13:192 Jin J, Zhang H, Zhang J, Liu P, Chen X, Li Z, Xu Y, Lu P, Cao P (2017) Integrated transcriptomics and metabolomics analysis to characterize cold stress responses in Nicotiana tabacum. BMC Genom 18:496 Kan L-Y, Xie G-S, Tao Z-L, Yang L-F, Cui Z-F (2009) Analysis on rubber tree cold injury in 2007/2008 winter in Hainan. Chin Agric Sci Bull 25:251–257 Kanehisa M, Goto S (2000) KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res 28:27–30 Kaplan F, Kopka J, Sung DY, Zhao W, Popp M, Porat R, Guy CL (2007) Transcript and metabolite profiling during cold acclimation of Arabidopsis reveals an intricate relationship of cold-regulated gene expression with modifications in metabolite content. Plant J 50:967–981 Karp PD, Ouzounis CA, Moore-Kochlacs C, Goldovsky L, Kaipa P, Ahrén D, Tsoka S, Darzentas N, Kunin V, López-Bigas N (2005) Expansion of the BioCyc collection of pathway/genome databases to 160 genomes. Nucleic Acids Res 33:6083–6089 Khaledian Y, Maali-Amiri R, Talei A (2015) Phenylpropanoid and antioxidant changes in chickpea plants during cold stress. Russ J Plant Physiol 62:772–778 Knight MR, Knight H (2012) Low-temperature perception leading to gene expression and cold tolerance in higher plants. New Phytol 195:737–751 Lee BH, Henderson DA, Zhu JK (2005) The Arabidopsis coldresponsive transcriptome and its regulation by ICE1. Plant Cell 17:3155–3175 Li J, Yang S, Gai J (2017) Transcriptome comparative analysis between the cytoplasmic male sterile line and fertile line in soybean (Glycine max (L.) Merr.). Genes Genom 39:1117–1127 Lin M-X, Yang H-J (1994) Physiological responses of Hevea brasiliensis during the chilling injury. Chin J Trop Crops 15:7–11 Lindemose S, O’Shea C, Jensen MK, Skriver K (2013) Structure, function and networks of transcription factors involved in abiotic stress responses. Int J Mol Sci 14:5842–5878 Liu Q, Kasuga M, Sakuma Y, Abe H, Miura S, Yamaguchi-Shinozaki K, Shinozaki K (1998) Two transcription factors, DREB1 and DREB2, with an EREBP/AP2 DNA binding domain, separate two cellular signal transduction pathways in drought- and low temperature responsive gene expression, respectively, in Arabidopsis. Plant Cell 10:1391–1406 Liu J-P, Xia Z-Q, Tian X-Y, Li Y-J (2015) Transcriptome sequencing and analysis of rubber tree (Hevea brasiliensis Muell.) to discover putative genes associated with tapping panel dryness (TPD). BMC Genom 16:398 Liu J-P, Zhuang Y-F, Guo X-L, Li Y-J (2016) Molecular mechanism of ethylene stimulation of latex yield in rubber tree (Hevea brasiliensis) revealed by de novo sequencing and transcriptome analysis. BMC Genom 17:257 Luo P, He J-J, Yao Y-L, Dai X-H, Cheng R-X, Wang L-F (2012) Differential responses of two rubber tree clones to chilling stress. Afr J Biotech 11:13466–13471 Mai J, Herbette S, Vandame M, Kositsup B, Kasemsap P, Cavaloc E, Julien J, Ameglio T, Roeckel-Drevet P (2009) Effect of chilling on photosynthesis and antioxidant enzymes in Hevea brasiliensis Muell. Arg Trees 23:863–874 Matthews L, Gopinath G, Gillespie M, Caudy M, Croft D, de Bono B, Garapati P, Hemish J, Hermjakob H, Jassal B et al (2009)

13

Reactome knowledgebase of human biological pathways and processes. Nucleic Acids Res 37:D619–D622 Mizoi J, Shinozaki K, Yamaguchi-Shinozaki K (2012) AP2/ERF family transcription factors in plant abiotic stress responses. Biochim Biophys Acta 1819:86–96 Niu Y, Figueroa P, Browse J (2011) Characterization of JAZ-interacting bHLH transcription factors that regulate jasmonate responses in Arabidopsis. J Exp Bot 62:2143–2154 Park S, Lee CM, Doherty CJ, Gilmour SJ, Kim Y, Thomashow MF (2015) Regulation of the Arabidopsis CBF regulon by a complex low-temperature regulatory network. Plant J 82:193–207 Persak H, Pitzschke A (2014) Dominant repression by Arabidopsis transcription factor MYB44 causes oxidative damage and hypersensitivity to abiotic stress. Int J Mol Sci 15:2517–2537 Petrussa E, Braidot E, Zancani M, Peresson C, Bertolini A, Patui S, Vianello A (2013) Plant flavonoids—biosynthesis, transport and involvement in stress responses. Int J Mol Sci 14:14950–14973 Priyadarshan PM (2011) Biology of hevea rubber. CAB International, Wallingford, pp 1–163 Provart NJ, Gil P, Chen W, Han B, Chang HS, Wang X, Zhu T (2003) Gene expression phenotypes of Arabidopsis associated with sensitivity to low temperatures. Plant Physiol 132:893–906 Pushparajah E (1983) Problems and potentials for establishing Hevea under difficult environmental conditions. Planter 59:242–251 Ramamoorthy R, Jiang SY, Kumar N, Venkatesh PN, Ramachandran S (2008) A comprehensive transcriptional profiling of the WRKY gene family in rice under various abiotic and phytohormone treatments. Plant Cell Physiol 49:865–879 Roberts A, Pimentel H, Trapnell C, Pachter L (2011) Identification of novel transcripts in annotated genomes using RNA-SEq. Bioinformatics 27:2325–2329 Roy S (2016) Function of MYB domain transcription factors in abiotic stress and epigenetic control of stress response in plant genome. Plant Signal Behav 11:e1117723 Schaefer CF, Anthony K, Krupa S, Buchoff J, Day M, Hannay T, Buetow KH (2009) PID: the pathway interaction database. Nucleic Acids Res 37:D674–D679 Schulz E, Tohge T, Zuther E, Fernie AR, Hincha DK (2016) Flavonoids are determinants of freezing tolerance and cold acclimation in Arabidopsis thaliana. Sci Rep 6:34027 Seo J, Sohn H, Noh K, Jung C, An J, Donovan C, Somers D, Kim D, Jeong S-C, Kim C-G et al (2012) Expression of the Arabidopsis AtMYB44 gene confers drought/salt-stress tolerance in transgenic soybean. Mol Breed 29:601–608 Sharma M, Laxmi A (2016) Jasmonates: emerging players in controlling temperature stress tolerance. Front Plant Sci 6:1129 Solecka D (1997) Role of phenylpropanoid compounds in plant responses to different stress factors. Acta Physiol Plant 19:257–268 Song S, Qi T, Fan M, Zhang X, Gao H, Huang H, Wu D, Guo H, Xie D (2013) The bHLH subgroup IIId factors negatively regulate jasmonate-mediated plant defense and development. PLoS Genet 9:e1003653 Stockinger EJ, Gilmour SJ, Thomashow MF (1997) Arabidopsis thaliana CBF1 encodes an AP2 domain-containing transcription activator that binds to the C repeat/DRE, a cis-acting DNA regulatory element that stimulates transcription in response to low temperature and water deficit. Proc Natl Acad Sci USA 94:1035–1040 Strickler SR, Bombarely A, Mueller LA (2012) Designing a transcriptome next generation sequencing project for a nonmodel plant species. Am J Bot 99:257–266 Tang C, Yang M, Fang Y, Luo Y, Gao S, Xiao X, An Z, Zhou B, Zhang B, Tan X et al (2016) The rubber tree genome reveals new insights into rubber production and species adaptation. Nat Plants 2:16073

13

Genes & Genomics Theocharis A, Clément C, Barka EA (2012) Physiological and molecular changes in plants grown at low temperatures. Planta 235:1091–1105 Thomas PD, Campbell MJ, Kejariwal A, Mi H, Karlak B, Daverman R, Diemer K, Muruganujan A, Narechania A (2003) PANTHER: a library of protein families and subfamilies indexed by function. Genome Res 13:2129–2141 Thomashow MF (1999) Plant cold acclimation freezing tolerance genes and regulatory mechanisms. Annu Rev Plant Physiol 50:571–599 Tian Y-H, Yuan H-F, Xie J, Deng J-W, Dao X-S, Zheng Y-L (2016) Effect of diurnal irradiance on night-chilling tolerance of six rubber cultivars. Photosynthetica 54:374–380 Trapnell C, Roberts A, Goff L, Pertea G, Kim D, Kelley DR, Pimentel H, Salzberg SL, Rinn JL, Pachter L (2012a) Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks. Nat Protoc 7:562–578 Trapnell C, Hendrickson DG, Sauvageau M, Goff L, Rinn JL, Pachter L (2012b) Differential analysis of gene regulation at transcript resolution with RNA-sEq. Nat Biotechnol 31:46–53 Van Buskirk HA, Thomashow MF (2006) Arabidopsis transcription factors regulating cold acclimation. Physiol Plant 126:72–80 Vannini C, Locatelli F, Bracale M, Magnani E, Marsoni M, Osnato M, Mattana M, Baldoni E, Coraggio I (2004) Overexpression of the rice Osmyb4 gene increases chilling and freezing tolerance of Arabidopsis thaliana plants. Plant J 37:115–127 Varshney RK, Nayak SN, May GD, Jackson SA (2009) Next-generation sequencing technologies and their implications for crop genetics and breeding. Trends Biotechnol 27:522–530 Vogel JT, Zarka DG, Van Buskirk HA, Fowler SG, Thomashow MF (2005) Roles of the CBF2 and ZAT12 transcription factors in configuring the low temperature transcriptome of Arabidopsis. Plant J 41:195–211 Wang Y-R, Liu H-X, Guo Z-Y (1978) Effects of chilling temperature of plant metabolism of Hevea brasiliensis. Acta Bot Sin 20:44–53 Wang X-J, Li W-G, Gao X-S, Wu C-T, Zhang X-F (2012) Physiological characteristics of Hevea brasiliensis in response to low temperature stress and its regulation mechanisms. Plant Physiol J 48:318–324 Wang L, Zhu W, Fang L, Sun X, Su L, Liang Z, Wang N, Londo JP, Li S, Xin H (2014) Genome-wide identification of WRKY family genes and their response to cold stress in Vitis vinifera. BMC Plant Biol 14:103 Wang H, Wang H, Shao H, Tang X (2016) Recent advances in utilizing transcription factors to improve plant abiotic stress tolerance by transgenic technology. Front Plant Sci 7:67 Wang P, Ma L, Li Y, Wang S, Li L, Yang R (2017a) Transcriptome analysis reveals sunflower cytochrome P450 CYP93A1 responses to high salinity treatment at the seedling stage. Genes Genom 39:581–591 Wang X, Zhao M, Wu W, Korir NK, Qian Y, Wang Z (2017b) Comparative transcriptome analysis of berry-sizing effects of gibberellin (GA3) on seedless Vitis vinifera L. Genes Genom 39:493–507 Wei Y, Shi H, Xia Z, Tie W, Ding Z, Yan Y, Wang W, Hu W, Li K (2016) Genome-wide identification and expression analysis of the WRKY gene family in cassava. Front Plant Sci 7:25 Wu ZG, Jiang W, Chen SL, Mantri N, Tao ZM, Jiang CX (2016) Insights from the cold transcriptome and metabolome of Dendrobium officinale: global reprogramming of metabolic and gene regulation networks during cold acclimation. Front Plant Sci 7:1653 Xie C, Mao X, Huang J, Ding Y, Wu J, Dong S, Kong L, Gao G, Li CY, Wei L (2011) KOBAS 2.0: a web server for annotation and identification of enriched pathways and diseases. Nucleic Acids Res 39:W316-322 Yang QS, Gao J, He WD, Dou TX, Ding LJ, Wu JH, Li CY, Peng XX, Zhang S, Yi GJ (2015) Comparative transcriptomics analysis

Genes & Genomics reveals difference of key gene expression between banana and plantain in response to cold stress. BMC Genom 16:446 Zeng C, Chen Z, Xia J, Zhang K, Chen X, Zhou Y, Bo W, Song S, Deng D, Guo X et al (2014) Chilling acclimation provides immunity to stress by altering regulatory networks and inducing genes with protective functions in cassava. BMC Plant Biol 14:207 Zhan X, Yang L, Wang D, Zhu JK, Lang Z (2016) De novo assembly and analysis of the transcriptome of Ocimum americanum var. pilosum under cold stress. BMC Genom 17:209 Zhang Y, Yu H, Yang X, Li Q, Ling J, Wang H, Gu X, Huang S, Jiang W (2016) CsWRKY46, a WRKY transcription factor from

cucumber, confers cold resistance in transgenic-plant by regulating a set of cold-stress responsive genes in an ABA-dependent manner. Plant Physiol Biochem 108:478–487 Zhao C, Lang Z, Zhu JK (2015) Cold responsive gene transcription becomes more complex. Trends Plant Sci 20:466–468 Zhou Q-Y, Tian A-G, Zou H-F, Xie ZM, Lei G, Huang J, Wang CM, Wang HW, Zhang JS, Chen SY (2008) Soybean WRKY type transcription factor genes, GmWRKY13, GmWRKY21, and GmWRKY54, confer differential tolerance to abiotic stresses in transgenic Arabidopsis plants. Plant Biotech J 6:486–503

13