Key differences identified between actinic keratosis and cutaneous squamous cell carcinoma by transcriptome profiling.

BACKGROUND: Cutaneous squamous cell carcinoma (cSCC) is one of the most common malignancies in fair-skinned populations worldwide and its incidence is increasing. Despite previous observations of multiple genetic abnormalities in cSCC, the oncogenic process remains elusive. The purpose of this study was to elucidate key molecular events associated with progression from premalignant actinic keratoses (AKs) to invasive cSCC by transcriptome profiling.
METHODS: We combined laser capture microdissection with the Affymetrix HGU133 Plus 2.0 microarrays to profile 30 cSCC and 10 AKs.
RESULTS: We identified a core set of 196 genes that are differentially expressed between AK and cSCC, and are enriched for processes including epidermal differentiation, cell migration, cell-cycle regulation and metabolism. Gene set enrichment analysis highlighted a key role for the mitogen activated protein kinase (MAPK) pathway in cSCC compared with AK. Furthermore, the histological subtype of the tumour was shown to influence the expression profile.
CONCLUSION: These data indicate that the MAPK pathway may be pivotal to the transition from AK to cSCC, thus representing a potential target for cSCC prevention. In addition, transcriptome differences identified between cSCC subtypes have important implications for future development of targeted therapies for this malignancy.

Skin cancers are by far the most common malignancies in fair-skinned populations, with an incidence now reaching epidemic proportions (Madan ; Rogers ). 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 ; 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 ; Kraus ). 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 ; Harwood , 2013).

An estimated 65% of cSCC arise from precursor lesions termed actinic keratoses (AKs) (Criscione ). 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 ; Memon ). 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 ). 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 ).

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 ). 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 ; Wang ). 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 ; Popp ; Ashton ; Clausen ; Purdie , 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 ). 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 ). 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 ; Kathpalia ; Nindl ).

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

Sample name	AK/SCC	Patient no.	Differentiation status	Immune status	Gender	Location	Age (years)
AK1	AK	8	NA	IC	M	Leg	52
AK2	AK	10	NA	IS	M	Scalp	65
AK3	AK	9	NA	IS	M	Hand	69
AK4	AK	11	NA	IC	M	R hand	89
AK5	AK	2	NA	IS	F	Hand	71
AK6	AK	1	NA	IS	M	Hand	63
AK7	AK	4	NA	IS	M	Hand	57
AK8	AK	5	NA	IS	M	Hand	57
AK9	AK	6	NA	IS	M	Scalp	72
AK10	AK	7	NA	IC	M	Hand	80
SCC1	SCC	1	WD	IS	M	Elbow	63
SCC2	SCC	2	WD	IS	F	Hand	69
SCC3	SCC	3	WD	IC	M	Calf	74
SCC4	SCC	4	WD	IS	M	Upper arm	56
SCC5	SCC	5	WD	IS	M	Hand	49
SCC6	SCC	9	WD	IS	M	Scalp	70
SCC7	SCC	10	WD	IS	M	Temple	64
SCC8	SCC	11	WD	IC	M	Temple	88
SCC9	SCC	14	WD	IC	M	Temple	77
SCC10	SCC	15	WD	IC	M	Pinna	NA
SCC11	SCC	16	WD	IS	M	Neck	30
SCC12	SCC	22	WD	IS	M	Scalp	56
SCC13	SCC	23	WD	IS	F	Finger	47
SCC14	SCC	30	WD	IC	M	Groin	67
SCC15	SCC	31	WD	IS	M	Forearm	46
SCC16	SCC	6	MD	IS	M	Pinna	73
SCC17	SCC	8	MD	IC	M	Chest	52
SCC18	SCC	17	MD	IC	M	Temple	77
SCC19	SCC	21	MD	IS	F	Shin	66
SCC20	SCC	24	W-MD	IS	M	Pinna	80
SCC21	SCC	25	MD	IS	M	Pinna	58
SCC22	SCC	26	MD	IS	M	Neck	41
SCC23	SCC	27	MD	IS	M	Chest	51
SCC24	SCC	28	MD	IC	M	Temple	69
SCC25	SCC	12	M-PD	IC	F	Calf	86
SCC26	SCC	13	M-PD	IS	F	Finger	59
SCC27	SCC	18	M-PD	IS	F	Foot	73
SCC28	SCC	19	PD	IS	M	Cheek	75
SCC29	SCC	20	PD	IS	M	Pinna	68
SCC30	SCC	29	M-PD	IS	M	Neck	67

Abbreviations: AK=actinic keratosis; F=female; IC=immunocompetent; IS=immunosuppressed; M=male; MD=moderately differentiated; NA=not available; PD=poorly differentiated; SCC=squamous cell carcinoma; WD=well differentiated.

