Temporal Gene Expression Kinetics for Human Keratinocytes Exposed to Hyperthermic Stress

The gene expression kinetics for human cells exposed to hyperthermic stress are not well characterized. In this study, we identified and characterized the genes that are differentially expressed in human epidermal keratinocyte (HEK) cells exposed to hyperthermic stress. In order to obtain temporal gene expression kinetics, we exposed HEK cells to a heat stress protocol (44 °C for 40 min) and used messenger RNA (mRNA) microarrays at 0 h, 4 h and 24 h post-exposure. Bioinformatics software was employed to characterize the chief biological processes and canonical pathways associated with these heat stress genes. The data shows that the genes encoding for heat shock proteins (HSPs) that function to prevent further protein denaturation and aggregation, such as HSP40, HSP70 and HSP105, exhibit maximal expression immediately after exposure to hyperthermic stress. In contrast, the smaller HSPs, such as HSP10 and HSP27, which function in mitochondrial protein biogenesis and cellular adaptation, exhibit maximal expression during the “recovery phase”, roughly 24 h post-exposure. These data suggest that the temporal expression kinetics for each particular HSP appears to correlate with the cellular function that is required at each time point. In summary, these data provide additional insight regarding the expression kinetics of genes that are triggered in HEK cells exposed to hyperthermic stress.

in the human genome. Second, many previous studies made gene expression measurements at a single time point; as a result, temporal gene expression kinetics could not be ascertained. Finally, the majority of past studies examined the response of skin cells to ultraviolet stress; however, few studies have specifically analyzed the transcriptional response of HEK to hyperthermic stress [10].
In this study, we sought to identify and characterize the genes that are differentially expressed in HEK exposed to hyperthermic stress. To obtain temporal gene expression kinetics, we exposed HEK to a heat stress protocol (44 °C for 40 min) and used mRNA microarrays at three time points: 0 h, 4 h and 24 h post-heat treatment, as outlined in Figure 1A. Then, we used bioinformatics software to characterize the chief biological processes and canonical pathways associated with these genes. These data provide valuable new insights that give us a much clearer picture of the genes and intracellular signaling pathways that are triggered in human cells exposed to hyperthermic stress. (B) Cellular MTT viability for keratinocytes exposed to hyperthermic stress. Viability percentage data (relative to untreated samples) was measured for wells with concentration densities ranging from 10,000 to 100,000 cells per well. Viability was measured at 0, 4 and 24 h post-heat treatment. (C) Microarray global gene expression profiles. The data shows that 762, 1,422 and 1,478 total genes were differentially expressed at 0, 4 and 24 h post-exposure, respectively. (D) Gene expression for minimal stress protein, heat shock protein 70 (HSPA6) using qRT-PCR. The mRNA expression fold values were measured for sham and heat treatment groups. For each time point, three samples (S1, S2 and S3) were assessed, each in triplicates. Values were calculated in relation to β-actin and normalized to a separate RNA calibrator.

Keratinocytes Induce an Appreciable Stress Response When Exposed To Hyperthermic Stress
We conducted an initial set of experiments to investigate the impact that our selected hyperthermia stress protocol had on cellular viability. For these studies, HEKs were exposed to hyperthermic stress and cellular metabolic activity was measured at 0, 4 and 24 h post-exposure using MTT assays ( Figure 1B). For all cell densities tested, the data shows that greater than 90% of the cells survive after being exposed to our hyperthermic stress protocol. This data confirms that the HEKs cells were viable post-heat exposure.
After confirming that the HEK cells survived post-exposure, we then sought to test whether our selected hyperthermic stress protocol induced an appreciable transcriptional stress response To investigate the transcriptional stress response, we exposed HEK cells to hyperthermic stress and used microarray gene chips to quantify the expression level of genes at 0, 4 and 24 h post-exposure ( Figure 1A). Since a signature feature of an appreciable stress response is the rapid and marked upregulation of minimal stress genes [12], we conducted microarray and PCR analyses to examine the genes that HEKs express when exposed to hyperthermia ( Figure 1C,D). Using Ingenuity Pathway Analysis (IPA) software, we filtered the microarray data using the following cut-offs: (i) an absolute value of expression magnitude (log 2 ratio relative to control) greater or equal to 1.5; and (ii) a p-value less than or equal to 0.05. Of the 54,675 probe sets tested on each gene chip, we found that 762, 1,422 and 1,478 genes were differentially expressed at the 0, 4 and 24 h time points, respectively ( Figure 1C). The number of upregulated genes increased in a logarithmic fashion from 466 at 0 h, 872 at 4 h and up to 1,127 at 24 h. Interestingly, roughly 60% of the transcripts regulated at 0 and 4 h were upregulated, but by 24 h, ~75% of the differentially expressed genes were upregulated.
Specifically, we found that the heat-treated groups expressed transcripts for many heat shock proteins (Hsps). We also found that the transcripts encoding for heat shock protein 70 (HSPA6) exhibited the greatest increase in expression and the highest level of statistical significance. To validate these microarray results, we then conducted qRT-PCR analyses for HSPA6 (Hsp70) ( Figure 1D). For all of the samples tested, we found that the heat-treated treatment group exhibited statistically significant increases in HSPA6 expression, ranging between 10-and 25-fold compared to the sham group. These results were consistent with the microarray data, which showed that HSPA6 expression increased between 4-and 32-fold.

