Whole-transcriptome gene expression profiling in an epidermolysis bullosa simplex Dowling-Meara model keratinocyte cell line uncovered novel, potential therapeutic targets and affected pathways

Abstract

Background

To be able to develop effective therapeutics for epidermolysis bullosa simplex (EBS), it is necessary to elucidate the molecular pathomechanisms that give rise to the disease’s characteristic severe skin-blistering phenotype.

Results

Starting with a whole-transcriptome microarray analysis of an EBS Dowling-Meara model cell line (KEB7), we identified 207 genes showing differential expression relative to control keratinocytes. A complementary qRT-PCR study of 156 candidates confirmed 76.58 % of the selected genes to be significantly up-regulated or down-regulated (p-value <0.05) within biological replicates. Our hit list contains previously identified genes involved in epithelial cell proliferation, cell-substrate adhesion, and responses to diverse biological stimuli. In addition, we identified novel candidate genes and potential affected pathways not previously considered as relevant to EBS pathology.

Conclusions

Our results broaden our understanding of the molecular processes dysregulated in EBS.

Background

KRT5 and KRT14 are the main stress-absorbing keratins in basal keratinocytes of the epidermis and related stratified epithelia. These rod-shaped proteins form heterodimeric units that interact to build up the cytoskeletal intermediate filament (IF) network, a resilient yet adaptable scaffold that maintains cellular structural stability and, in turn, normal skin integrity and function. Disruption of IF-organization as a consequence of keratin mutation is the basis of a number of inherited skin fragility syndromes [1, 2]. In the case of epidermolysis bullosa simplex (EBS), which exhibits several clinical variants [3], specific phenotypes can be largely correlated with the positions of missense mutations in structurally sensitive portions of either KRT5 or KRT14 [4]. Aberrant IF organization results in fragile epidermal basal cells that readily lyse following mild mechanical trauma or minor traction, leading to intraepithelial fluid accumulation and recurrent blister formation [5–8]. At the molecular level, this cytoskeletal collapse manifests as aggregates of misfolded keratins, along with activation of stress-response cascades [9, 10]. Detailed elucidation of the underlying pathomechanisms in EBS is an important prerequisite for developing innovative therapeutics; however, relatively few studies have focused on expression-profiling of mutant epidermal cells (summarized in Additional file 1: Table S1). In 2007, Lu et al. [11] reported the expression profile of KRT5^−/− EBS mouse epidermis, concentrating mainly on the regulation of inflammatory cytokines. Liovic et al. [12] investigated the response of the EBS keratinocyte cell lines KEB4 (mild EBS-loc phenotype, KRT14 mutation V270 M) and KEB7 (severe EBS Dowling-Meara (EBS-DM) phenotype, K14 mutation R125P) to hypo-osmotic stress in comparison to the wild-type cell line NEB1, and found dual-specificity phosphatases and their downstream targets ERK and p38 to be differentially regulated in EBS cells. A subsequent study by the same group identified differences in the expression profiles of cell-junction components in EBS versus wild-type cell lines [13]. More recent profiling studies reported aberrant expression of genes involved in keratinization, cell growth, proliferation, the immune response, and fatty acid metabolism in EBS [14]. Because there is limited overlap of the genes identified in those studies, further efforts seem necessary to fully elucidate EBS-relevant genes to better understand EBS pathology. We recently described the gene expression profile of KEB7 cells after applying suppression subtractive hybridization, and found dysregulated genes involved in keratinocyte differentiation, migration and wound healing [15]. Here, we follow-up that analysis with a more expansive expression profiling study, combining whole-transcriptome microarray examination with bioinformatics-assisted functional clustering, complementary qRT-PCR validation, and western blotting analysis of selected hits. Using this approach we were able to verify candidate genes previously described as being differentially expressed in EBS-DM, and to discover other differentially regulated genes not previously implicated in EBS pathology. A graphical summary of our experimental design is shown in Fig. 1.

Fig. 1

Experimental design of our study

Methods

Cell lines and culture conditions

The immortalized keratinocyte cell lines KEB7 (EBS-DM severe phenotype caused by the R125P mutation in KRT14) and NEB1 (wild-type control cell line) [16] were cultured in standard DMEM (HyClone Laboratories, GE Healthcare, South Logan, UT, USA) supplemented with 25 % Ham’s F12 nutrient mixture, 10 % fetal bovine serum (FBS), 1.8 × 10⁻⁴ M adenine, 0.4 µg/ml hydrocortisone, 5 µg/ml transferrin, 2 × 10⁻¹¹ M liothyronine, 5 µg/ml insulin and 10 µg/ml EGF. Unless otherwise stated, all chemicals and reagents were obtained from Life Technologies™ (Karlsruhe, Germany) or Sigma-Aldrich Inc. (Taufkirchen, Germany). Cell lines were cultivated in standard 25 cm² flasks (Techno Plastic Products AG, Trasadingen, Switzerland) at 37 °C under a 5 % CO₂ atmosphere without fibroblast feeder cells.