Expression microarray analysis

Laser capture microdissection and RNA extraction were performed as previously described (Lambert ), 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 ; Hamoudi ). 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 ).

Table 2A

Differentially expressed genes downregulated between AK and cSCC, as identified by both ANOVA and eBayes

Gene name	Log2 FC	Adj P-value	Gene name	Log2 FC	Adj P-value
FLG2	−6.0	2.96E−02	LOC283070	−2.2	3.36E−03
KRT9	−5.8	8.31E−04	SEPP1	−2.1	1.57E−02
KRT77	−5.4	8.31E−04	SERTAD4	−2.0	1.34E−03
FLG	−5.3	3.30E−02	GLDN	−1.9	2.87E−02
LOR	−5.2	1.73E−02	UBL3	−1.9	3.34E−03
LCE1B	−5.2	3.67E−02	LYPLAL1	−1.9	2.41E−02
SERPINB12	−5.2	1.55E−02	ARL5A	−1.8	8.24E−03
AZGP1	−4.9	6.23E−04	HNRPLL	−1.8	3.14E−02
CD36	−4.7	5.44E−04	BOC	−1.8	7.15E−03
AADACL2	−4.1	1.43E−02	GGTA1	−1.8	2.91E−02
HPGD	−3.8	3.22E−03	ELOVL6	−1.7	4.14E−03
SERPINA12	−3.6	4.57E−02	SPATA6	−1.7	1.07E−02
TNFRSF19	−3.3	4.82E−02	C1orf96	−1.7	2.05E−03
MUC15	−3.3	6.07E−03	IDE	−1.7	1.44E−02
MFAP3L	−3.3	1.63E−02	CREBL2	−1.7	1.55E−03
ATP6V1C2	−3.2	4.90E−03	PGRMC2	−1.6	3.43E−02
ELOVL4	−3.1	3.34E−03	KRT10	−1.6	1.51E−02
BPIL2	−3.1	3.46E−03	SLC30A1	−1.6	7.91E−04
TGFBR3	−3.0	5.44E−04	ZDHHC23	−1.5	1.55E−02
EDNRB	−2.9	1.82E−05	TMEM45A	−1.5	3.36E−03
PDZD2	−2.9	1.55E−03	ATP7A	−1.5	1.45E−02
ITM2A	−2.9	1.67E−02	C5orf41	−1.5	3.86E−02
MARCH3	−2.8	7.15E−03	KIAA1370	−1.5	7.15E−03
LAMB4	−2.8	4.69E−02	OLFM2	−1.4	3.86E−02
EGR3	−2.7	7.76E−03	PARM1	−1.4	2.25E−02
ZNF682	−2.7	1.53E−02	TCP11L2	−1.4	8.24E−03
METTL7A	−2.6	1.21E−03	SECISBP2L	−1.3	1.64E−02
RORA	−2.6	5.18E−04	MIPEP	−1.3	3.66E−02
QPCT	−2.5	1.19E−02	OXCT1	−1.3	2.87E−02
LGR6	−2.5	2.19E−02	ATAD2B	−1.3	3.72E−02
ID4	−2.5	3.28E−02	CNOT6L	−1.2	2.32E−02
HLF	−2.5	4.57E−03	SLC26A11	−1.2	1.45E−02
GRAMD1C	−2.4	2.35E−03	RNLS	−1.2	3.01E−02
LONRF1	−2.3	8.31E−04	ACVR2A	−1.2	2.05E−03
PPM1L	−2.3	2.20E−02	CEP68	−1.2	2.05E−02
GLRX	−2.2	2.43E−03	DYNC1LI1	−1.1	6.22E−03
ITGBL1	−2.2	2.18E−02	TANC1	−1.1	1.34E−02
FAM13C	−2.2	1.07E−02	ATG2B	−1.1	1.21E−03
MATN2	−2.2	3.13E−02	WDFY3	−1.1	2.97E−02
CADM1	−2.2	1.36E−02	KIAA1012	−1.0	1.57E−02
GAS7	−2.2	9.34E−04	NARG1L	−1.0	2.57E−02
TRAM1L1	−2.2	4.33E−02

Abbreviations: AK=actinic keratosis; ANOVA=analysis of variance; cSCC= cutaneous squamous cell carcinoma; FC=fold change.

Table 2B

Differentially expressed genes upregulated between AK and cSCC, as identified by both ANOVA and eBayes

