Transcriptome Profiles of Populus euphratica upon Heat Shock stress.

Heat stress, which strongly affects plant performance and often results in reduced vegetative growth and yields depression, has become an increasingly serious global problem. Populus euphratica Oliv. which has been considered as a tree model for the study of plant response to abiotic stresses, could be resistant to an extremely wide environmental temperature range (-40 °C to 45 °C). Previous study is mainly focused on its gene regulation upon drought and salt stress. However, little is known about gene regulation at the global transcriptome level upon heat stress. To understand the gene network controlling heat stress in P. euphratica, a transcriptome sequencing using Illumina Hiseq 2000 was performed to generate a 10 gigabases depth for each sample in the tissue of leaf. 119,573 unigeneswere generated with an average length of 474 bp. Approximately 49,605 (41.49%) unigenes exhibited significantly different expressions between two libraries. Among these unigenes, 11,165 (9.34%) were upregulated and 38,440 (32.15%) were down regulated. Heat shock proteins classified as molecular chaperones showed a significant percentage (1.13%) in the up regulated group. Heat responsive genes, such as polyubiquitins, were over expressed in heat treated sample. GO enrichment analysis revealed that the Go terms for differentially expressed unigenes were significantly enriched in hormone-mediated signal, biological process regulation and metabolic process regulation. Our data revealed a global transcriptome picture of P. euphratica in response to heat shock. The identified potential heat stress-related transcripts can be used to infer the gene regulation networks underlying the molecular mechanisms of heat response in P. euphratica.

INTRODUCTION

Previous reports have predicted that greenhouse gases, such as CO2, methane, chlorofluorocarbons, and nitrous oxides, gradually increase the global ambient temperature [1]. Heat stress is a serious global threat to plant growth and production because of high ambient temperatures [2]. When exposed to heat stress, severe cellular injury and even cell death may take place in a very short time; this adverse condition can be attributed to the disruption of cellular organization [3]. The changes at molecular level take place immediately after stress happened and result in altering gene expression and transcript accumulation. Thus, stress-related proteins are synthesized to induce stress tolerance [4]. Previous studies have established an important adaptive strategy, in which almost all of the organisms respond to heat shock by synthesizing heat shock proteins (HSPs) [1]. HSPs induce tolerance and improve physiological processes, including photosynthesis, assimilate partitioning, and membrane stability [1, 5].

Although HSP stress response strategy is conserved in prokaryotes and eukaryotes [6], only a few plant species or genotypes share abilities to cope with heat stress. High variations in heat stress-coping strategies among plants provide opportunities to study the mechanism of heat stress tolerancewithin and between species. In recent studies, deep sequencing technologies have been widely used [7, 8]. Using this method, scholars can capture expressed sequence tags (ESTs) and identify novel transcripts in a specific tissue at a particular point at a whole genome level without bias [9, 10].

Populus euphratica Oliv. is naturally distributed in semiarid areas and play important roles in maintaining local arid ecosystems [11-14]. As a tree species model widely used to elucidate abiotic resistance mechanisms [15-17], P. euphratica can tolerate temperatures as high as 45 °C. To obtain heat tolerance genes and investigate the mechanisms involved in heat stress response of this species, we presented a de novo assembly transcriptome of P. euphratica exposed to stress at 45 °C by using Solexa data. The acquired information may contribute to the development of strategies and facilitate improvement of heat tolerance of trees.

MATERIALS AND METHODS

Sample Preparation

Plant materials were collected from two-year-old P. euphratica seedlings obtained from the Xinjiang Autonomous Region of China and planted in conditions as previously reported [13]. These potted plants were well irrigated at an interval of 3 d according to the evaporation rate two months before treatment. Heat treatment was designed as described in (Fig. ). Plants were acclimated by heating the chamber at a constant, non-lethal temperature (37 °C, 2 h); the temperature was then increased linearly from 37 °C to 45 °C for 2 h and then stay at 45 °C for 3 h. At the end of the treatment, the leaves from four different plants were harvested, frozen immediately in liquid nitrogen, and stored at –80 °C for RNA extraction [13].

H2O2 Content Assay

The level of H2O2 was measured by monitoring the A415 of the titanium–peroxide complex according to the method reported by [7]. Absorbance values were calibrated to a standard curve generated with known concentrations of H2O2.

RNA Isolation and Illumina cDNA Library Construction

Total RNA was extracted using the CTAB [18]. After treated with RNAase free DNAase, the A260/A280 ratios were examined by NanoDrop 2000 and ratios ranging from 1.9 to 2.1 were selected [13]. The integrity of the RNA samples was examined using an Agilent 2100 Bioanalyzer and their RNA integrity number (RIN) values were >8.0 [13].

To construct the Illumina Hiseq 2000 libraries, 25 μg RNA sample with a concentration of ≥750 ng/μl was used for each cDNA library construction. Poly (A) mRNA was initially enriched with oligo(dT) and subsequently fragmented into small pieces of 200 bp to 700 bp by using divalent cations at an elevated temperature. Based on these cleaved RNA fragments, random hexamer-primer and reverse transcriptase (Invitrogen) were used to synthesize cDNA [13, 19]. Two paired-end cDNA libraries with an insert size of 200 bp were constructed and subsequently sequenced using Illumina Hiseq™ 2000 [13, 19].

De Novo Assembly and Assessment