Total RNA isolation from cultivated KEB7 and NEB1 cells

Total RNA from cells grown to approximately 80 % confluence was isolated using an RNeasy Mini Kit (QiagenGmbH, Hilden, Germany) according to the manufacturer’s instructions. RNA was resuspended in nuclease-free water and quantified spectrophotometrically at 260 nm (DS-11 Spectrophotometer, DeNovix, Wilmington, DE, USA). RNA preparations were considered as suitable only if samples were of sufficient yield, exhibited intact bands on agarose gels, and displayed no spurious peaks or RNA degradation artifacts on the UV absorption spectrum. RNA samples were stored in aliquots at −80 °C.

Whole-transcriptome microarray analysis

Total RNA, isolated from a pool of different passages of KEB7 versus NEB1 cells (p25 and p28–37), was used for probing a high-precision, whole-transcript Human Gene 1.0 ST array (Affymetrix, Santa Clara, CA, USA). The entire experimental pipeline, from sense target labeling to hybridization, washing, array scanning and final raw data capture, was performed at the Center of Excellence for Fluorescent Bioanalysis KFB, University of Regensburg, Germany, using standard Affymetrix protocols, reagents and instrumentation. Four microarrays per cell line were processed to account for technical and/or biological variability. The raw data were checked for quality, background-adjusted, quantile-normalized (by imposing the same distribution of gene signal intensities for each array used under the condition that the expression of most genes is relatively unaltered across mutant and wild-type cells), and statistically analyzed by the linear model for microarray data (LIMMA) method [17] using the Integromics Biomarker Discovery^® platform (IBD, Integromics, Granada, Spain). Assuming experimental deviations/errors of less than a log2-fold change of +1/−1, only candidates showing regulation greater than +2/−2 were subjected to further analysis. The respective statistical parameters were also used for the subsequent enrichment and qRT-PCR studies.

Gene ontology enrichment and functional clustering of differentially expressed genes

Affymetrix IDs of differentially regulated genes were subjected to comprehensive bioinformatical analysis embedded within IBD. Gene ontology (GO) classification was performed to confirm enrichment of microarray-derived genes and to assign them to specific biological themes and functions. To ascertain GO-term enrichment of genes statistically overrepresented in the candidate hit list (log2-fold change of ≥ +1 and ≤ −1, adjusted p-value of <0.05), the software tool calculated a one-sided hypergeometric p-value (identical to the one-tailed version of Fisher’s exact test), so that terms with values below 0.05 can be considered to be significantly enriched.

cDNA synthesis, primer design, and complementary qRT-PCR

