Main

Skin cancers are by far the most common malignancies in fair-skinned populations, with an incidence now reaching epidemic proportions (Madan et al, 2010; Rogers et al, 2010). There are ∼700 000 cases of cutaneous squamous cell carcinoma (cSCC) diagnosed each year in the United States and the frequency is rising worldwide (de Vries et al, 2005; Skin Cancer Foundation, 2010). While the overall survival of patients with cSCC is high, these tumours are associated with considerable morbidity and pose a substantial financial burden for health-care systems. When cSCC metastasises, the prognosis is poor, with a 5-year survival rate of <25% (Rowe et al, 1992; Kraus et al, 1998). Immunosuppressed patients such as organ transplant recipients (OTRs) and individuals with chronic lymphocytic leukaemia are at significantly increased risk for cSCC and are prone to multiple and often aggressive tumours (Euvrard et al, 2003; Harwood et al, 2006, 2013).

An estimated 65% of cSCC arise from precursor lesions termed actinic keratoses (AKs) (Criscione et al, 2009). Actinic keratosis is associated with epidermal atypia and is thought to represent a continuum of progression from dysplastic keratinocytes to cSCC. Actinic keratoses affect up to 18% of people in the United Kingdom over 60 years of age, rising to 64% in Australia (Frost et al, 2000; Memon et al, 2000). They are the strongest independent risk factor for cSCC development, and often present as multiple lesions, producing ‘field cancerisation’ if confluent. The precise rate of progression is unknown, however, a study from the United States that prospectively followed 169 patients with a total of 7784 AK estimated the risk of progression for an individual lesion was 2.57% at 4 years (Criscione et al, 2009). Evidence for progression of AK to cSCC is also provided by genetic studies that report AKs have a similar karyotypic profile to cSCC, but display a reduced degree of complexity, consistent with an earlier stage of tumour development (Ashton et al, 2003).

Despite the frequency of cSCC, their underlying molecular pathogenesis is poorly characterised, especially the changes involved in progression from AK to cSCC. Recent exome-level sequencing of cSCC revealed a huge mutational burden of ∼1 per 30 000 base pairs of coding sequence (Durinck et al, 2011). This makes cSCC the most highly mutated human malignancy, compounding the difficulty in defining ‘driver’ molecular events underlying their development. Evidence to date supports important tumour suppressor roles for TP53 and NOTCH, but the contribution of additional genes and pathways is unclear (Durinck et al, 2011; Wang et al, 2011). Earlier studies demonstrated considerable karyotypic complexity in cSCC, with frequent gains on chromosomes 3q, 8q and 20q and losses of 3p, 4p, 9p, 13q, 17p and 17q (Quinn et al, 1994; Popp et al, 2002; Ashton et al, 2003; Clausen et al, 2006; Purdie et al, 2007, 2009). Our previous genome-wide analysis using single-nucleotide polymorphism (SNP) microarrays revealed that well-differentiated (WD) tumours display a distinct genetic profile from moderately (MD) and poorly differentiated (PD) tumours, indicating that these may represent a separate subgroup of cSCC (Purdie et al, 2009). However, current published expression microarray studies have failed to corroborate this at the level of gene expression and have yielded little consensus on the genes that are differentially expressed in cSCC (Van Haren et al, 2009). This can be largely attributed to small sample sizes, non-microdissected tumour specimens and the absence of correlation with cSCC grade among the tumours analysed (Haider et al, 2006; Kathpalia et al, 2006; Nindl et al, 2006).

The purpose of our study was to identify potential ‘drivers’ in the transition from AK to cSCC by comparing the transcriptome of tumour cells isolated by laser capture microdissection from cSCC, with that of dysplastic keratinocytes from AKs. We have identified a series of differentially expressed genes (DEGs) in cSCC and have demonstrated key biological processes that distinguish cSCC from AK. Furthermore, an enrichment of DEGs has been demonstrated in key pathways such as the mitogen activated protein kinase (MAPK) pathway, implicating this pathway in the evolution of cSCC from precursor lesions.

Materials and methods

Tumour samples

Fresh-frozen biopsies of normal skin, AK and cSCC were obtained at the time of surgical excision and immediately snap-frozen in liquid nitrogen. All patients provided informed consent in accordance with ethical approval from the East London and City Health Authority local ethics committee. For microarray analysis, patient material was preferentially included where matched AK and cSCC lesions were available from the same patient. Subsequent samples were chosen to represent a range of histological diagnoses from both immunosuppressed and immunocompetent patients. In total, 10 AKs and 30 primary cSCC (9 matched to AK) were analysed. The cSCC comprised 15 WD, 9 MD (including 1 cSCC admixed with an additional WD component, W-MD) and 6 PD (including 4 cSCC admixed with an MD component, M-PD). The samples were derived from 10 immunocompetent patients, 20 OTRs and 1 patient with chronic lymphocytic leukaemia, reflecting the spectrum of patients presenting with cSCC to our institution. Patient and sample characteristics are detailed in Table 1.