Gene name	Log2 FC	Adj P-value	Gene name	Log2 FC	Adj P-value
MMP1	4.6	3.46E−02	CDC42BPB	1.4	6.22E−03
MMP10	4.4	3.77E−04	S100A6	1.4	4.56E−02
CXCL1	4.1	7.91E−04	IL4R	1.4	1.33E−02
INHBA	3.7	6.07E−03	SPATS2	1.4	3.72E−02
SPP1	3.7	7.91E−04	ZCCHC10	1.4	1.33E−02
PTHLH	3.5	7.81E−04	SRM	1.4	4.39E−02
SH2D5	3.5	2.91E−03	JAG1	1.4	4.14E−02
ALDH1A3	3.3	1.85E−02	JOSD1	1.4	1.43E−02
LAMC2	3.1	2.35E−03	ZYX	1.3	2.76E−02
GPRC5A	2.9	3.88E−02	NOLC1	1.3	2.87E−02
PLAUR	2.8	3.77E−04	GLT25D1	1.3	3.03E−02
GPR68	2.7	2.17E−02	CTSB	1.3	2.97E−02
STC1	2.3	4.33E−02	YKT6	1.3	2.90E−02
ERO1L	2.2	3.46E−03	HIST1H2BI	1.3	2.37E−02
ANXA2	2.2	1.97E−02	HN1L	1.3	4.53E−02
MTHFD1L	2.2	1.43E−02	FLNB	1.3	2.95E−02
SLC16A3	2.2	1.33E−03	MAP3K13	1.3	3.75E−02
EIF4EBP1	2.2	1.64E−02	PLEC1	1.3	1.73E−02
PLAU	2.2	6.31E−03	RPS6KA4	1.3	2.05E−03
FLNA	2.1	1.08E−03	IFRD2	1.3	1.19E−02
IGF2BP2	2.1	2.35E−03	ANKLE2	1.2	2.97E−02
ITGA5	2.1	8.31E−04	TSC2	1.2	3.32E−02
HRH1	2.1	4.20E−02	SLC9A1	1.2	5.44E−04
KLF7	2.1	8.19E−03	TES	1.2	2.16E−02
FAM83A	2.0	2.16E−02	LRRC8A	1.2	3.01E−02
TRIP10	2.0	7.15E−03	N4BP1	1.2	6.07E−03
NRP2	2.0	1.85E−02	DNAJC8	1.2	7.11E−03
PDIA5	2.0	2.43E−03	PLOD1	1.2	1.96E−03
SERPINB1	2.0	1.45E−02	S100A13	1.2	2.88E−02
SRGAP1	1.9	1.43E−02	PDXK	1.2	3.30E−02
TYMP	1.9	2.29E−02	TCF3	1.1	3.34E−03
SESTD1	1.8	2.88E−02	PVR	1.1	3.63E−03
RASIP1	1.8	1.57E−02	VPS72	1.1	1.08E−02
C11orf17	1.8	1.72E−02	PAK2	1.1	1.22E−02
NOP56	1.8	1.40E−02	CAD	1.1	2.26E−02
C16orf57	1.8	1.05E−02	RUVBL1	1.1	4.06E−02
TGFA	1.7	1.19E−02	TICAM1	1.1	3.63E−03
MET	1.7	5.00E−03	B4GALT1	1.1	1.57E−02
CARHSP1	1.7	1.05E−02	GPR153	1.1	8.57E−04
ZAK	1.7	5.44E−04	HNRNPU	1.1	2.97E−02
HN1	1.7	2.11E−02	ZNF697	1.1	7.15E−03
FTL	1.6	1.84E−02	CDC27	1.0	1.55E−03
TNFRSF12A	1.6	7.15E−03	MARK3	1.0	3.52E−02
CDCP1	1.6	7.66E−03	CEBPD	1.0	4.02E−02
RFTN1	1.6	3.30E−02	AP2A2	1.0	3.34E−03
KLF6	1.6	4.74E−02	ZBED4	1.0	3.98E−02
TRIO	1.6	4.06E−02	CTSLL3	1.0	1.69E−02
EHD2	1.5	4.73E−02	BRD9	1.0	3.86E−02
TCOF1	1.5	4.57E−02	SF3A2	0.9	1.33E−03
NRIP3	1.5	3.34E−03	HEATR2	0.9	3.75E−02
CDH3	1.5	2.60E−02	RFT1	0.9	2.06E−02
JUN	1.5	1.50E−02	SOCS1	0.8	1.28E−02
CDC20	1.5	6.14E−03	GTSE1	0.8	3.77E−04
RRBP1	1.5	2.17E−02	GALE	0.7	2.37E−02
BOP1	1.4	1.73E−02	ATP13A1	0.7	4.89E−02
SDC4	1.4	4.89E−02	VPS37B	0.6	3.72E−02
VEGFA	1.4	4.12E−02

Abbreviations: AK=actinic keratosis; ANOVA=analysis of variance; cSCC= cutaneous squamous cell carcinoma; FC=fold change.

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