Temporal Gene Expression Profiles for Keratinocytes Exposed To Hyperthermic Stress
In order to visualize the global distribution of the gene expression profiles, we then created volcano plots for all genes at each time point. In these plots, the level of statistical significance (p-value) for each gene is plotted versus the fold-change in expression level relative to untreated sham ( Figure 2). Genes with the highest level of statistical significance appear on the right side of the plot. In addition, upregulated genes that display large-magnitude fold-changes appear towards the top of the plot (denoted with green triangles), whereas downregulated genes appear towards the bottom (denoted with red triangles). We found that considerably more genes were differentially expressed at 4 h and 24 h than at the 0 h time point. We also found that the transcripts expressed at the later time points also exhibited lower p-values and larger-magnitude changes in expression level than the genes at 0 h. Collectively, our results show that hyperthermic stress triggers a rapid, robust and extensive transcriptional response in HEK cells. The temporal kinetics expression data also show that this stress response is immediately activated in response to stress and is sustained for many hours after exposure.

Figure 2.
Volcano plots of the gene expression profiles at each time point. The magnitude of differential expression (log 2 fold-change) is plotted versus the level of statistical significance (p-value) of all genes in the microarray for 0, 4 and 24 h post-heat exposure. Upregulated and downregulated genes with absolute log 2 ratios 1.5 and p-values ≤0.05 are depicted with green and red triangle symbols, respectively. Insignificant genes are denoted with light gray symbols.

Biological Functions of Heat-Shock-Induced Genes
After identifying the genes expressed in response to hyperthermic stress, we then used IPA software to identify each gene's primary biological functions. Figure 3 contains pie charts of the total number of genes associated with each cellular and molecular function. We found that 23 biological functions were common to all three time points. Interestingly, only nine functions in this group were either unique to one or shared by two time points (Figure 3). In general, the 0 h and 24 h time points share many genes with functions required for free radical scavenging, whereas only the 4 h time point contains genes with functions dedicated primarily to protein regulation (i.e., protein trafficking and protein folding) ( Figure 3B). Next, we identified the biological function groups that had the greatest number of differentially expressed genes. For each time point, we generated color-coded pie charts containing the top five biological functions and the number of genes in each group ( Figure 4). We found that two biological functions (i.e., cellular growth and proliferation and cell death) appeared in the top five at each time point. In addition to these functions, the 0 h and 4 h time points also shared two more top functions, which are associated with cellular development and gene expression. The only major difference between these time points was that the 0 h time point had more genes associated with molecular transport, whereas the 4 h group had more genes involved in cellular movement functions. In contrast, despite exhibiting a similar number of regulated genes, the 4 h and 24 h groups only shared three top functions: cellular movement, cell death and cell growth. The primary difference between these time points is that the 4 h group had more genes dedicated to cellular development and gene expression, whereas the 24 h group exhibited genes with cellular function and cell cycle processes.

