Introduction
Basal cell carcinoma (BCC) is the most common malignancy among Caucasians, with a rising estimated yearly incidence of 2.75 million cases worldwide (1). BCC is more common among elderly men, with a peak incidence after the age of eighty. Individuals with a fair phenotype, including red or blonde hair and light eyes, are particularly at risk. While BCC can develop on any skin surface, the majority of lesions appear on the head and neck. Based on their histological and clinical features, BCC can be classified into one of several types including nodular, superficial, morpheaform, nevoid, and pigmented (2). Nodular BCC is the most common subtype, comprising 21% of all cases (3). Although BCC is slow-growing and rarely metastasize, it is locally invasive and thus may cause extensive damage to surrounding tissue. Despite effective treatment via local excision, tumor recurrence is relatively common at 1–10% (4).
Despite its high prevalence, the etiopathogenesis of BCC is still unclear. Previous studies have indicated a multifactorial, polygenic basis for disease. The current model for BCC pathogenesis maintains that UV radiation causes DNA damage in exposed cells. If this damage goes unrepaired, the resulting oncogene-activating or tumor suppressor-inactivating mutations allow cells to bypass cell cycle regulation and thus undergo uncontrolled proliferation. The role of UV damage in BCC pathogenesis is further indicated by the predominant location of BCC on sun-exposed surfaces, as well as the presence of ‘UV signature’ mutations (i.e., T → C transversions) in affected cells (5,6). Other exposures that can predispose to carcinogenesis are those to arsenicals, polyaromatic hydrocarbons, immunosuppression, and psoralen therapy (2).
Studies have identified several somatic mutations associated with the development of BCC (2,5,7–9). Of note, the presence of mutations in the Patched-1 (PTCH-1) gene has helped elucidate the role of the sonic hedgehog (SHH) signaling pathway in pathogenesis. PTCH-1 is an inhibitor of the protein Smoothened (SMO), whose role is to activate the transcription of cell cycle regulators like WNT5A. When PTCH-1 activity is lost, SMO is constitutively activated, allowing uncontrolled cell proliferation to take place (5). A loss of function mutation in P53 has also frequently been associated with BCC as well as many other cancers (6,10). P53 encodes the protein p53, an established tumor suppressor. This protein causes cell cycle arrest in the presence of DNA damage, thereby preventing the replication of damaged genetic material and allowing either DNA repair or apoptosis to take place.
Genome-wide association studies (GWAS) have been instrumental in identifying several loci that confer susceptibility to BCC. These studies report loci associated with BCC and pigmentation genes (SLC45A2, TYR, MC1R, ASIP) as well as loci associated with BCC alone (PADI6, RHOU, KLF14, KRT5, TERT/CLPTM1L). These results support the conclusion that both pigmentation-independent and pigmentation-dependent pathways exist in the development of BCC. Although genomic studies and linkage analyses provide a framework for identifying putative loci, they do not address the gene expression that underlies disease pathogenesis.
In this study, we employed microarray analysis to determine a molecular profile for BCC which was then used to i) classify samples based on phenotype via hierarchical clustering methods, ii) identify significantly enriched pathways important to BCC pathogenesis, and iii) integrate genetic information with our transcriptional data in order to identify potential genetic risk factors. Specifically, we obtained a list of differentially expressed genes (DEGs) from a comparison of lesional vs. site-matched non-lesional skin samples from patients with a confirmed diagnosis of nodular basal cell carcinoma. Pathway analysis of the resulting DEGs identified multiple dysregulated functional pathways, including those involved in PPAR-γ signaling, TGF-β signaling, and lipid metabolism. A comparison of our list of DEGs to molecular profiles in published studies identified several overlapping genes and pathways of interest in BCC. We further compared the chromosomal locations of our list of DEGs with reported genomic susceptibility loci to focus the search for genes of pathogenetic significance. Furthermore, our analysis identified potential transcriptional ‘hot spots’ in which there is an enhanced correlation of significantly altered gene expression at particular chromosomal locations, areas which may be of particular interest for future genetic studies. By integrating transcriptional data with genomic information, our study reveals potential susceptibility loci/regions associated with BCC in terms of altered gene expression. Detailed characterization of the genetic and transcriptional alteration associated with BCC development may lead to novel therapeutic modalities based on specific genomic and/or transcriptional targets.
Results
Hierarchical clustering separates samples by disease status
A total of 8 skin biopsies from patients with nodular BCC, lesional (n=4) and non-lesional (n=4), were analysed using Affymetrix Human Genome U95A2 microarray chips (Affymetrix) to generate gene expression profiles (see Table I for details on samples and patients). We performed the comparison of expression profiles between the tumor and normal groups. At the false discovery rate of 0.1, we identified a total of 331 differentially expressed genes (DEGs) with at least 2-fold expression changes. 144 DEGs are upregulated (fold changes ranging from +2.0 to +53.1) in the tumor group while 187 genes are down-regulated (fold changes ranging from −2.0 to −32.7). Hierarchical clustering of obtained DEGs was performed, which separates the 8 samples into their corresponding phenotype groups (Fig. 1).
Functional analysis reveals dysregulation of genes involved in multiple pathways
To explore the key biological processes altered in the tumor vs non-tumor control samples, we performed enrichment tests to identify the significantly overrepresented canonical pathways among the differentially expressed genes. Functional annotation and pathway analysis were performed using the database for annotation, visualization, and integrated discovery (DAVID) and Pubmed literature searches (Table II). The pathways enriched in the differentially expressed genes include cell-cell interactions such as focal adhesion (19 genes) and ECM-receptor interaction (16 genes), PPAR signaling pathway (11 genes) and TGF-β signaling pathway (7 genes), terpenoid backbone biosynthesis (6 genes), and fatty acid metabolism (5 genes). Our analysis revealed, with the exception of one gene (MMP-1), the downregulation of DEGs falling within the PPAR-γ signaling pathway. The downregulation was also seen in DEGs involved in fatty acid metabolism, and terpenoid backbone biosynthesis. Additionally, the DEGs pertaining to focal adhesion and cell-cell interactions are mostly upregulated, with the exception of one gene (ITGA7) (Fig. 2).
Comparison of DEGs across microarray studies reveals overlapping genes/pathways of interest
To evaluate potential consensus genes associated with BCC pathogenesis, we compared our list of DEGs to four previously published microarray studies regarding gene expression in BCC (17–20). In our analysis, we excluded Yu et al due to significant methodological differences. Specifically, the authors compared molecular profiles between different BCC subtypes and not between tumor (without subtype distinction) and normal skin. We first examined data from Howell et al, as their study used site-matched non-lesional samples as controls (Table III, second left-most column). Comparing our data to a similarly conducted study allowed us to minimize the presence of potentially confounding experimental design and technical variances. Twenty-six DEGs were found to overlap between our study and Howell et al; 8 genes were upregulated in the same direction (MDK, LUM, COL4A1, CDH11, DUSP10, COL5A2, STAT1, and SDC2), 14 genes were downregulated in the same direction (NR4A1, CYB5A, APOC1, DHCR24, PLA2G2A, FDPS, PPARG, ADH1B, HMGCR, DUSP1, PLA2G7, LPL, FABP4, and ALDH1A1) and 4 genes were differentially expressed in opposite directions (UBE2D1, KRT7, KRT18 and DAPK1). Pathway analysis revealed genes that were differentially expressed in processes such as PPAR-γ signaling, cell-cell interaction, terpenoid backbone biosynthesis, and MAPK signaling.
We also compared our data more broadly with combined data from Howell et al(18) as well as with O’Driscoll et al(19), and Asplund et al(17), disregarding certain methodological differences. A total of 1842 DEGs were compared and 193 genes were found to overlap with at least one other study (Table III). 186 DEGs were dysregulated in the same direction (100 upregulated and 86 downregulated) while 7 genes were differentially expressed in opposite directions (FOSB, CYR61, DICER1, UBE2D1, KRT7, KRT18 and DAPK1). Pathway analysis of these DEGs revealed a dysregulation of the genes involved in focal adhesion, extracellular matrix-receptor interaction, terpenoid backbone biosynthesis, and steroid biosynthesis, which overlapped with a large number of pathways derived from analysis of our DEGs. Thirteen DEGs overlapped across three studies, including our own (Table IV). Functional annotation and pathway analysis revealed that 5 of the 13 DEGs were involved in either cell-cell interaction or terpenoid backbone biosynthesis.
DEGs and transcriptional ‘hot spots’ map to several genetic susceptibility loci/regions
We have mapped chromosomal locations for the top 20 upregulated and downregulated DEGs identified in our study (Table V). We then compared chromosomal locations of DEGs with putative BCC susceptibility loci previously reported in genome-wide association and linkage studies, as well as regions where somatic mutations, determined either in human subjects or murine models, have been implicated in the pathophysiology of BCC (Table VI). A total of 62 DEGs, 25 upregulated and 37 downregulated, were mapped to these regions.
Next, we examined our gene expression data to identify chromosomes with a significant enrichment of DEGs. A statistically significant over-representation of DEGs was found on chromosome 5 (odds ratio=1.65, P=0.025) (Fig. 3). The location of DEGs within chromosome 5 includes 5p13–15.33, 5q11.2–14.3, 5q22.1–23.3, and 5q31–35.3 (Fig. 4), which contain 22 DEGs (12 upregulated, 10 downregulated). The transcriptional ‘hot spots’ on 5p13 and 5p15 were also identified in recent BCC genome-wide association studies. The 5q11.2–14.3, the 5q22.1–23.3 and the 5q31–35.3 regions represent novel chromosomal locations with potential relevance to BCC that may warrant further genetic investigation.
Discussion
Although BCC has been investigated at both the genetic and transcriptional levels, many details related to pathogenetic events remain unknown. To help elucidate the links between the genetic factors and altered gene expression that impact tumor development, we performed gene expression analysis to define a BCC associated transcriptional profile.
Our analysis was conducted using 8 samples, with 4 lesional samples and 4 site-matched non-lesional controls. We utilized a site-matched, pair-wise study design to eliminate variance due to differences in gene expression between individuals. This approach allowed us to focus directly on transcriptional alterations between tumor and non-tumor tissues without the aforementioned confounding variables. Further studies investigating differences between non-lesional skin in patients with BCC and skin of healthy individuals may provide additional insight into underlying genetic contributions to BCC. While differences in sample preparation, sample type, microarray platform, and analytical software present complicating factors in comparative analysis, our results do overlap with several DEGs reported in three previous microarray studies (Tables III and IV).
Hierarchical clustering of obtained DEGs revealed that samples could be distinguished by disease status, with disease and control samples separating into discrete groups. The establishment of a BCC tumor expression signature may be useful for the development of molecular diagnostic modalities in BCC, and extended transcriptional profiles could allow classification of newly diagnosed BCC into specific subtypes (i.e. nodular, morpheaform, and superficial) as an addition to standard clinical and pathological criteria. In the future, it may be possible to classify patients into subgroups according to predicted therapeutic efficacy or risk of recurrence, thereby improving patient outcomes.
Functional analysis of our data revealed transcriptional dysregulation in multiple pathways affecting PPAR-γ signaling, lipid metabolism, TGF-β signaling, cell-cell interactions, as well as many others, with intriguing implications regarding cancer pathogenesis and localizing potential therapeutic targets. For example, recent research into the PPAR-γ signaling has elucidated its role in a variety of cellular processes. PPAR-γ is a receptor whose ligands include steroid, hormones, and retinoids (21–23). Once activated, PPAR-γ dimerizes with the retinoid X receptor to activate transcription of genes involved in lipid metabolism and differentiation. In our study, there was a downregulation of PPAR-γ as well as its transcriptional target genes (ADIPOQ, FABP4, PLIN1, LPL, ACS, NR4A1, FADS2, HMGCS1), indicating aberrant PPAR-γ signaling in our samples. Our analysis also revealed transcriptional dysregulation in pathways regarding lipid metabolism, unsaturated fatty acid biosynthesis, steroid biosynthesis as well as terpenoid backbone biosynthesis. Whether these processes are a cause or consequence of PPAR-γ signaling dysregulation remains to be determined.
Interestingly, the PPAR-γ signaling pathway has been implicated in a variety of cancers, including, but not limited to, bladder cancers, colon cancers, squamous cell carcinomas and melanomas (21,22,24–26). PPAR-γ has been shown to be down-regulated in certain tumors (27,28); however, other studies have shown an upregulation of PPAR-γ in other malignancies (29–31). These discrepancies warrant further investigation on the role of PPAR-γ in specific tumors. PPAR-γ activation has been implicated in reducing cell proliferation and/or inducing apoptosis in a wide array of cancer cell lines. In vitro studies on lung cancer cell lines have shown that PPAR-γ activation acts to reduce cell growth by promoting differentiation (32). Some transcriptional target genes have been identified as prognostic factors for certain tumors. Of note, fatty acid binding protein 4 (FABP4) expression has been shown to correlate with bladder cancer progression (33). FABP4 plays a role in signal transduction, affects glucose and lipid metabolism, and potentiates apoptosis (34). Decreased levels of FABP4 denoted a worse prognosis in patients with bladder tumors. Our study revealed a decreased expression of FABP4. Future studies may indicate if transcriptional levels of FABP4 can be used as an indicator of cancer progression in basal cell carcinoma.
Although transcriptional dysregulation of the PPAR signaling pathway has been established in a variety of tumors, our study represents aberrant signaling in basal cell carcinomas. The PPAR-γ signaling pathway represents an intriguing therapeutic target, as PPAR-γ activation via pharmaceuticals has been used in studies to investigate cancer treatment. In vitro use of thiazolidinediones (TZD), such as rosiglitazone and troglitazone, demonstrated anti-proliferative, pro-apoptotic, and differentiation-promoting effects (21,26). Future experiments using PPAR-γ agonists on mouse models as well as BCC cell lines may elucidate more clues on the pathogenesis of the disease and help develop novel therapeutic options for patients with basal cell carcinoma.
Our data confirm certain well-established genetic mechanisms underlying BCC pathogenesis. The sonic hedgehog (Shh) signaling pathway is often cited as an important player in disease development (2,5,7,8). Several genes in the Shh pathway can be found amongst the DEGs in our list. Our analysis revealed the upregulation of PTCH-1, a primary mediator of the Shh pathway. When not bound by Shh, PTCH-1 inhibits Smoothened, thus preventing signal transduction and transcription of downstream target genes. Among these target genes are PTCH-1 for negative feedback, GLI1 for positive feedback, WNT5A, a gene involved in differentiation, and MYCN, a gene associated with the development of neuroblastomas. When inhibition of Smoothened is lost due to mutation, transcription of target genes can occur constitutively and therefore promote disease initiation and progression (35). Given that our results reveal upregulation of PTCH1, WNT5A, and MYCN, this constitutive activation seems likely. Additionally, several DEGs share the same chromosomal locations as established genes of interest in BCC pathogenesis such as PTCH1, SUFU, PTCH2, and Gli2 (Table VI), thus reinforcing their putative role as genetic susceptibility loci.
P53, a tumor suppressor protein responsible for cell cycle arrest in the presence of DNA damage has been implicated in a multitude of cancers, including BCC (9,10,36,37). Our results show that four downregulated DEGs share the same chromosomal location as p53 (Table VI), indicating that 17p13 may in fact be an important susceptibility locus for BCC. Additionally, two other DEGs confirm the importance of p53 in BCC pathogenesis. The first, CDKN1A, encodes p21, a cyclin dependent kinase inhibitor and major mediator of the p53 pathway. Because p21 is tightly regulated by p53, the downregulation of CDKN1A may indicate the presence of a mutation in the tumor suppressor gene (38). P53 also plays an important role in the regulation of MMP1, a degradative enzyme family member that breaks down the extracellular matrix during tissue development, remodeling, and repair (37). In our list of DEGs, MMP1 was upregulated. Given that p53 downregulates MMP1, this may again reflect a gene mutation. MMP1 is also important in tumor progression through its role in the stimulation of tumor-induced angiogenesis and local tissue invasion. The concomitant downregulation of TIMP4, a family member of MMP inhibitors whose down-regulation has been associated with excessive ECM degradation (39), may also contribute to tumor growth and invasion.
Several physiological mechanisms appear to be activated in order to counteract the pathological state in BCC. Our analysis reveals an upregulation of DAPK1 (DAPK1), a candidate tumor suppressor whose mechanism includes inhibition of ERK. This upregulation affects both the Ras-MAPK and TGF-β pathways that may support tumor suppressive changes. The Ras-MAPK pathway plays an important role in many cellular functions, including proliferation, differentiation, migration and cell survival (40). Constitutive activation of this pathway, either via mutation or dysregulation, is associated with many cancers, including melanoma, breast, pancreatic, head and neck, and colon cancers (41–45). In this tightly regulated pathway, the inhibition of ERK by DAPK1 results in the downregulation of FOS and MYC, two oncogenes associated with uncontrolled proliferation and thus tumorigenesis. The TGF-β pathway is similarly involved in proliferation, differentiation, growth and cell death (46). Inhibition of ERK in this pathway allows for the upregulation of the transcription factor E2F5 (E2F5). E2F5 has an established role in the inhibition of MYC (47), therefore its upregulation may be another mechanism for tumor suppression. Further study on the role of E2F5 in BCC is warranted given recent evidence suggesting that E2F5 may contribute to tumorigenesis (48,49).
A recently published meta-analysis from our lab revealed overlapping DEGs across no more than 2 BCC microarray studies (50). We extended this analysis here with the most currently available literature and discovered that 13 DEGs overlapped across 3 studies (Table IV). Functional annotation of these genes revealed transcriptional dysregulation of processes involved in lipid/steroid metabolism and of components of the extracellular matrix. The dysregulated genes involved in lipid metabolism may point to underlying aberrations leading to BCC tumorigenesis and represent potential biomarkers for the development of BCC. To understand if transcriptional dysregulation of lipid metabolism can act as a diagnostic indicator, future studies may include linking fine-tuned time-course measurements with changes in gene expression. Also, we noted an upregulation in extracellular matrix-related genes. Previous studies have indicated that BCC samples that were less invasive typically contained a dense matrix surrounding the cells (51). It has been postulated that this stroma precludes cellular proliferation and tumor metastasis by introducing a physical barrier to migrating cells and helps explain the slow-growing properties of BCC.
Numerous genome-wide association studies have identified putative susceptibility loci that are associated with basal cell carcinoma (52–62). However, consensus risk loci have not been established and these studies do not shed light on the causal relationships between genes and phenotype. Here, we combined information from established genetic linkage studies on BCC susceptibility with gene expression data from our microarray analysis to draw insights on the development of basal cell carcinoma. We expect that subsets of genes differentially expressed in BCC are due to genetic alterations. Thus, disease associated DEGs are likely to represent an enriched pool of candidate risk genes. Moreover, we found 26 BCC associated DEGs that mapped to eight previously reported BCC susceptibility identified loci (Table VI). Four DEGs mapped to locations previously associated with BCC and pigmentation genes related to eye color, hair color, and/or skin color. The remaining DEGs mapped to putative susceptibility loci that were associated with BCC only and not pigmentation. It should be noted that none of the 26 DEGs that mapped to loci were related to pigmentation alone, indicating that any primary genetic associations are likely to be linked to tumor development.
A recent GWAS identified KRT5 as a gene of interest associated with BCC at the 12q11–13 locus (61). Our analysis revealed 9 downregulated DEGs that mapped to the 12q13 region. This pool of 9 DEGs at 12q13 warrant further study to investigate which of these represent true risk loci. Thus, the strategy of merging genetics and transcriptional datasets can be leveraged to refine the search for susceptibility loci, particularly those with functional consequence.
Our study identified chromosome 5 as a chromosome with significantly enriched DEGs. We further found four regions on chromosome 5 that might serve as transcriptional ‘hot spots’ for BCC. In particular, the 5p13–15.33 region overlapped with two susceptibility loci (5p13 and 5p15) that were previously identified in GWAS and linkage studies. Our study also revealed novel chromosomal locations at the 5q11.2–14.3, 5q22.1–23.3 and 5q31–35.3 regions. 5q11 was recently identified in a GWAS as a susceptibility locus in patients with esophageal squamous cell carcinoma (63). The 5q35 region was previously associated with prostate cancer development, however, the causal genes at play were not identified (64,65). In genome-wide association and linkage studies, 5q21.1 and 5q31–33 loci were associated with atopic dermatitis (66,67). These results underscore the potential pathogenetic significance of the identified chromosomal locations, suggesting regions that may be prioritized for rigorous genetic association studies.
We expect our approach of integrating available genetic information with transcriptional data will facilitate future investigations to pinpoint susceptibility loci with greater precision to better illuminate the causative links between genetic alteration, transcriptional dysregulation, and disease initiation and progression in BCC. Taken together, such information could be used to improve on current diagnostic and prognostic modalities, and further our understanding of disease mechanisms in order to develop enriched targets for therapy.