Table 1 Sample characteristics
Full size table

Expression microarray analysis

Laser capture microdissection and RNA extraction were performed as previously described (Lambert et al, 2012), such that all samples were estimated to be a minimum 90% enrichment for tumour or dysplastic cells. RNA quality and concentration was measured using the Agilent Bioanalyser (Agilent, Berkshire, UK). Ten nanograms of RNA was used as a template for whole transcriptome amplification and cDNA synthesis using the NuGen WT-ovation Pico RNA Amplification System according to the manufacturer’s instructions (NuGen, San Carlos, CA, USA). Five micrograms of amplified cDNA was labelled with the FL-ovation Biotin Kit (NuGen) and hybridised to the Affymetrix HGU133 Plus 2.0 microarrays (Affymetrix, Santa Clara, CA, USA), comprising 54 675 features. All procedures were performed according to the manufacturer’s recommended protocol. Raw data for all samples are available through the Gene Expression Omnibus (GEO: http://www.ncbi.nlm.nih.gov/geo/), accession number GSE45216.

Statistical analysis

Quality control analysis was performed as recommended for Affymetrix microarrays using the R/Bioconductor statistical programming environment (Tumor Analysis Best Practices Working Group, 2004). All subsequent data analysis procedures were also performed in R. Data were normalised using GCRMA and filtered using a variance of ⩾0.1. Differential gene expression analysis was performed using the normalised and filtered set of probes with ANOVA and eBayes. An adjusted P⩽0.05 after correction for multiple testing (Benjamini–Hochberg) was considered as significantly different between the two groups. Gene set enrichment analysis (GSEA) was performed as previously described using the C1, C2, C3 and C5 gene sets (Subramanian et al, 2005; Hamoudi et al, 2010). Leading edge analysis was performed to identify the most biologically important genes within sets that were identified as differentially regulated between AK and SCC. Gene ontology (GO) was performed using the online bioinformatic tool DAVID. To assess the effect of differentiation status and immune status on the expression profile, unsupervised hierarchical clustering of the GCRMA normalised and filtered probes (n=22 926) was performed using Ward average linkage. A Fisher’s exact test was used to measure the compactness of clustering for each variable.

Quantitative real-time PCR

Quantitative real-time PCR (qRT-PCR) was performed on an independent set of biopsies. Total RNA was extracted from fresh-frozen biopsies by homogenisation under liquid nitrogen. In total, RNA was extracted from five normal skin samples, nine AKs, nine WD SCC and 9 MD SCC using the RNeasy Mini kit with on-column DNAse digestion (Qiagen, Manchester, UK) according to the manufacturer’s instructions. First-strand cDNA was synthesised from 500 ng RNA using Superscript III First-Strand Synthesis SuperMix, with a combination of 2.5 μ M oligo dT and random primers (2.5 ng μl−1) (Invitrogen, Paisley, UK). cDNA was diluted 1 : 10 before real-time quantitative PCR (QPCR) analysis. Quantitative PCRs were run in triplicate on the ABI 7500 Real time PCR machine (Invitrogen) using 1 μl diluted cDNA and Brilliant II SYBR Green QPCR Mastermix (Stratagene, Wokingham, UK) with the two-step cycling protocol recommended by the manufacturer. The PCR primers were intron spanning with the exception of those for JUN, which has a single exon (primer sequences are available in Supplementary methods). Relative expression data were calculated by the equation 2ΔΔCt, using the endogeneous control gene RPS9 as the internal reference. Data were normalised to the mean Ct value from five normal (non-lesional) skin samples. For each sample, cDNA synthesis was also performed in the absence of reverse transcriptase to exclude the possibility of genomic contamination. Statistical analysis was performed using a two-tailed t-test, with a P<0.05 considered as statistically significant.

Results

Genes involved in epidermal differentiation, angiogenesis, taxis, proliferation and adhesion are differentially expressed between AK and cSCC

To identify DEGs between AK and cSCC, expression microarray data from 30 cSCC and 10 AKs were analysed by two different statistical algorithms: ANOVA and eBayes (Supplementary Table 1). The probes that were identified as differentially expressed by both methods were considered as a robust set of DEGs between AK and cSCC (adjusted P<0.05). In total, this included 239 probes (131 were upregulated in the cSCC and 108 were downregulated), which corresponded to 196 annotated genes (Tables 2A, 2b). This probe list notably included two upregulated genes that have been consistently identified as overexpressed in cSCC from previous microarray studies (MMP1 and CDH3), as well as downregulation of the differentiation-specific keratin, KRT10, which is also consistently reported in cSCC (Van Haren et al, 2009).

Table 2a Differentially expressed genes downregulated between AK and cSCC, as identified by both ANOVA and eBayes
Full size table
Table 2b Differentially expressed genes upregulated between AK and cSCC, as identified by both ANOVA and eBayes
Full size table

Clustering of the samples using these probes revealed two main groups, one containing the majority of cSCC (20 out of 30) and the second containing the AK and remaining 10 cSCC (Figure 1). The cluster containing only cSCC included 73% (11 out of 15) of the MD–PD tumours, compared with 40% (6 out of 15) of the WD tumours, although this difference did not reach statistical significance using a two-tailed Fisher’s exact test. There were no other significant clinical or pathological features relating to age at diagnosis, gender, immune status or tumour location (sun exposed vs non-sun exposed) that accounted for separation of the cSCC across the two groups. Careful laser capture microdissection of tissue samples used for RNA extraction makes it unlikely that stromal contamination was responsible for the clustering of some cSCC with the AK samples. Rather, it may indicate that there is a spectrum of tumour development. AK10 was clustered on an outlying clade within the same cluster as the other AK samples. There were no outstanding clinical or pathological features of this case (male, aged 80, immunocompetent, sun-exposed site) that accounted for its position within the clustering, again suggesting that it could represent natural variation among the samples.

Figure 1
figure 1

Heat map of 239 overlapping differentially expressed probes identified by ANOVA and eBayes analysis of cSCC vs AK. Characterisation bars beneath the dendrogram highlight key clinicopathological variables. Abbreviations: AK= actinic keratosis; IC= immunocompetent; IS= immunosuppressed; MD= moderately differentiated; PD= poorly differentiated; SCC= squamous cell carcinoma; WD= well differentiated.

Full size image

To identify biological processes represented by the DEGs, GO was performed using the online bioinformatic tool, DAVID. This revealed terms associated with angiogenesis, epidermal development and differentiation, taxis, proliferation and adhesion (Supplementary Table 2). In keeping with these ontology terms, some of the genes that showed the largest fold-change downregulation in the cSCC were those involved in terminal epidermal differentiation, including filaggrin (FLG), filaggrin family member 2 (FLG2), loricrin (LOR) and late cornified envelope 1B (LCE1B). Those that were strongly upregulated in cSCC included the matrix metallopeptidases MMP1 and MMP10, and a set of genes associated with migration and focal adhesion (including filamin A (FLNA), FLNB, integrin alpha 5 (ITGA5), jun proto-oncogene (JUN), met proto-oncogene, laminin-5 gamma2 chain (LAMC2), p21 protein (Cdc42/Rac)-activated kinase 2 (PAK2), vascular endothelial growth factor A (VEGFA) and Zyxin (ZYX).

qRT-PCR validates microarray differential gene expression analysis

To validate the findings of the microarray analysis, eight genes identified by both statistical methods (ACVR2A, ID4, INHBA, JUN, MET, MMP10, PAK2 and PTHLH) were assessed by qRT-PCR in an independent series of 27 lesions (9 AKs, 9 WD SCC and 9 MD SCC). In addition, a further three genes that were identified by the ANOVA only were validated (ADAM17, ANXA1 and PLK1), as well as MAPK1 that was identified by eBayes only. In all cases, the qRT-PCR data confirmed significant upregulation or downregulation between the AK and either the WD SCC, MD SCC or combined values (WD+MD), suggesting good reliability for the genes identified by this microarray study (Figure 2).

Figure 2
figure 2

qRT-PCR analysis of 12 differentially expressed genes between AK and cSCC. Solid black squares indicate outlier samples. For the purposes of comparison, the cSCC has been split into WD and MD categories, as well as overall values (SCC). Abbreviations: AK= actinic keratosis; MD= moderately differentiated; SCC= squamous cell carcinoma; WD= well differentiated.

Full size image

GSEA identifies overrepresentation of members of the MAPK pathway genes in cSCC

GSEA was performed to investigate the expression of genes within a priori established gene sets, to identify whether a group of genes that share a defined category (such as biological function, pathway or chromosomal location) show collective differences in expression between AK and cSCC. This is a useful tool for extrapolating biologically relevant information from microarray studies, and is particularly robust as it focuses on a group of genes that share a category, rather than individual genes that may be subject to signal-to-noise background. In total, 28 significantly enriched gene sets were identified (P⩽0.05, FDR <0.25). These could be broadly classified into four categories: signal transduction or apoptosis, cell cycle or replication, metabolic processes or cell motility and miscellaneous (which included gene sets related to the expression profile of other cancer types) (Table 3; Figure 3A and B). Analysis of the leading edge genes underlying the enrichment of each individual gene set revealed that many were consistently represented, suggesting that they strongly influenced the expression pattern in SCC. These included multiple members of the MAPK family (including JUN, FOS, MET, MAPK1, MAP2K2 and MAP3K5) and the apoptotic response (including PAK2, BAX, DEDD, DAXX and DFFA), suggesting that the interplay of the MAPK pathway and apoptotic response may be a critical determinant of AK to cSCC progression. No gene sets were identified from chromosomal and cytogenetic locations or conserved cis-regulatory motifs.

Table 3 GSEA analysis
Full size table
Figure 3
figure 3

Enrichment of differences in the MAPK signalling pathway (HSA04010) between AK and cSCC as identified by GSEA, and cluster analysis of cSCC. (A) Enrichment of genes in the MAPK signalling pathway. The blue line indicates the enrichment score (ES) and the black vertical lines beneath correspond to individual genes within the set, ranked according to their enrichment. If no enrichment was present, then the genes would be distributed equally from left to right. (B) Heat map of the enriched genes, showing their level of expression across AK and cSCC. Red indicates high expression and green indicates low expression. Each row corresponds to a gene within the MAPK gene set, while each column corresponds to an individual sample. (C) Cluster dendrogram revealing two main clusters of cSCC, separated predominantly by differentiation status (WD vs M/PD). Abbreviations: AK= actinic keratosis; IC= immunocompetent; IS= immunosuppressed; MD= moderately differentiated; PD= poorly differentiated; SCC= squamous cell carcinoma; WD= well differentiated.

Full size image

WD tumours display a distinct expression profile to MD and PD cSCC

To investigate the effect of different variables on the expression profile of the SCC, unsupervised hierarchical clustering of the cSCC only was performed. This revealed two main clusters: cluster 1 contained predominantly WD tumours and cluster 2 contained mostly MD and PD tumours (Figure 3C). The clusters show that WD tumours form the most significant cluster (P=0.002) with 13 out of 15 (86.7%) WD tumours in cluster 1, clearly separating them from the majority of MD and PD tumours (Figure 3D). Moderately differentiated and PD tumours themselves form weaker clusters (P=0.05 and P=0.06, respectively), however when combined they show a highly significant clustering pattern (P=0.002). This is in keeping with the histological heterogeneity of higher grade SCC, which often shows areas of varying differentiation rather than a consistent differentiation status across the whole tumour. Immune status showed no influence on clustering across all tumours, with around half of each status present in both clusters (P=0.55 for both).

Differential gene expression analysis (eBayes and ANOVA) identified no DEGs between WD cSCC and their MD/PD counterparts after correction for multiple testing, likely due to the small number of samples in each subgroup. Similarly, no genes were identified as differentially expressed between tumours from immunosuppressed vs immunocompetent patients (examined within the category of WD cSCC only to avoid bias from the differentiation status of the tumours).

Discussion

Cutaneous SCC is the most highly mutated of all cancers, and consequently displays a complex genetic background (Durinck et al, 2011). Our study identified a set of 196 DEGs between AK and cSCC, with enrichment for genes involved in loss of differentiation (including downregulation of FLG, FLG2, LOR, LCE1B, KRT9 and KRT10 in cSCC), and gain of invasive properties such as extracellular matrix remodelling and cell migration. We have identified the MAPK pathway as pivotal to many of these processes, and specifically found the oncogenes JUN and MET to be overexpressed in cSCC. This is particularly significant as the contribution of oncogenes to cSCC development has remained elusive, despite the identification of a number of key tumour suppressor genes such as TP53, CDKN2A and, recently, NOTCH (Jonason et al, 1996; Brown et al, 2004; Durinck et al, 2011; Wang et al, 2011). Furthermore, we have shown that WD cSCC display a distinct expression profile to MD or PD tumours, which has implications for future targeted therapies.

Actinic keratoses represent a precancerous stage in cSCC genesis, and provide an opportunity to characterise processes relevant to the progression of skin carcinogenesis. Many previous expression microarray studies of the transition from AK to cSCC have identified few (<10), or no, DEGs between these two disease states (Nindl et al, 2006; Padilla et al, 2010; Ra et al, 2011). In a systematic review of six microarrays studies, it was suggested that the low consensus for differential gene expression in cSCC resulted from clinical and methodological differences, including small sample sizes, tissue processing, different baseline comparators (normal skin or AK) and varied bioinformatic approaches (Van Haren et al, 2009). Importantly, the findings from our study are consistent with recently published work that combined laser capture microdissection with reverse phase protein microarray analysis (Einspahr et al, 2012). This study also identified activation of the MAPK pathway in cSCC compared with AK and normal skin, highlighting the power of stringent laser capture microdissection for improved consistency across studies. In addition, a further study that focused on a smaller series of AK and cSCC from OTRs also identified the MAPK pathway to be activated in cSCC (Hameetman et al, 2013). Although we cannot be sure the individual AK used across these studies had the potential to progress, the pattern of differential expression supports consistent differences in the biology of these tumours with a shift in invasive properties between AK and cSCC.

MAPK pathway in cSCC

The oncogenes MET, JUN and PAK2 were all overexpressed in cSCC compared with AK. Both MET and JUN are direct components of the MAPK pathway, which was identified by GSEA to be a key difference between AK and cSCC. MAPK signalling has a central role in regulating growth and survival of cancer cells, and the RAF/MEK/ERK pathway is dysregulated in approximately one-third of human cancers (Dhillon et al, 2007). Consequently, a number of small molecule inhibitors have been designed to target specific steps in MAPK signalling (including MET) and are now in clinical trials (Tu et al, 2010). The potential importance of this pathway in cSCC has recently been highlighted by the rapid development of cSCC in a proportion of patients treated with the multi-kinase inhibitor, Sorafenib, or the selective BRAF V600E inhibitors, Vemurafenib and Dabrafenib (Ribas and Flaherty, 2011). There is evidence that this is caused by paradoxical activation of the MAPK pathway, which is proposed to cooperate with pre-existing somatic UV-induced mutations in key oncogenes and tumour suppressors such as H-RAS and TP53 (Hatzivassiliou et al, 2010; Heidorn et al, 2010; Poulikakos et al, 2010; Arnault et al, 2012).

We also observed strong upregulation of the collagenase MMP1 and the stromeolysin MMP10 in cSCC, both of which are associated with invasion and metastasis via degradation of ECM proteins (Overall and Kleifeld, 2006). MMP1 upregulation is one of the few consistent findings across multiple cSCC microarray studies and, similar to targets in the MAPK pathway, has therapeutic potential for inhibition (Van Haren et al, 2009; Gialeli et al, 2011). In addition, pro-metastatic molecules such as PTHLH and osteopontin (SPP1) were upregulated in our data set, whereas putative tumour suppressors such as RORA, PDZD2 and AZGP1 were found to be downregulated (Tam et al, 2008; Yip et al, 2011; Wang et al, 2012).

Transcriptome profiles reflect histological subtypes of cutaneous SCC

Histological subtypes of cSCC are associated with distinct clinical characteristics and behaviour, with poor differentiation an independent prognostic factor for metastatic potential (Mourouzis et al, 2009; Breuninger et al, 2012; Jambusaria-Pahlajani et al, 2013). We have previously shown that WD cSCC display a distinct genetic profile to MD and PD tumours, with fewer genomic aberrations overall (Purdie et al, 2009). In addition, single gene-specific studies have reported differential expression or genomic aberration patterns between cSCC of different histological grades, including higher expression of LRIG-1 in WD cSCC, an association of PTPRD deletions with PD cSCC and overexpression of Survivin in PD tumours (Lo Muzio et al, 2001; Tanemura et al, 2005; Lambert et al, 2012). Our transcriptome data support the genomic profiling by finding that WD tumours display a distinct expression profile to MD and PD cSCC. While we were unable to identify specific genes that are differentially expressed between the subgroups (most likely due to small sample numbers), our data imply that the histological subtypes of cSCC should be defined when addressing the molecular pathogenesis, prevention and treatment of these tumours.

In conclusion, our study is the first to describe widespread transcriptome changes between precancerous AK and cSCC in both immunocompetent and immunosuppressed individuals. These data highlight several known oncogenes, and reveal a broad spectrum of disrupted cellular processes that are altered in cSCC. The MAPK pathway shows pivotal changes, and offers a new approach to targeted preventative and therapeutic strategies. Finally, our data indicate that histological subtypes of cSCC represent distinct entities at the transcriptome level, and may in the future require individually tailored therapeutic approaches.