Identification of Shared and Unique Genes
We next sought to determine whether any of the differentially expressed genes were observed at all time points. To answer this question, we performed a biomarker filter analysis. This tool allowed us to identify the genes that were unique and shared by each group. The Venn diagram in Figure 5 shows the total number of genes that were shared and unique to each time point. A total of 492, 1,066 and 1,142 genes were unique to the 0 h, 4 h and 24 h evaluation groups, respectively. We found that 121 genes were common to both the 0 h and 4 h, 107 genes to the 0 h and 24 h and 183 genes to the 4 h and 24 h. Notably, only 50 genes were shared among all three recovery time points. A comprehensive list of these genes is provided as supplementary material (Supplementary Tables 1-6). After identifying the 50 genes that were expressed at each time point, we then performed a series of analyses to better understand the function of each gene. For these analyses, we used the following web-based resources: IPA, GeneCards Human Gene Database, HGNC and Gene Ontology (GO). The pie chart in Figure 6 contains the symbols for each of the 50 common genes. This pie chart is also further divided by the family of each gene and the number of genes per family, which is denoted by the number towards the center of each sector of the pie chart. The gene families are color coded and are specified in the legend. The following eight families are provided in the figure: ligand-dependent nuclear receptor, peptidase, phosphatase, transporter, ion channel, kinase, transcription regulator, enzymes and others. The twenty-four genes included in the "other" group consist of genes that were not classified by IPA. Figure 6. Pie chart of 50 shared biomarkers. Each sector of the pie represents the number of genes associated with each function. The symbols and number of genes are also provided for each function. Genes were filtered for an absolute value log 2 ratio 1.5 and a significance value of p ≤ 0.05.
Next, we sought to examine the temporal expression kinetics for the genes with known biological functions. We found that only 44 of the 50 genes identified had well-characterized biological processes and GO annotations. Figure 7 provides the magnitude of the fold-change in expression for each of these genes as a function of time. The asterisks denote genes that have multiple key biological functions. The genes include the chaperones that function in response to stress and protein unfolding. These genes encode for the following heat shock 70 kDa proteins: HSPA6, HSPA1L and HSPA1A/A1B.  Figure 7. The temporal expression kinetics for the genes with known gene ontology (GO). Asterisks indicate the genes with more than one biological process.
Additional genes responsible for the response to stress or heat are CDHR1 (cadherin-related family member 1), TUB (tubby protein homolog), NOX1 (NADPH oxidase 1) and IRF1 (interferon regulatory factor 1). Two genes-PSMB2 (proteasome subunit beta type 2) and UBA6 (ubiquitin-like modifier activating enzyme 6)-play roles in protein ubiquitination. CDHR1, ARHGEF40 (Rho guanine nucleotide exchange factor (GEF) 40), CD9 (CD9 molecule), DST (dystonin), LRRN2 (leucine rich repeat neuronal 2) and PCDH11X/Y (protocadherin 11 Y-linked) support cell adhesion. ALPK1 (alpha-kinase 1), AGK (acylglycerol kinase), CAMK2B (calcium/calmodulin-dependent protein kinase II beta), MYLK (myosin light chain kinase) and TLK2 (tousled-like kinase 2) assist in phosphorylation, whereas, GRID1 (glutamate receptor delta-1 subunit), KCNK10 (potassium channel subfamily K member 10), ARL4C (ADP-ribosylation factor-like 4), HBA1/HBA2 (hemoglobin alpha 1), PITPNC1 (cytoplasmic phosphatidylinositol transfer protein 1), CAMK2B and NOX1 support the transport function. In addition, CADM2 (cell adhesion molecule 2), HIST1H2AG (histone cluster 1, H2ag), NEBL (nebulette) and IGHA1 play roles in junction organization, nucleosome assembly, actin filament length and immune response, respectively. Furthermore, as indicated in Figure 7, some of the above mentioned genes and the remaining ones are involved in processes related to homeostasis response. These include regulation of transcription, apoptosis, metabolism and signal transduction related-genes. The gene families related to the unique genes or the genes shared between each two time points, as well as the number of genes per each family are illustrated in the pie charts in Figure 8. The gene families included those mentioned for the 50 common genes in addition to other gene families, such as G-protein coupled receptor, growth factor, cytokine, translation regulator and transmembrane receptor family. The genes classified by IPA as "other" were omitted in this representation. A list of all these genes and fold-change following heat-shock in the different interval of recovery is added as tables in the supplemental materials. In Figure 9, we highlight a set of genes from the unique genes list that are recognized to play roles in cellular stress response. These include genes involved as heat shock proteins and genes supporting DNA damage and repair, inflammation and apoptosis.  Figure 9. Gene expression profiles for a set of unique genes involved in several well-characterized pathways of the cellular stress response.