Raw data generated from Solexa sequencing were preprocessed to remove non-sense sequences, including adapters, sequences with numerous unknown bases (>5%), and sequences with low-quality bases (>50% of the bases with a quality score ≤5), by using an in-house Perl script [19]. The preprocessed sequences were then assembled by using Trinity (Version: r2011-08-20) program [20]. Reads were first combined with certain overlap lengths to form fragments known as contigs,and then these contigs were further realigned to construct unigenes by Trinity. To fill the intra-scaffold gaps, we used the paired-end information to retrieve read pairs that contained one well-aligned read to the contigs and another read located in the gap region [21]; the collected reads were locally assembled. After the gap was closed, we constructed a non-redundant unigene set from the two assembled datasets by TGICL program [22]. To decide the sequential orientation of each unigene, we performed a set of sequential BLASTx alignment (E < 1e-5) against the non-redundant (NR) database of GenBank, Swiss-Prot protein database (http://www.expasy.ch/sprot), Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database [23], and Clusters of Orthologous Groups of proteins (COG) database (http://www.ncbi.nlm.nih.gov/COG/). For unigenes that cannot be aligned to any of these databases, ESTScan was used to determine the sequence orientation [24].

Unigene Annotation and Function Classification

All assembled unique sequences were searched against the databases of NR, NT, Swiss-Prot, KEGG, and poplar transcripts version3 by BLASTn (E < 1e-5) to gain the most descriptive annotation. The protein with the highest sequence similarity was retrieved. Pathway mapping analysis was performed based on KEGG database alignment. According to Nr annotation, GO information was obtained using the Blast2GO server [25]. GO functional classification was performed using the WEGO software with plants categories defined by molecular function, cellular component, and biological process ontologies [26].

Gene Expression Analysis

The expression level of P. euphratica unigene was normalized according to the number of RPKM [27, 28]. After normalization was conducted, we determined the number of reads in each coding region in the control and treated libraries. The ratio of the reads in the two libraries was then calculated. The statistical significance of the differential expression value of each unigene was determined according to a previously described method [29]; the results of the statistical tests were corrected for multiple testing by using the Benjamini and Hochberg false FDR correction [30]. To eliminate the effect of the highly sensitive Solexa/Illumina sequencing method, we defined DE unigenes as those with an absolute value of log2Ratio ≥ 2, FDR < 0.001 as well as those expressed at ≥ 3 RPKM in ≥1 sample. The unigenes that failed to contain a minimum of three RPKM in at least one sample were removed [31].

To obtain the significantly enriched GO terms of DE unigenes on the genome background, we annotated all of the DE unigenes to the GO database by using a hypergeometric test [32]. The cutoff p value after correction was 0.05 on rigorous Bonferroni correction method. GO terms meeting these standards were defined as significantly enriched GO terms. The KEGG pathway enrichment analysis of DE unigenes was performed to determine the main biochemical pathways and signal transduction pathways involving DE unigenes. Corrected p ≤ 0.05 was used as a threshold to identify the overrepresented pathways.

Quantitative PCR Analysis

Quantitative PCR (qPCR) was conducted using a power SYBR Green PCR kit (ABI) in a MicroAmp™96-well plate with a StepOnePlus™ Real-time PCR System (ABI). The relative quantification value was calculated according to the 2-∆∆CT method [33]. Before qPCRanalyses, we illustrate the suitability of several reference gene including RPL17, TUB and 60S and finally chose PeActin (GenBank Accession Number: EF148840) as an internal control [34]. Reactions were prepared using a total volume of 20 µl containing 10 µl of 2× SYBR premix, 2 µl of cDNA template, and 0.5 µl of each specific primer to obtain a final concentration of 200 nM. The reactions were performed under the following conditions: 94 °C for 2 min; 45 cycles of 94 °C for 20 s; 60 °C for 35 s; and 68 °C for 1 min. The gene-specific primers used in the qPCR analysis are listed in (Table ). For PCR analyses, three independent samples were used which were different from those used for RNA-seq. Three technical replicates were performedfor each sample as shown in previous report [33].

RESULTS AND DISCUSSION

Solexa Sequencing Results

The total RNA obtained from approximately 30 leaves of four treated or non-treated plants was used to create biologically independent pools to capture heat-responsive genes. After the adapters, low-quality sequences, and ambiguous reads were removed [12], a total of 64.6 million and 68.9 million clean reads with a mean length of 90 bp were generated in the control group (CK) and heat shock-treated sample (HS4), respectively (Table ). The raw data has been submitted to the NCBI Sequence Read Archive database with the accession number SRP029139. All of the trimmed reads were de novo assembled into contigs according to the Trinity method [20]. The average contig size was 135 and 134 bp in CK and HS4 libraries, respectively. Using paired-end information, we joined the contigs to obtain a total of 119,573 unigenes with an average length of 474 bp and a N50 of 548 bp (Table ). Length distribution analysis results showed that 9.74% (11,643) of these unigenes were longer than 1,000 bp (Fig. ).

Functional Annotation of the All Unigenes

Using BLASTx or BLASTn with an e-value of <1.0E–5, we annotated all of the unigenes according to the public databases, including NR database, nucleotide sequence database, the poplar genome (v3), Swiss Prot database, COG and KEGG databases. Among the119,573 high quality unique sequences, 93,214 (78%) showed at least one significant match to an existing gene model in the BLAST search. The COG classification and functional classifications GO terms of all unigenes were shown in (Fig. ). In the biological process category of GO terms, the five largest groups included cellular process, metabolic process, stimulus response process, pigmentation process, and biological regulation process. In the molecular function category, the unigenes with binding and catalytic activities formed the largest groups (Fig. ).

Differential Expression Unigenes Identification

According to our applied criteria (FDR < 0.001, absolute value of log2Ratio ≥ 1, and RPKM ≥ 3 in ≥ 1 sample), approximately 41.49% (49,605 unigenes) of the total unigenes were identified as differentially expressed (DE) unigenes, in which 11,165 were up regulated and 38,440 were down regulated (Fig. ). Among these DE unigenes, 43,827 had a gene model in poplar database (Table ).

To identify DE unigenes accurately, we selected the top 100 most up regulated transcripts and the top 100 most down regulated transcripts with good annotation (FDR < 0.001, Table ). Among the top 100 most upregulated unigenes, 98 transcripts could match a transcript in the genome database of P. trichocarpa. Those encoding HSPs (e.g., Unigene119328, Unigene105373, and Unigene119280), heat shock factors (HSFs, e.g., Unigene119554, Unigene103826, and Unigene100969), dehydration responsive element binding protein (DREB) transcription factors (e.g. Unigene 109904, Unigene118327, Unigene117802, and Unigene 118202), and ubiquitin/polyubiquitin (e.g., Unigene110702, Unigene118782, Unigene111816, Unigene117269, and Unigene117112) occupied 62% of the most up regulated interpretable transcripts in the HS4 sample (Table ). This result indicates a concentrated function annotation. Furthermore, 49 chaperone proteins such as HSP and DnaJ were included in these 100 unigenes. The high percentage of HSPs and HSFs in the top 100 most up regulated gene list suggested that HSPs were the primary factors affecting the heat tolerance of P. euphratica, and also indicated the efficiency of our treatments.

Distinct from those of up regulated unigenes, the unigenes with diverse functions were identified when the top 100 most downregulated analysis was performed. Some zinc finger proteins, cytochrome P450 genes, AP2/ERF-domain containing transcription factors, and protein kinases were listed in the top 100 most down regulated unigenes. The remaining clusters of unigenes primarily encoded abiotic stress-related transcription factors, general metabolism involving enzymes, and auxin-induced proteins (Table ).

We further performed qPCR analysis to validate the reliability of our sequencing results. Ten members were randomly selected from the top 100 most up regulated and the top 100 most down regulated DE unigene list. The qPCR results indicated that all of the DE unigenes exhibited similar expression kinetics to the deep sequencing results (Table ).

KEGG Mapping and Gene Ontology Enrichment Analyses for DE Unigenes

Based on a comparison against the KEGG database, twenty-five pathways changed significantly (p ≤ 0.05) under heat shock stress, including the pathways involved in carbohydrate, energy, lipid, vitamin, hormone, and pyruvate metabolisms (Table ). Plant pathogen interaction, N-glycan biosynthesis, and zeatin biosynthesis pathways were three of the top five most differentially expressed pathways that seem to play important roles in heat shock response.

GO enrichment analysis showed biological processes related to biological regulation and hormone-mediated signaling were significantly enriched (p ≤ 0.05, Table ). This result suggested that these processes are important in heat stress response. Unigenes related to endogenous stimulus and chemical stimulus responses were also enriched for DE unigenes, indicating a comprehensive change in gene expression in heat-treated P. euphratica. We also identified the enriched GO terms for upregulated DE unigenes based on genomic background. GO terms including stress response process, temperature stimulus response process, oxidative stress response process, and reactive oxygen species (ROS) response were enriched (Fig. ). The increased transcripts in these systems may indicate their important roles in protection of P. euphratica under conditions of heat stress by removal of or response to ROS.

Heat Shock Proteins and Heat Shock Factors in Response to Heat Shock Stress

HSPs that function as molecular chaperones are usually induced by heat shock at any developmental stage of higher plants [35]. HSPs are considered as candidate genes providing protective traits in plants against abiotic stresses because HSPs function in stabilizing proteins, repairing protein structures, and assisting in protein refolding [36]. HSPs are classified into five families according to molecular size: Hsp100, Hsp90, Hsp70, Hsp60, and small HSPs (sHSP) [37, 38]. According to our dataset, numerous HSPs were significantly induced after heat shock stress treatment. Besides 49 HSPs presented in the top 100 most up regulated unigene list (Table ), 93 members designated as heat shock-related genes were up regulated while 28 members were down regulated from the whole database (Table ). Among these unigenes, nine sequences encoding HSFs, which are known ascentral regulators of heat stress responsive genes [39], were up regulated. Eight sHSPs were up regulated, and none was down regulated, indicating a vital role played by sHSP in heat shock stress response. A highly conserved HSP gene named HSA32 has been reported in different plant species. However, only slightly increase of this gene has been found according our dataset. Close inspection of this gene revealed a higher abundance after 37 °C treatment, indicating its early response to heat in P. euphratica (data not shown). Heat shock cognate proteins (HSCs), DnaJ (HSP40) homologs and other HSPs were often classified as chaperone proteins. HSC70 is required to activate HSF1, which appears to control the early response of target genes in heat stress. DnaJ exhibits the same capability as HSP70 in reactivating proteins in bacteria and may respond to various environmental stresses independently as a molecular chaperone [40]. In our study, three sequences encoding HSC proteins were found and two of them were up regulated. In contrast to HSCs, DnaJ-like proteins comprise more down regulated members than up regulated ones, indicating that P. euphratica may initially utilize HSP70 and then DnaJs to resist heat stress. All these results indicated chaperone proteins function as master players for heat stress tolerance.

DREB Transcription Factors in Response to Heat Shock Stress

Some drought-responsive transcription factors, such as DREB, were overrepresented (FDR < 0.001) in the HS4 sample. Seven sequences with DREB annotation were identified in our dataset. Among these sequences, six were included in the top 100 most up regulated list. Based on Pfam database searches of predicted open reading frames and NCBI blast program, two genes, i.e., PeDREB2 (GenBank Accession No. EF137176.1; encoded by unigene 110992, unigene 109904 and unigene 117802) and PeDREB2L (GenBank Accession No. EF567422.1; encoded by unigene 118203, unigene 118327 and unigene 118202) were identified. DREB2-type transcription factors have been considered as important regulon not only in drought and salt stress signaling but also in heat stress response signaling in Arabidopsis [41, 42]. DREB2A can bind to HsfA3 promoter directly and function as a regulator of HsfA3 [43]. Here, we showed that PeDREB2 and PeDREB2L can be significantly induced by heat shock (Tables ). Furthermore, two DREB2A-interacting protein genes that encode DRIP2 (unigene 52587 and unigene 38695), which were supposed to function negatively in response to abiotic stress, were dramatically down regulated in the HS4 sample (Table ). Considering the salt response function of PeDREB2 and drought response function of PeDREB2L we have reported [34, 41], we concluded that a complex cross-talk between different abiotic stress responses is involved in P. euphratica and that the DREB2-cellular signal transduction is crucial to the heat adaptation of this tree species. However, further studies should be conducted to determine whether these two factors or only one of them can participate in heat shock response.

Heat Stress Increases Intracellular Accumulation of H2O2

Many physiological damages occur when plants exposing to heat stress, e.g., causing ROS accumulation. DREB2A and HSFA2, which were highly induced after heat shock stress in our transcript profile, were also known to be involved in hydrogen peroxide (H2O2)-response signaling or strongly induced by H2O2 treatment [44]. HSPs, which encode cytosolic proteins for the detoxification of H2O2 and function as targets of HSFA2, were also overrepresented after heat shock stress. These results indicated that many H2O2-response genes may be implicated in heat shock stress. We further examined the H2O2 level before and after heat shock treatment. As shown in (Fig. ), an increase intracellular accumulation of H2O2 has been investigated for heat stress treated plants, suggesting heat stress have induced oxidative production. Heat shock-induced H2O2 may facilitate transcripts accumulation of heat shock responsive genes in P. euphratica.

We also investigated oxidative stress-responsive Ascorbate Peroxidases (APXs) that may be involved in heat stress response. Previously, a barley APX1 has been found to play roles in heat tolerance. A functional heat shock elements have been found in the 5’-promoter region of some Arabidopsis APX genes, indicating their potential role in heat stress response [37]. However, no APX gene has yet been investigated to be unregulated in our study. We supposed that some heat-responsive gene may be induced in different stages or in different tissues.

Calcium/calmodulin-mediated Signal Network Involved in Heat Stress

As a second universal messenger involved in abiotic stress responses in plants, calcium/calmodulin-mediated signal network functions as an important messenger in thermomemory associated signal transduction [45]. Calmodulin has been reported to be involved in heat shock signal transduction in wheat [46]. Calcium signals are perceived by calcium sensors such as calmodulins (CaMs) and calcium-dependent protein kinases (CDPK) that relay the signals into downstream targets in different stimuli [47].

As a regulatory protein involved in cellular calcium-dependent signaling pathways, CDPKs are monomeric proteins containing a CaM-like domain with four EF-hand motifs [48]. The CDPK family is a large family with various functions participating in different stress signal pathways. In the present study, four unigenes (Unigene24235, Unigene24346, Unigene20655, and Unigene113411) encoding different CDPK proteins were significantly up regulated after heat shock treatment, suggesting that CDPK might be an important component in heat shock response of P. euphratica. CaMs are small proteins that transmit the Ca2+ signal by interacting with target proteins and regulating their activity. Several CaMs, CMLs and dozens of CaM-binding protein genes were showed upregulated versus the control group. Particularly, some CaM-binding HSPs were found to be induced over 10-fold by heat shock treatment (Table ). Thus, we propose that CaM might play important roles during the expression of HSPs through binding directly to cytoplasmic HSPs. To define the functional significance of these Ca2+ sensors, loss of-function and gain-of-function analyses for each of these genes under heat shock stress remains to be studies in the future.

Polyubiquitin in P. euphratica Heat Shock Response

Ubiquitin, a highly conserved protein of 76 amino acids present in eukaryotic cells, is involved in protein degradation via a multicatalytic proteasome complex [49, 50]. The ubiquitin system has been implicated in immune response, development, and programmed cell death [49]. As one of the two ubiquitin classes, polyubiquitins may have important functions in heat shock response in P. euphratica because of their overexpression (FDR < 0.001) in the top 100 up regulated gene list (Table ). A yeast polybiquitin gene ubi4 is induced in response to environmental stresses, such as heat shock, starvation, and oxidative stress [50]. Ubi1 is another polybiquitin gene from an entophytic fungus that shows increased transcription at high temperatures; this result indicates a stress reaction of hyphal cells after heat shock treatment [49]. Although we cannot explain the mechanism of the overexpression of unigenes encoding polyubiquitin proteins in the heat shock-treated library precisely, the most probable explanation is that heat shock stress treatment induced damage to some proteins, which are targeted for proteolysis. In higher plants, the best-characterized outcome of ubiquitination is the process mediating target protein degradation via 26S proteasome, in which polyubiquitin is essential for the 26S proteasome recognition [51]. Evidence suggests that ubiquitination may have a critical function in regulating plant responses to abiotic stresses and has prompted researchers to conduct further studies to identify ubiquitin genes that mediate plant tolerance of abiotic stress. However, the mechanism by which polyubiquitin functions in response to adverse environmental conditions remains unknown. A recent study on the genome-wide transcriptome analyses of Populus simonii also established the crucial role that DREB played in heat response. However, few up-regulated polyubiquitin genes have been found in their study [52]. Basing on our sequencing results, we identified three polyubiquitin genes with full coding sequences against poplar genome database. Considering that all of these genes contain a common heat shock element sequence (CNNGAANNTTCNNG) in their upstream promoter region [53], we concluded that the up regulated polyubiquitin genes in P. euphratica play important roles in heat shock response.

CONCLUSION

In this work, we have analyzed global transcriptome profiling of P. euphratica upon heat shock stress.

Heat responsive genes, such as heat shock proteins, DREB transcription factors and polyubiquitins, were significantly up regulated in heat treated sample. Go terms for differentially expressed unigenes were significantly enriched in hormone-mediated signal, biological process regulation and metabolic process regulation. Our results provide a global picture of P. euphratica in response to heat shock stress at transcriptomic level by using Illumina Hiseq 2000. Detailed characterization of these heat shock responsive genes will be studied in future.

Table 1.

Overview of the sequencing and assembly.

	HS4	CK	All
Total Clean Reads	64,611,114	68,888,892
Total Clean Nucleotides	5,815,000,260	6,200,000,280
Q20 percentage	96.30%	96.20%
N percentage	0.01%	0.01%
GC percentage	44.78%	44.59%
Contig
Total Number	738,481	620,759
Total Length(nt)	99,461,549	82,873,556
Mean Length(nt)	135	134
N50	116	116
Unigene
Total Number	130,421	109,624	119,573
Total Length(nt)	47,236,810	39,455,066	56,628,880
Mean Length(nt)	362	360	474
N50	410	406	548

Table 2.

Top 100 most upregulated unigenes of P. euphratica by heat shock stress treatment.

Gene ID	Poplar v3 Model	Log2 (Fold Change)	Annotation
Unigene113351	Potri.001G042700.1	12.52	heat shock protein 70 cognate
Unigene119328	Potri.008G062300.1	11.44	18.2 kDa class I heat shock protein
Unigene118411	Potri.001G042700.1	11.41	heat shock protein 70 cognate
Unigene100860	Potri.004G187400.1	11.32	low molecular weight heat shock protein
Unigene117613	Potri.001G042700.1	11.28	heat shock protein 70 cognate
Unigene119280	Potri.010G195700.1	11.13	18.2 kDa class I heat shock protein
Unigene100276	Potri.004G187400.1	11.12	heat-shock protein
Unigene26446	Potri.009G049800.1	11.12	17.5 kd heat shock protein GmHSP17.6L
Unigene117236	Potri.009G049900.1	11.05	Hsp20.1 protein
Unigene118998	Potri.004G073600.1	11.04	heat shock protein
Unigene110823	Potri.010G195700.1	11.00	18.2 kDa class I heat shock protein
Unigene109379	Potri.009G049900.1	10.92	heat shock protein
Unigene100392	Potri.009G147900.1	10.84	17.5 kDa class I heat shock protein
Unigene36975	Potri.010G053400.1	10.78	chloroplast small heat shock protein
Unigene119447	Potri.006G223900.1	10.69	heat shock protein 17.7 - garden pea
Unigene119010	Potri.013G089200.1	10.61	22.0 kDa class IV heat shock protein
Unigene106116	Potri.019G081200.1	10.59	17.5 kd heat shock protein GmHSP17.6L
Unigene106960	Potri.003G167500.1	10.45	Bcl-2-associated athanogene-like protein
Unigene119446	Potri.003G109200.2	10.33	LMW heat shock protein
Unigene115745	Potri.012G022400.1	10.30	heat shock protein
Unigene36846	Potri.006G093500.1	10.15	18.2 kDa class I heat shock protein
Unigene116444	Potri.003G071100.1	9.94	17.4 kDa class III heat shock protein
Unigene116523	Potri.017G084000.1	9.93	FtsH protease
Unigene115042	Potri.004G187200.1	9.75	cytosolic class I small heat shock protein type 2
Unigene110702	Potri.004G211600.1	9.68	ubiquitin-like protein 5
Unigene106090	Potri.005G259900.1	9.61	calmodulin-like protein
Unigene106270	Potri.017G084000.1	9.32	FtsH protease
Unigene117683	Potri.017G084000.1	9.20	FtsH protease; putative
Unigene115386	Potri.003G071100.1	9.17	17.4 kDa class III heat shock protein
Unigene26606	Potri.019G081200.1	9.11	17.5 kd heat shock protein Gmhsp17.6L
Unigene105373	Potri.008G054000.1	8.98	heat shock protein 70
Unigene101586	Potri.010G206600.1	8.89	heat shock protein 70
Unigene113828	Potri.010G206600.1	8.84	heat shock protein 70
Unigene118352	Potri.015G056900.1	8.79	endopeptidaseClp (EC 3.4.21.-) ATP-binding chain SB100
Unigene105915	Potri.015G057000.1	8.59	heat shock protein
Unigene106538	Potri.010G206600.1	8.57	heat shock protein 70
Unigene118982	Potri.002G166300.1	8.43	BAG6; calmodulin binding / protein binding
Unigene107019	Potri.015G004800.1	8.37	Avr9/Cf-9 rapidly elicited protein 65
Unigene119573	Potri.002G166300.1	8.25	BAG6; calmodulin binding / protein binding
Unigene111455	Potri.010G088600.1	7.92	heat shock protein 70 family protein
Unigene119411	Potri.017G146600.1	7.80	heat shock protein
Unigene117980	Potri.009G039200.1	7.77	17.6 kDa class I heat shock protein;Hsp20.0
Unigene34332	Potri.012G021600.1	7.69	Avr9/Cf-9 rapidly elicited protein 65
Unigene117779	Potri.011G057600.1	7.67	dnaJ subfamily B member 5
Unigene119554	Potri.006G226800.3	7.56	AtHSFA2
Unigene118622	Potri.013G018000.1	7.51	heat shock protein 70 cognate
Unigene101449	Potri.004G034400.2	7.50	F-box family protein; late embryogenesis abundant protein
Unigene49867	Potri.005G214800.2	7.31	heat stress transcription factor A-6b
Unigene103826	Potri.005G214800.2	7.18	heat stress transcription factor A-6b
Unigene100924		7.10	NADH dehydrogenase subunit 6
Unigene100969	Potri.005G214800.2	7.01	heat stress transcription factor A-6b
Unigene31571	Potri.017G146600.1	7.01	Heat shock protein 83
Unigene118961	Potri.004G213400.1	6.97	Rubisco subunit binding-protein alpha subunit
Unigene118390	Potri.005G183300.1	6.96	calcium-binding protein
Unigene115732	Potri.002G191600.1	6.89	putative galactinol synthase
Unigene118782	Potri.017G135600.3	6.88	polyubiquitin
Unigene117269	Potri.017G135600.1	6.88	polyubiquitin
Unigene119548	Potri.001G182100.1	6.80	Probable pyridoxin biosynthesis PDX1-like protein 2
Unigene26267	Potri.013G018000.1	6.71	heat shock protein 70 cognate
Unigene111816	Potri.017G135600.1	6.52	polyubiquitin
Unigene100739	Potri.001G289700.1	6.52	EGY3, ethylene-dependent gravitropism-deficient and yellow-green-like 3
Unigene117112	Potri.017G135600.1	6.45	polyubiquitin
Unigene39330		6.43	nad6;ubiquinone oxidoreductase chain 6;NADH dehydrogenase subunit 6
Unigene102654	Potri.002G048200.1	6.42	AtHSFA7A
Unigene102474	Potri.004G034400.1	6.36	F-box family protein; late embryogenesis abundant protein
Unigene119268	Potri.001G289700.1	6.09	EGY3, ethylene-dependent gravitropism-deficient and yellow-green-like 3
Unigene103979	Potri.001G303600.1	6.06	phenylpropanoid:glucosyltransferase 1
Unigene118877	Potri.010G168300.1	6.04	zinc finger (C3HC4-type RING finger) family protein
Unigene26854	Potri.003G216100.1	6.04	oxidoreductase/ transition metal ion binding protein
Unigene116435	Potri.011G085100.2	6.03	SGS domain-containing protein; calcyclin binding protein
Unigene110992	Potri.008G073600.2	6.03	Putative dehydration responsive element binding protein 2H
Unigene25167	Potri.011G122500.1	6.01	WD-40 repeat family protein
Unigene111646	Potri.010G168300.1	5.98	zinc finger (C3HC4-type RING finger) family protein
Unigene113318	Potri.014G044300.2	5.87	chaperonin 10
Unigene110294	Potri.016G003400.3	5.83	HSP80
Unigene4724	Potri.011G085100.1	5.80	SGS domain-containing protein; calcyclin binding protein
Unigene109904	Potri.008G073600.2	5.75	dehydration responsive element binding protein 2H
Unigene118305	Potri.014G116800.3	5.73	putative galactinol synthase
Unigene110176	Potri.014G116800.3	5.72	putative galactinol synthase
Unigene26341	Potri.017G130700.2	5.70	peroxisomal small heat shock protein
Unigene112432	Potri.014G056600.1	5.69	four F5 protein-related / 4F5 protein-related
Unigene118203	Potri.010G183700.1	5.67	dehydration responsive element binding protein
Unigene111511	Potri.001G286700.1	5.66	HSP80
Unigene35342	Potri.004G213400.1	5.63	Rubisco subunit binding-protein alpha suunit
Unigene116440	Potri.007G136700.1	5.56	DNAJ chaperone C-terminal domain-containing protein
Unigene114084	Potri.002G191600.1	5.54	putative galactinol synthase
Unigene119368	Potri.018G075200.1	5.52	ethylene-responsive transcriptional coactivator
Unigene118327	Potri.010G183700.1	5.47	Dehydration responsive element binding protein
Unigene114471	Potri.016G003400.2	5.46	heat shock protein 90-2
Unigene117802	Potri.008G073600.2	5.41	Dehydration responsive element binding protein
Unigene116777	Potri.008G101600.1	5.41	gibberellin 2-oxidase
Unigene25884	Potri.001G381000.1	5.38	5-azacytidine resistance protein -related
Unigene107042	Potri.014G184200.1	5.38	Retrotransposon gag protein
Unigene34276	Potri.016G120000.1	5.37	DNAJ heat shock N-terminal domain-containing protein
Unigene113416	Potri.006G073500.1	5.32	low molecular weight heat-shock protein
Unigene109968	Potri.014G116800.3	5.24	putative galactinol synthase
Unigene44055	Potri.009G099000.1	5.21	glucosyltransferase
Unigene43602	Potri.001G154200.1	5.18	AP2/ERF domain-containing transcription factor
Unigene33231	Potri.019G083800.1	5.17	ubiquitin conjugating-like enzyme
Unigene118202	Potri.010G183700.1	5.12	dehydration responsive element binding protein 2H

Table 3.

Top 100 most downregulated unigenes of P. euphratica by heat shock stress treatment.

Gene ID	Poplar v3 Model	log2 (Fold Change)	Annotation
Unigene47078	Potri.011G150100.1	-7.21	oxidoreductase, 2OG-Fe(II) oxygenase family protein
Unigene16607	Potri.019G131600.2	-6.10	RPT2
Unigene46509	Potri.012G090000.1	-5.89	cytochrome P450
Unigene35116	Potri.016G071900.1	-5.84	aspartyl protease family protein
Unigene43437	Potri.004G148400.1	-5.67	transducin family protein / WD-40 repeat family protein
Unigene118552	Potri.007G121100.1	-5.64	DNA binding protein
Unigene27890	Potri.008G061800.1	-5.59	protease inhibitor/seed storage/LTP family protein
Unigene51359	Potri.004G162600.1	-5.53	salt tolerance-like protein
Unigene48152	Potri.015G027700.1	-5.52	UDP-glucose:isoflavone 7-O-glucosyltransferase
Unigene15296	Potri.001G118400.1	-5.50	cytochrome P450
Unigene3678	Potri.002G039100.1	-5.42	AP2/ERF domain-containing transcription factor
Unigene24729	Potri.012G091500.2	-5.41	similar to Probable pectatelyase 22 precursor
Unigene26237	Potri.012G083300.1	-5.41	alliinase family protein
Unigene41914	Potri.002G144200.2	-5.38	GRAS family transcription factor
Unigene1323	Potri.001G071000.1	-5.34	xyloglucanendotransglucosidase
Unigene42138	Potri.011G147900.1	-5.31	zinc finger (C3HC4-type RING finger) family protein
Unigene25314	Potri.006G232400.2	-5.27	41 kD chloroplast nucleoid DNA binding protein (CND41)
Unigene47358	Potri.014G020500.1	-5.20	cytochrome P450
Unigene118351	Potri.010G120100.1	-5.17	leucine-rich repeat receptor-like protein kinase
Unigene46286	Potri.008G065000.1	-5.15	transferase family protein
Unigene15583	Potri.004G152300.1	-5.15	polyol transporter
Unigene49328	Potri.007G121100.1	-5.14	zinc finger (B-box type) family protein;
Unigene4718	Potri.005G195600.1	-5.13	peroxidase 4
Unigene23436	Potri.016G071900.1	-5.12	aspartyl protease family protein
Unigene47056	Potri.001G052300.1	-5.09	Similar to pectatelyase.
Unigene18168	Potri.019G119300.2	-5.00	leucine-rich repeat family protein
Unigene22853	Potri.015G068200.1	-5.00	Germin-like protein subfamily T member 2
Unigene4523	Potri.001G083900.1	-4.99	cytochrome P450
Unigene19863	Potri.001G079000.1	-4.98	Ferric reductase-like transmembrane component
Unigene110599	Potri.004G049100.1	-4.98	HSL1; ATP binding protein serine/threonine kinase
Unigene29860	Potri.001G124200.1	-4.96	glycerol 3-phosphate permease
Unigene41440	Potri.010G042100.3	-4.91	polygalacturonase
Unigene17878	Potri.002G085100.3	-4.91	nodulin family protein
Unigene33652	Potri.016G038000.1	-4.88	type-a response regulator
Unigene104488	Potri.016G083900.1	-4.87	MYB transcription factor
Unigene19269	Potri.013G103800.1	-4.86	phenylcoumaranbenzylic ether reductase –like protein
Unigene47697	Potri.016G083900.1	-4.86	MYB transcription factor
Unigene25745	Potri.001G254100.2	-4.85	gibberellin induced protein
Unigene108845	Potri.001G399600.8	-4.82	Protein SPA1-RELATED 3
Unigene12317	Potri.003G150700.1	-4.79	AP2/ERF domain-containing transcription factor
Unigene24847	Potri.017G076800.1	-4.78	nitrate transporter NRT1-5
Unigene114796	Potri.T126700.1	-4.72	wall-associated receptor kinase-like 10-like, partial
Unigene35916	Potri.001G348100.1	-4.71	glycosyltransferase UGT90A7
Unigene109572	Potri.017G134900.1	-4.71	similar to Putative receptor protein kinase TMK1 precursor
Unigene20711	Potri.009G053800.1	-4.69	DiT1 (dicarboxylate transporter 1); oxoglutarate:malateantiporter
Unigene50973	Potri.007G097500.2	-4.68	RING1; protein binding / ubiquitin-protein ligase/ zinc ion binding
Unigene41578	Potri.005G117100.2	-4.68	zinc finger (B-box type) family protein
Unigene47294	Potri.004G027300.1	-4.66	Very similar to receptor-like protein kinase
Unigene31087	Potri.011G151200.2	-4.66	subtilase family protein
Unigene26479	Potri.003G151000.1	-4.65	AP2/ERF domain-containing transcription factor
Unigene115429	Potri.002G065300.1	-4.63	peroxidase
Unigene110741	Potri.016G016800.1	-4.62	Anthocyanidin 3-O-glucosyltransferase 6
Unigene42213	Potri.009G169900.2	-4.59	3-hydroxy-3-methylglutaryl coenzyme A reductase
Unigene115367	Potri.004G057700.2	-4.59	cysteine protease
Unigene25400	Potri.001G399600.8	-4.59	SPA4 (SPA1-RELATED 4); protein binding / signal transducer
Unigene30850	Potri.002G248500.1	-4.58	basic helix-loop-helix (bHLH) family protein;
Unigene44683	Potri.003G207200.2	-4.57	basic helix-loop-helix (bHLH) family protein;
Unigene24456	Potri.015G068200.1	-4.57	germin-like protein
Unigene44617	Potri.003G101800.1	-4.56	lipase
Unigene34224	Potri.017G097800.1	-4.56	FERONIA receptor-like kinase
Unigene48675	Potri.001G211000.7	-4.56	F-box family protein;
Unigene12711	Potri.009G097300.2	-4.54	cyclase-like protein
Unigene25153	Potri.011G110200.1	-4.54	SIEP1L protein
Unigene48872	Potri.014G149700.1	-4.53	pectinesterase family protein
Unigene11709	Potri.019G119300.2	-4.53	leucine-rich repeat family protein
Unigene115331	Potri.014G094800.1	-4.53	similar to arabinogalactan-protein (AGP20);
Unigene103692	Potri.012G031100.4	-4.52	SIGE (SIGMA FACTOR E); DNA binding / DNA-directed RNA polymerase
Unigene118975	Potri.011G076400.1	-4.51	GDSL-motif lipase/hydrolase-like protein
Unigene34284	Potri.008G166200.1	-4.51	AP2/ERF domain-containing transcription factor
Unigene26224	Potri.001G161400.6	-4.51	plasma membrane proton ATPase
Unigene32382	Potri.016G047500.1	-4.50	Hypothetical protein
Unigene51122	Potri.001G211000.7	-4.49	F-box family protein;
Unigene14480	Potri.009G113600.1	-4.48	polyol transporter
Unigene36234	Potri.006G086100.1	-4.47	expansin
Unigene20944	Potri.013G115800.1	-4.47	lectin protein kinase family protein
Unigene22167	Potri.002G039100.1	-4.47	AP2/ERF domain-containing transcription factor
Unigene2790	Potri.013G115800.1	-4.47	lectin protein kinase family protein
Unigene46536	Potri.019G054800.3	-4.47	Similar to thioredoxin f precursor - garden pea.
Unigene102005	Potri.010G236800.1	-4.47	GDSL-motif lipase/hydrolase family protein
Unigene111472	Potri.004G084000.1	-4.47	leucine-rich repeat receptor-like protein kinase
Unigene20679	Potri.014G037900.1	-4.46	cytochrome P450
Unigene19398	Potri.011G150400.1	-4.45	oxidoreductase, 2OG-Fe oxygenase family protein, expressed
Unigene105153	Potri.006G047500.1	-4.45	bifunctionalmonodehydroascorbate reductase and carbonic anhydrase nectarin-3
Unigene38315	Potri.004G097000.1	-4.44	serine/threonine protein kinase family protein
Unigene27298	Potri.008G152300.2	-4.43	with no lysine kinase
Unigene20610	Potri.018G009200.1	-4.43	ATCNGC17
Unigene47384	Potri.014G060500.1	-4.42	GRAS family transcription factor
Unigene25829	Potri.019G057500.4	-4.41	similar to Alpha-expansin 8 precursor;AtEXPA8
Unigene4810	Potri.010G046300.1	-4.41	low affinity inorganic phosphate transporter
Unigene3581	Potri.008G152300.2	-4.40	with no lysine kinase
Unigene21605	Potri.T079300.1	-4.40	cytochrome P450
Unigene2029	Potri.011G164200.1	-4.40	oxidoreductase, 2OG-Fe(II) oxygenase family protein
Unigene102830	Potri.010G112200.6	-4.40	ATMPK15
Unigene36053	Potri.001G079600.1	-4.39	AP2/ERF domain-containing transcription factor
Unigene27223	Potri.014G026700.1	-4.38	subtilisin-like protease preproenzyme
Unigene45648	Potri.011G066900.4	-4.37	cysteine protease
Unigene15515	Potri.001G210100.2	-4.37	protease inhibitor/seed storage/lipid transfer protein family protein
Unigene23053	Potri.019G085400.1	-4.37	C2-H2 zinc finger protein
Unigene14649	Potri.008G088300.2	-4.36	alpha-expansin 3
Unigene110047	Potri.012G087800.3	-4.36	ribitol dehydrogenase-like/short-chain dehydrogenase/reductase protein

Table 4.

Quantitative RT-PCR for 10 selective DEGs.

Gene ID	Poplar v3 Model	Annotation	45°C Abundance
Unigene117613	Potri.001G042700.1	Heat shock protein 70	8.97±0.29
Unigene117112	Potri.017G135600.1	Polyubiquitin UBQ14	8.64±0.68
Unigene118203	Potri.010G183700.1	DRE-binding protein	6.87±0.53
Unigene117779	Potri.011G057600.1	DnaJ subfamily B member 5	7.69±0.56
Unigene119554	Potri.006G226800.3	Heat shock transcription factor A2	7.32±1.39
Unigene 100969	Potri.005G214800.2	Heat shock transcription factor A6B	6.82±0.21
Unigene 115732	Potri.002G191600.1	Galactinol synthase 4	6.21±0.17
Unigene114471	Potri.016G003400.2	Heat shock protein 90-2	6.53±0.35
Unigene 111646	Potri.010G168300.1	Zinc finger family protein	5.59±0.89
Unigene113416	Potri.006G073500.1	Low molecular weight heat shock protein	5.23±0.10

All P. trichocarpa V3 gene models were obtain by blasting with a threthhold of 1e-5.qPCR was performed on 10 members randomly selected from top 100 upregulated unigenes of the HS4 DEG list.

Table 5.

Enrichment analysis for KEGG metabolic pathway of differentially expressed unigenes.

Pathway ID	Pathway	DE Unigenes (13611)		Corrected p value
		Up	Down
ko04626	Plant-pathogen interaction	199	1065	1.48E-10
ko00510	N-Glycan biosynthesis	1	129	3.45E-07
ko04141	Protein processing in endoplasmic reticulum	105	280	6.58E-07
ko00945	Stilbenoid, diarylheptanoid and gingerol biosynthesis	29	155	4.64E-05
ko00908	Zeatin biosynthesis	29	98	1.21E-04
ko03022	Basal transcription factors	12	115	1.21E-04
ko00100	Steroid biosynthesis	9	57	8.00E-04
ko00905	Brassinosteroid biosynthesis	2	31	8.78E-04
ko00402	Benzoxazinoid biosynthesis	13	46	1.25E-03
ko03040	Spliceosome	75	637	2.27E-03
ko00944	Flavone and flavonol biosynthesis	4	63	2.43E-03
ko00910	Nitrogen metabolism	11	103	2.60E-03
ko00941	Flavonoid biosynthesis	48	142	5.21E-03
ko00903	Limonene and pinene degradation	21	161	6.98E-03
ko00511	Other glycan degradation	5	64	1.03E-02
ko00966	Glucosinolate biosynthesis	8	51	1.56E-02
ko02010	ABC transporters	11	142	2.19E-02
ko00603	Glycosphingolipid biosynthesis - globo series	0	16	2.29E-02
ko00040	Pentose and glucuronate interconversions	17	99	2.38E-02
ko00531	Glycosaminoglycan degradation	3	39	2.68E-02
ko00604	Glycosphingolipid biosynthesis - ganglio series	1	23	3.62E-02
ko00600	Sphingolipid metabolism	4	69	3.75E-02
ko00940	Phenylpropanoid biosynthesis	44	260	3.92E-02
ko00450	Selenocompound metabolism	11	66	4.07E-02
ko04712	Circadian rhythm - plant	13	157	4.33E-02

The cutoff p value after correction was 0.05 on Q value correction method.

Table 6.

Enrichment analysis for GO terms of differentially expressed unigenes

Gene Ontology Term	Cluster Frequency	Genome Frequency of Use	Corrected p-value
Biological Process
hormone-mediated signaling	304 out of 10455 genes, 2.9%	458 out of 19477 genes, 2.4%	2.10E-05
regulation of biological process	1849 out of 10455 genes, 17.7%	3221 out of 19477 genes, 16.5%	0.00238
biological regulation	2185 out of 10455 genes, 20.9%	3842 out of 19477 genes, 19.7%	0.00652
response to endogenous stimulus	610 out of 10455 genes, 5.8%	1014 out of 19477 genes, 5.2%	0.01496
regulation of cellular process	753 out of 10455 genes, 7.2%	1270 out of 19477 genes, 6.5%	0.02336
regulation of metabolic process	1282 out of 10455 genes, 12.3%	2220 out of 19477 genes, 11.4%	0.03066
response to chemical stimulus	1133 out of 10455 genes, 10.8%	1954 out of 19477 genes, 10.0%	0.03993
Molecular Function
oligosaccharyltransferase activity	30 out of 12519 genes, 0.2%	34 out of 23520 genes, 0.1%	0.0092
monooxygenase activity	152 out of 12519 genes, 1.2%	229 out of 23520 genes, 1.0%	0.02015

The cutoff p value after correction was 0.05 on rigorous Bonferroni correction method.