Heat map of 239 overlapping differentially expressed probes identified by ANOVA and eBayes analysis of cSCC 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.

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

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.

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

GSEA category	Size	Source	ES	NES	NOM P-value	FDR q-value	Tag %	Gene %
Signal transduction or apoptosis
MAPK signalling pathway	121	HSA04010	−0.36389	−1.5201	0.0245	0.2425	0.322	0.176
TNFR1 pathway	20	na	−0.48995	−1.7601	0.0080	0.1775	0.25	0.0566
FAS pathway	21	na	−0.47668	−1.6076	0.0351	0.2479	0.381	0.156
ST FAS signalling pathway	36	na	−0.61198	−2.0475	0.0000	0.0167	0.361	0.156
Death pathway	20	na	−0.57202	−1.7337	0.0120	0.1189	0.4	0.156
MET pathway	25	na	−0.46801	−1.5399	0.0440	0.2416	0.4	0.181
Cell cycle or replication
M phase	65	GO:0000279	−0.62169	−1.8425	0.0197	0.1063	0.631	0.246
M phase of mitotic cell cycle	54	GO:0000087	−0.65812	−1.8475	0.0100	0.1466	0.704	0.258
Cell-cycle process	107	GO:0022402	−0.57127	−1.7439	0.0301	0.1498	0.598	0.246
Cell-cycle phase	96	GO:0022403	−0.56841	−1.7673	0.0275	0.1550	0.583	0.246
Mitotic cell cycle	92	GO:0000278	−0.59663	−1.7838	0.0235	0.1875	0.609	0.246
Mitosis	51	GO:0007067	−0.66288	−1.865	0.0080	0.2242	0.725	0.258
Condensed chromosome	20	na	−0.61929	−1.7196	0.01536	0.17401	0.6	0.226
Metabolic processes
Negative regulation of cellular protein metabolic process	20	GO:0032269	−0.60511	−1.7566	0.0037	0.1511	0.5	0.209
Negative regulation of protein metabolic process	21	GO:0051248	−0.5716	−1.7147	0.0019	0.1646	0.476	0.209
Ribonucleoprotein complex	48	na	−0.53241	−1.7659	0.01416	0.17705	0.583	0.275
Pyrophosphatase activity	99	GO:0016462	−0.39159	−1.6548	0.0000	0.1865	0.354	0.18
Hydrolase activity acting on acid anhydrides	101	GO:0016817	−0.38398	−1.6483	0.0000	0.1689	0.347	0.18
Nucleoside triphosphatase activity	93	GO:0017111	−0.40821	−1.7111	0.0019	0.2021	0.366	0.18
Protein kinase binding	30	GO:0019901	−0.52176	−1.643	0.0118	0.1559	0.3	0.113
Pyrimidine metabolism	33	na	−0.55439	−1.6614	0.0196	0.1865	0.515	0.248
Cell motility
Leading edge	29	GO:0031252	−0.54978	−1.5753	0.0301	0.1725	0.483	0.198
Ruffle	20	GO:0001726	−0.5811	−1.5859	0.0245	0.1805	0.5	0.185
Cell projection	51	GO:0042995	−0.46711	−1.5155	0.0447	0.2203	0.431	0.213
Miscellaneous
VHL RCC up	66		−0.51542	−1.8918	0.0000	0.2099	0.5	0.281
Renal cell carcinoma	42	HSA05211	−0.60243	−1.836	0.0000	0.1507	0.548	0.213
Bladder cancer	30	HSA05219	−0.53204	−1.7384	0.0018	0.2485	0.333	0.138
Acute myeloid leukaemia	31	HSA05221	−0.4781	−1.5209	0.0371	0.2498	0.323	0.138

Abbreviations: ES=enrichment score; FDR=false discovery rate; GO=gene ontology; GSEA=gene set enrichment analysis; MAPK=mitogen activated protein kinase; NES=nominal enrichment score; VHL RCC=Von Hippel-Lindau renal cell carcinoma..

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.

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 ). 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 ; Brown ; Durinck ; Wang ). 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 ; Padilla ; Ra ). 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 ). Importantly, the findings from our study are consistent with recently published work that combined laser capture microdissection with reverse phase protein microarray analysis (Einspahr ). 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 ). 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 ). 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 ). 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 ; Heidorn ; Poulikakos ; Arnault ).

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 ; Gialeli ). 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 ; Yip ; Wang ).

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 ; Breuninger ; Jambusaria-Pahlajani ). We have previously shown that WD cSCC display a distinct genetic profile to MD and PD tumours, with fewer genomic aberrations overall (Purdie ). 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 ; Tanemura ; Lambert ). 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.