Pathway Analysis of Microarray Analyzed Data
We next used IPA to analyze the metabolic and signaling pathways that are affected in response to heat stress following the three post-heat time points. We found that at 0 h, 33 of 247 canonical pathways were expressed with statistical significance, defined as [(−log(p-value) ≥ 1.3 or p-value ≤ 0.05] (    Subsequently, in order to highlight the pathways more relevant to cellular stress response, we further refined our analysis by excluding signaling pathway categories, such as cancer, cardiovascular and nervous system signaling, toxicity list and disease-specific pathways, in addition to non-skin cell-type-specific canonical pathways. Figure 10 shows a refined presentation of the top canonical pathways for each post-heat recovery time point. For the cell cycle, the G1/S checkpoint regulation pathway was affected in all three time points. The ceramide signaling pathway was affected in both 0 h and 4 h. The glucocorticoid receptor signaling pathway and role of CHK proteins in cell cycle checkpoint control were affected in both 0 h and 24 h. All the remaining canonical pathways involved were exclusive to a particular time point. Notably, although most of the canonical pathways in each post-heat time point were unique, they actually serve overlapping top cellular functions. The canonical pathways playing a role in cellular function and maintenance that were affected in 0 h post-heat consisted of MSP-RON signaling, Rac signaling pathways and two endocytosis pathways. Both endocytosis pathways, namely, macropinocytosis and clathrin-mediated endocytosis signaling pathways, play important roles in plasma membrane reorganization, internalization of extracellular molecules and in transduction of signals within the cell and between cells. Deregulation of seven and 13 genes in macropinocytosis signaling and clathrin-mediated endocytosis signaling, respectively, starting 0 h post-heat, implicate an immediate effect of heat shock on the plasma membrane and signal transduction. Macropinocytosis depends on signaling to the actin cytoskeleton, such as the Rac signaling pathway, which was also affected at 0 h post-heat, with eight deregulated genes, including the Rac activator, TIAM1 (T-lymphoma invasion and metastasis 1). Immediately, TIAM1 mRNA showed more than a 4-fold increase. Interestingly, TIAM1 has been reported to prevent cell death in heat-stressed wild-type murine keratinocytes, since the TIAM1 knock-out cells displayed increased apoptosis [20]. Of the canonical pathways affected at 4 h post-heat, tight junction signaling plays a role in cellular function and maintenance, whereas those affected at 24 h post-heat included amino sugar metabolism and cAMP-mediated signaling pathways. The latter is also important in controlling cellular assembly and organization in addition to post-translational modification. Notably, tight junction signaling also

Cell Culture Conditions and Hyperthermia Stress Protocol
HEK were cultured in EpiLife ® medium containing supplement S7, as per the manufacturer's recommendations (Life Technologies). For the heat-shock experiment, cells were put in 1.7 mL tubes (1.0 × 10 6 cells/tube) that were sealed with Para-film ® and placed in a circulating water bath at 44 °C for 40 min. No significant change in pH was observed. Following hyperthermia, cells were either harvested within 5 min or were first allowed to recover in the incubator at 37 °C /95% humidity/5% CO 2 for 4 h or for 24 h. Control cells harvested alongside the heat-treated cells were also put in 1.7 mL tubes (1.0 × 10 6 cells/tube), but were maintained at all time in an incubator at 37 °C /95% humidity/5% CO 2 .

Cellular Viability Assays
Viability was evaluated at different post-exposure time points using MTT (3-(4,5-dimethylthiazol-2-yl)-2,5-diphenyltetrazolium bromide) assays (ATCC, Manassas, VA, USA). In brief, HEKs were exposed to hyperthermic stress, then were incubated for 5 min (0 h post-heat), 4 h (4 h post-heat) or 24 h (24 h post-heat). Subsequently, MTT reagents and detergents were added, as per the manufacturer's instructions. Absorbance readings were measured at 570 nm with a Synergy HT Plate Reader (Biotech, Winooski, VT, USA). Absorbance readings were also measured for a ladder of serial dilutions (10 3 -10 6 cells) in culture medium. Absorbance values were plotted versus the cell number, and these curves were used to determine the number of viable cells in each well.

RNA Isolation and Normalization
RNA was extracted from the heat-treated and control cells using a RNeasy Mini-Kit (Qiagen), according to the manufactures instructions. RNA extraction was performed at three different time intervals as follows: within 5 min following heat-shock (0 h post-heat), after a 4 h recovery period at 37 °C (4 h post-heat) or after a 24 h recovery period at 37 °C (24 h post-heat). At each time point, the RNA concentration was assessed on a NanoDrop Spectrophotometer (NanoDrop Technologies), and the quality was measured on a 2100 Bioanalyzer™ (Agilent Technologies). RNA samples with a RNA Integrity Number greater than 9.5 were subjected to microarray analysis.

Microarray Data Analysis
Expression analysis was performed for each time point in triplicates (three control and three heat-shock treated) using the Affymetrix GeneChip ® Human Genome U133 (HG-U133) plus 2.0 Array that contains 54,675 probe sets. Briefly, two micrograms of RNA were used for preparation of biotin-labeled targets (cRNA) using MessageAmp™-based protocols (Ambion, Inc.). Labeled cRNA was fragmented (0.5 μg/μL per reaction) and was used for array hybridization and washing. The cRNA was mixed with a hybridization cocktail, heated to 99 °C for 5 min and then incubated at 45 °C for 5 min. Hybridization arrays were conducted for 16 h in an Affymetrix Model 640 hybridization oven (45 °C, 60 rpm). Arrays were washed and stained on an FS450 Fluidics station and were scanned on a GeneChip ® Scanner 3000 7G. Image signal data, detection calls and annotations were generated for every gene using the Affymetrix Statistical Algorithm MAS 5.0 (GeneChip ® Operating Software v1.3). A log 2 transformation was conducted and a Student's t-test was performed for comparison of the two groups (control and heat-shocked). We conducted multiple testing correction-Benjamini and Hochberg-to determine the false discovery rate, and statistical significant genes were identified using Bonferroni correction procedures (−log 10 p cutoff > 6.04) [21].
For interpretation of the results, the Ingenuity Pathways Analysis tool (IPA version 8.7, Ingenuity ® Systems Inc., Redwood City, CA, USA; [22]) was used. IPA is a web-based software application, which enables filtering and dataset comparisons, to identify biological mechanisms, pathways and functions most relevant to experimental datasets or differentially expressed genes. The cut-off criteria for our IPA analysis were: an absolute value of log 2 ratio ≥1.5 and a p-value ≤0.05. Other web-based resources, such as the GeneCards Human Gene Database [23], The HUGO Gene Nomenclature Committee (HGNC, [24]) and the Gene Ontology [25] were also used to further supplement the analysis. The targets identified in the microarray study were validated using qRT-PCR. Runs were performed using TaqMan ® RNA-to-CT™ 1-Step Kits and TaqMan ® Assays for: HSPA6 and GAPDH (Applied Biosystems). Calibrator RNA was used as the control (Cell Applications Inc.). All assays were performed by Asuragen Services.

Conclusions
In this study, we examined the temporal gene expression kinetics for HEK cells exposed to hyperthermic stress. Collectively, our results demonstrate that hyperthermic stress triggers a rapid, robust and extensive transcriptional response in HEK cells. The temporal kinetics expression data show that this cellular response is immediately activated in response to stress and is sustained for many hours after exposure. Specifically, we show that the number of upregulated genes increased in logarithmic fashion as a function of time, ranging from 466 at 0 h, 872 at 4 h and up to 1,127 at 24 h. This finding is consistent with previous in vivo studies, which show that HSP70 expression is bi-phasic and peaks at both 12 h and 24 h post-exposure [5,7,8]. However, this data is in contrast to previous in vitro reports, which show that the number of regulated genes is maximal roughly 4 to 12 h post-exposure [3][4][5]. Furthermore, in this study, we show that many of the genes activated at 24 h encode for chemokines that play central roles in inflammation (i.e., IL-8, IL-10, IL-22, IL-15 and JAK family kinases). This suggests that the second expression peak at 24 h is likely due to cellular inflammation.
We identified a large group of diverse heat-inducible genes in this work. The temporal kinetics data shows that each individual gene exhibits distinct expression profiles, which vary dramatically in magnitude, onset time and duration of expression. From the data, we determined that the heat-inducible genes can be grouped according to their temporal expression phases. The first group belongs to an immediate-expression phase consisting of genes that are rapidly expressed shortly (i.e., 5-20 min) after exposure to hyperthermic stress. The second group consists of genes expressed several minutes up to 4 h post exposure, and the third group can be considered an "adaption or recovery phase" that occurs 4 h to 24 h post exposure.
The majority of genes expressed during the immediate expression phase encode for proteins involved in molecular transport, cellular growth and heat shock response. One notable molecular transport gene is HBA1/HBA2, which functions to facilitate the delivery of oxygen. Another interesting finding is that only the well characterized HSPs (i.e., HSP40 and HSP70) are expressed at the early time points, while the smaller HSPs and chaperonins (HSPE1, HSP60; HSPE1, HSP10; HSPB2, HSP27) are expressed at roughly 24 h post-exposure. Furthermore, other HSPs, such as HSPH1 (HSP105), which function primarily to prevent protein aggregation, exhibit maximal expression at the 4 h time-point. This data strongly supports that each particular type of HSP exhibits expression kinetics that appear to be directly correlated with the function that is required at each time point. In summary, these data provide valuable new insights that give us a much clearer picture of the genes and intracellular signaling pathways that are triggered in human epithelial keratinocytes exposed to hyperthermic stress.