Total RNA isolated from one passage of NEB1 and KEB7 cells (each cell line in technical triplicates) was treated with DNase I (Sigma-Aldrich, Taufkirchen, Germany) and then reverse transcribed into cDNA using an iScript™ cDNA Synthesis Kit (Bio-Rad Laboratories Inc., Hercules, CA USA) according to the manufacturer’s instructions. The resultant cDNAs served as template for complementary quantitative real-time PCR (qRT-PCR) in a 96-well plate format using GoTaq^®qPCR Master Mix (Promega Corporation, Madison, WI, USA) and a CFX96™ instrument (Bio-Rad Laboratories). Using Batch Entrez (http://www.ncbi.nlm.nih.gov/sites/batchentrez) and WIBR UTR extractor (http://jura.wi.mit.edu/bioc/tools/utrs/), we extracted the coding sequences (highlighted within the downloaded mRNAs) of the 79 most promising up-regulated and the 79 most promising down-regulated genes uncovered in the microarray analysis to automatically design oligonucleotides with similar properties (uniform GC content, comparable melting temperature, optimal length and similar product size; BATCH Primer 3 program, http://probes.pw.usda.gov/cgi-bin/batchprimer3/batchprimer3.cgi). Forward and reverse primer stocks (100 pmol/µl, see Additional file 2: Table S2) were purchased in 96-well format from Life Technologies GmbH (Karlsruhe, Germany). Differential expression of selected candidate genes was investigated by performing three independent PCR-runs (biological and technical replicates) in an in-house developed and established 96-well qRT-PCR array format. Seven constitutively expressed housekeeping genes (ACTB, ANXA1, B2 M, GAPDH, HPRT1, RPL13A, and TUBB) served as internal controls for data normalization as well as determination of experimental variance. Fold inductions and statistical significance were assessed via RealTimeStatMiner^® from Integromics^® that calculates the relative quantity (RQ) and log2-fold change using ΔΔCt (cycle threshold) and efficiency correction [18]. A parametric LIMMA moderate t-test [17] and a Benjamini–Hochberg false discovery rate (FDR) p-value correction [19], included as standard default parameters of the software (log2-fold change of +1/−1, p-value <0.05), were applied in the data analysis. The meaning of an FDR-adjusted p-value is that, if all genes with a p-value below a threshold of 0.05 are picked as differentially expressed, then the expected proportion of false discoveries in the selected group is controlled to be less than the threshold value, in this case 5 %.

KEB7/NEB1 total protein lysates and Western blotting

For preparation of KEB7 versus NEB1 protein lysates (using three 75-cm² flasks per cell line), the culture medium was removed, the cells were rinsed twice with sterile PBS, and then lysed in RIPA buffer (Radio-Immunoprecipitation Assay buffer; 150 mM NaCl, 1 % IGEPAL^® CA-630, 0.5 % deoxycholate, 0.1 % SDS, 50 mM Tris–HCl, pH 8.0, Sigma-Aldrich) supplemented with 1× Protease Inhibitor Cocktail (GE Healthcare, USA) for 10 min at RT (room temperature) and a further 20 min at −80 °C. Cells were detached with a cell scraper, and the resultant total protein lysates transferred to 1.5-ml tubes and centrifuged at maximum speed for 15 s to remove cellular debris. Protein in the supernatants was photometrically quantified by a standard Bradford assay (Bio-Rad, Hercules, CA, USA), and the samples were frozen at −80 °C for long-term storage. 10 µg of protein lysate per slot were subjected to SDS-PAGE (sodium dodecyl sulfate polyacrylamide gel electrophoresis) and then transferred onto Amersham™ Hybond™ ECL nitrocellulose membranes (GE Healthcare). After blocking of the membrane in 3 % nonfat milk in Tris-buffered saline (TBS, pH 8.0, Sigma-Aldrich) containing 0.05 % Tween-20 for 1 h at RT, target-specific primary antibodies diluted 1:1000 in TBS with 0.05 % Tween-20 were added overnight at 4 °C with gentle agitation (Wnt-5a antibody, sc-365370, mouse monoclonal, Santa Cruz Biotechnology, Inc.; Cytokeratin 19 antibody, 11-120-C020, mouse monoclonal, EXBIO Praha a.s.; SOX2 antibody, D1951-1NB-1-1R1M1M3/1C16_110120, mouse monoclonal, ABMART; GAPDH antibody, TA302944-100, goat polyclonal, OriGene; antibodies were pre-tested for adequate performance in different human cell lysates; data not shown). The next day the membranes were incubated in horseradish-peroxidase (HRP)-labeled secondary antibodies (1:10,000, polyclonal rabbit anti-mouse IgG/HRP, polyclonal rabbit anti-goat IgG/HRP, Dako Denmark A/S, Glostrup, Denmark) for 1 h at RT. Blots were developed using Immobilon™ Western Chemiluminescent HRP Substrate (Merck/Millipore, Darmstadt, Germany) and signals detected with a ChemiDoc XRS system (Bio-Rad Laboratories). Band intensities were densitometrically quantified by using the image-processing program ImageJ (http://rsb.info.nih.gov/ij/).

Results

Global transcriptional profiling uncovers 207 differentially expressed genes in the EBS-DM cell line KEB7

Transcriptional profiling using the Affymetrix Human Gene 1.0 ST gene chip platform offers coverage of nearly the entire transcriptome, defined by 28,869 annotated open reading frames/genes. We applied this platform to study the expression pattern of genes influenced by the cytoskeletal collapse caused by the R125P mutation in KRT14 in KEB7 cells. Table 1 summarizes the results of our microarray profiling (four technical replicates each) using a pool of different passages (biological replicates) of KEB7 cells relative to NEB1 wild-type control cells. A LIMMA statistical test comparing mutant versus wild-type cells identified 106 up-regulated and 101 down-regulated genes with significantly altered expression levels (log2-fold change ≥+1 and ≤−1, adjusted p-value <0.05) out of 382 differentially expressed genes (failed at log2-fold change and p-value set).

Table 1 Genes regulated in the EBS-DM model keratinocyte cell line KEB7 identified by microarray analysis

Bioinformatic enrichment and clustering suggest previously undescribed pathways and functions are affected in EBS-DM

We next assigned cellular themes to the differentially expressed genes by using IBD software. Functional annotation clustering allowed us to group highly related annotation terms into enriched functional categories, as shown in Table 2. The enrichment p-values of each GO term are presented as the ordered p-value of a one-sided hypergeometric test. Significance of enrichment is indicated by a p-value of <0.05. Apart from already well-known clusters such as the regulation of epithelial cell proliferation [20], cell adhesion and response to stress [12, 13], response to retinoic acid and UV [21–23], immune response [11] or fatty acid metabolism [14], our analysis also identified previously undescribed functional classes potentially implicated in EBS-DM. The latter include genes involved in morphogenesis, and genes that regulate BMP- and canonical Wnt-signaling. In addition, two enrichment clusters associated with protein tyrosine phosphatases and symporter activities, and one further cellular component cluster linked to neurofilaments, were uncovered.

Table 2 Gene ontology (GO) functional annotation clustering

A complementary qRT-PCR array approach confirms differential regulation of the majority of selected candidates

Typically in microarray investigations, just a few randomly selected candidates are re-analyzed by qRT-PCR analysis for confirmation of the first hits, but this can lead to a significant number of undetected false positives. In the present study we used an array-based qRT-PCR analysis to investigate 158 of the 207 (76 %) microarray hits. To reduce the number of transcripts to be able to fit our in-house 96-well qRT-PCR array arrangement, we examined 79 randomly selected up-regulated and 79 down-regulated genes (plus controls) in three independent PCR runs (biological and technical replicates). Furthermore, we based our fold-change cut-off values on the average of seven housekeeping genes (rather than on two, as has been customary in the literature [24, 25], so at to improve the reliability of the data obtained). Out of the 158 genes analyzed, significant (p-value <0.05) differences in expression (fold change set ≥+2 for up-regulated and ≤−2 for down-regulated genes) of 121 genes (~76.58 %) was confirmed, as highlighted in Table 1. Eleven genes (~6.96 %) turned out not to be significantly modulated, and 26 genes (~16.45 %) completely failed in all independent PCR runs, probably due to poorly designed primers; this would potentially give an underestimate of the number of reproducible hits.

Western blot analysis confirms differential expression at the protein level of 3 selected candidates

Western blot analysis using crude lysates of KEB7 and NEB1 cells confirmed differential regulation of some of the selected candidates at the protein level after normalization with glyceraldehyde-3-phosphate-dehydrogenase (GAPDH), as shown in Fig. 2. Sex determining region Y-box 2 (SOX2, a log2-fold change of −1.71 by microarray analysis, no data by qRT-PCR) and cytokeratin 19 (KRT19, log2-fold change of −5.73 by microarray analysis, 1633-fold down-regulation by qRT-PCR) were both found to be significantly down-regulated in KEB7 cells, whereas wingless-type MMTV integration site family member 5a (WNT5a) was found to be significantly up-regulated (log2-fold change of 1.94 by microarray analysis, 6.9-fold up-regulation by qRT-PCR). Based on densitometric evaluation of the protein blots, SOX2 is 25-fold down-regulated, KRT19 61.7-fold down-regulated, and WNT5a 123-fold up-regulated in KEB7 cells.

Fig. 2

Western blot analysis for confirmation of differential expression at the protein level of 3 selected candidates—KRT19, WNT5a and SOX2. GAPDH occurred as normalization protein

Discussion

Global transcriptional profiling uncovers novel differentially expressed genes as well as affected pathways in the EBS-DM cell line KEB7

The present study employed global high-throughput microarray analysis together with bioinformatics-assisted functional clustering and enrichment to identify 207 differentially expressed genes in the EBS-DM cell line KEB7. A subset of 158 genes was subjected to further investigation by qRT-PCR, of which three quarters (76.58 %) could be validated as being differentially regulated. Parenthetically, the used human gene 1.0 ST Array from Affymetrix also has the ability to detect alternative mRNA isoforms. The mRNA source of the array transcripts arise from the reference sequence database (RefSeq). Each of these transcripts is read out by 26 different oligonucleotides (25mere), which are spread about the whole length of the transcription unit. Hence, a averaged probe set consists of 26 independent measure points per transcript, which are summated in only one signal. Because we have not taken into account this feature of the array to detect these alternative mRNA isoforms or splice variants we designed for qPCR validation only primer pairs for one transcript variant. Therefore, the number of our reproducible hits could potentially be underestimated.

In addition to verifying genes from previous studies (and thereby supporting their potential role in EBS pathology), a novel contribution of the present work is its uncovering of genes involved in the BMP-signaling pathway, the canonical Wnt-receptor signaling pathway, and the response to retinoic acid, suggesting that dysregulation of these processes as well may play a role in disease pathology. It is well demonstrated that keratinocyte differentiation is closely linked to hormonal action, particularly to calcitriol- [26, 27] and retinoic acid-signaling [21, 23]. In the present study, we identified four differentially regulated genes in KEB7 cells having roles in retinoic acid signaling: three genes, HSD17B2 (hydroxysteroid (17-β) dehydrogenase 2), SOX2 and MEST (mesoderm specific transcript), are down-regulated, and DKK1 (dickkopf Wnt signaling pathway inhibitor 1) is up-regulated. These potential disturbances to retinoic acid signaling may at least partially account for the observations by Wagner et al. [15] of aberrant differentiation processes in KEB7 cells. Generally, retinoids are considered to play a role in normalization of keratinocyte differentiation by down-regulating desmosomal proteins, exerting anti-proliferative effects, and regulating lipid-synthesis, growth factors, and cytokines [22, 23, 28, 29]. An interesting point in this regard is the ability of retinoic acid to act as an antagonist of Jun N-terminal kinase (JNK)-signaling, mediated by activating protein 1 (AP1) [30], since JNK/mitogen activated protein kinase (MAPK) signal transduction has previously been shown to be dysregulated in EBS-DM model keratinocytes [9, 31, 32]. Furthermore, we recently advanced a positive feedback model, suggesting that mutant KRT14 leads to an increase in AP1-dependent expression of KRT14, which in turn amplifies the level of aberrant JNK stress-signaling [33]. Concerning the canonical Wnt-receptor signaling pathway, enrichment clustering applied in our present study identified four significantly dysregulated genes: GPC3 (glypican3) and SOX2 are both down-regulated, and WNT5a and PRICKLE1 (prickle homolog 1) are up-regulated. Interestingly, WNT5a is one of the so-called non-canonical Wnt ligands. During normal development, WNT5a is secreted and directs the migration of target cells along concentration gradients. However, deregulated WNT5a signaling facilitates invasion by multiple tumor types into contiguous tissues. EB per se is often associated with carcinoma development, mostly in EB dystrophicans (DEB) and junctionalis (JEB) [34–36], but rare cases of carcinoma in EBS have also been described [36]. To date, the roles of WNT5a in cutaneous squamous cell carcinoma (SCC) and basal cell carcinoma (BCC) as well as the effect of WNT5a on keratinocyte migration have not been fully investigated, although Pourreyron et al. [37] recently demonstrated up-regulation of WNT5a in SCC/BCC and its localization to the leading edge of tumors as well as in tumor-associated fibroblasts. Regarding BMP signaling, enrichment clustering showed positive regulation of this pathway in EBS by highlighting three genes, GPC3, MSX2 (Msh homeobox 2) and ILK (Integrin-linked kinase), as being significantly down-regulated. BMPs are secreted signaling polypeptides belonging to the transforming growth factor-β (TGF-β) superfamily that act as important multifactor players in the development of vertebrate skin (and appendages). BMPs thereby not only regulate hair follicle (HF) morphogenesis and keratinocyte differentiation in both the HF and epidermis, but are, along with other interacting growth factor family members, also involved in normal postnatal tissue remodeling and homeostasis [38, 39]. Furthermore, BMP signaling acts as a potent tumor suppressor in the skin, inhibiting mainly epidermal- and HF-derived tumor formation [38, 40]. Consequently, dysregulation of this pathway can lead to abnormal skin development and tumor formation. Progress in this research area could be interesting to get a better knowledge of the disease pathology.

KRT19—a novel biomarker in EBS-DM?

A further example from our list of dysregulated genes is the smallest known acidic keratin, KRT19 (significantly down-regulated in our microarray analysis, independently confirmed by both qRT-PCR and western blotting). In our previous study, in which we applied a similar whole-transcriptome microarray methodology to calcitriol-stimulated human primary keratinocytes, we [41] found a significantly increased expression level of KRT19 after treatment with the indicated secosteroid. Hence, irrespective of the exact role of KRT19 in the disease mechanism, it could constitute another useful biomarker in EBS.

Conclusions

In summary, our study employed global high-throughput microarray analysis plus bioinformatics-assisted functional clustering and enrichment to uncover 207 differentially expressed genes in the EBS-DM cell line. Thereof, 158 genes were subjected to further investigation by a complementary qRT-PCR analysis, of which 76.58 % could be validated as being differentially regulated. We here present a pool of novel candidate genes and potential affected pathways, which may play a role in the disease pathology.