Molecular profiling of cutaneous squamous cell carcinomas and actinic keratoses from organ transplant recipients

Background The risk of developing cutaneous squamous cell carcinoma (SCC) is markedly increased in organ transplant recipients (OTRs) compared to the normal population. Next to sun exposure, the immunosuppressive regimen is an important risk factor for the development of SCC in OTRs. Various gene mutations (e.g. TP53) and genetic alterations (e.g. loss of CDKN2A, amplification of RAS) have been found in SCCs. The aim of this genome-wide study was to identify pathways and genomic alterations that are consistently involved in the formation of SCCs and their precursor lesions, actinic keratoses (AKs). Methods To perform the analysis in an isogenic background, RNA and DNA were isolated from SCC, AK and normal (unexposed) epidermis (NS) from each of 13 OTRs. Samples were subjected to genome-wide expression analysis and genome SNP analysis using Illumina’s HumanWG-6 BeadChips and Infinium II HumanHap550 Genotyping BeadChips, respectively. mRNA expression results were verified by quantitative PCR. Results Hierarchical cluster analysis of mRNA expression profiles showed SCC, AK and NS samples to separate into three distinct groups. Several thousand genes were differentially expressed between epidermis, AK and SCC; most upregulated in SCCs were hyperproliferation related genes and stress markers, such as keratin 6 (KRT6), KRT16 and KRT17. Matching to oncogenic pathways revealed activation of downstream targets of RAS and cMYC in SCCs and of NFκB and TNF already in AKs. In contrast to what has been reported previously, genome-wide SNP analysis showed very few copy number variations in AKs and SCCs, and these variations had no apparent relationship with observed changes in mRNA expression profiles. Conclusion Vast differences in gene expression profiles exist between SCC, AK and NS from immunosuppressed OTRs. Moreover, several pathways activated in SCCs were already activated in AKs, confirming the assumption that AKs are the precursor lesions of SCCs. Since the drastic changes in gene expression appeared unlinked to specific genomic gains or losses, the causal events driving SCC development require further investigation. Other molecular mechanisms, such as DNA methylation or miRNA alterations, may affect gene expression in SCCs of OTRs. Further study is required to identify the mechanisms of early activation of NFκB and TNF, and to establish whether these pathways offer a feasible target for preventive intervention among OTRs.


Background
Non-melanoma skin carcinomas, comprising cutaneous squamous cell carcinoma (SCC) and basal cell carcinoma (BCC), are the most common malignancies in fair-skinned Caucasians. The risk of developing cutaneous SCC is markedly increased (100-fold) in organ transplant recipients (OTRs) compared to the normal population [1,2] and is associated with the immunosuppressive therapy needed to prevent graft rejection [3]. Moreover, some studies have suggested that SCCs from OTRs display a more aggressive behavior [4].
Like other cancers, the development of cutaneous SCC is thought to be a multi-step process, involving sequential acquisition of genetic changes. Sun (UVR) exposure is the principal carcinogen inducing both immunologic tumortolerance and DNA damage that can lead to mutations in oncogenes and tumor suppressor genes (TSGs) [5,6]. Other risk factors include impaired immune surveillance, as seen in OTRs, and human papillomavirus (HPV) infections [7]. However, the oncogenic potential of HPV in cutaneous SCCs remains unclear [8]. In a recent transcriptome sequencing study no evidence of active viral gene expression was found in large group of SCCs [9].
Cutaneous SCCs are thought to develop from precursor lesions, actinic keratoses (AKs). Even though most AKs (80-99 %) do not progress to SCC, there is still both histological and molecular evidence to support this hypothesis (reviewed by C.J. Ko [10]). UVR-related mutations inactivating the TP53 tumor suppressor gene are very common in SCCs and AKs, both in immunocompetent patients and OTRs [5,11]. TP53 mutations are thought to occur very early in the development of skin cancer, since they are already present in microscopic clusters of keratinocytes in sun-exposed human skin [12]. These clusters are thought to precede tumor formation [13] and are more frequent in OTR compared to immunocompetent patients [14]. However, not every SCC contains a mutated TP53 gene [5].
Other common genetic changes found in cutaneous SCC development are the inactivation of the CDKN2A tumor suppressor gene (via mutation, promoter methylation and/or chromosomal loss) and mutation and/or amplification of RAS genes (reviewed in Boukamp et al. [15]). Both oncogenic changes are associated with the progression of AK to SCC. Also, (over-) expression of RAS downstream proteins (MAPKs and cyclins) has been reported in a subset of SCCs [16].
At the chromosomal level, several studies have shown that cutaneous SCCs can display complex karyotypes with large numbers of allelic imbalances [17][18][19]. Moreover, widespread gains and losses of chromosomal fragments have been reported to be already present in AKs [17,20].
The aim of our study was to identify genes and/or pathways that are consistently involved in the development of cutaneous SCC and their precursor lesions, AKs in OTRs. Thus far, the genome-wide transcriptional profiling studies on SCCs are diverse in nature, show little overlap in differentially expressed genes (DEGs), have limited sample size, use cell lines and do not integrate genomic and expression data [21][22][23][24]. This genome-wide study assesses in parallel both RNA and DNA alterations in a sizeable well-defined group of OTRs presenting with both AKs and SCCs. A better understanding of the pathological mechanisms leading to skin carcinomas in these patients could lead to novel molecular targets for effective intervention.

Patient material
Patients were selected from the large group of OTRs that visit the dermatology clinic of the LUMC at a regular basis to provide a well-documented sample series for our study. Moreover from these OTRs it was possible to select both SCCs and AKs as well as normal skin (NS) to circumvent inter-patient variability. Fifteen OTRs (7 females and 8 males) participated in this study ( Table 1). The average age of the patients when the selected SCC was diagnosed was 59 years (range 47-74 years). Sidematched normal skin in OTRs often already contains a fair amount of UV damage, a phenomenon described as 'field cancerisation' [25,26]. Therefore, (UV-unexposed) skin from the lower back was chosen as normal skin and served as internal negative control. Of patient 42, only SCC and AK samples were available, and for patient 15 only a SCC sample could be obtained. Samples of these two patients were only included in the Genome-wide expression analysis (GWEA). All patients received immunosuppressive drugs (Table 1).

Genome-wide SNP analysis of SCCs and AKs
Thirteen SCC samples and eleven AKs were subjected to genome-wide SNP array analysis together with the corresponding normal DNA from peripheral blood. Results are summarized in Table 2 and Figure 1. Six of thirteen SCCs showed chromosomal aberrations, ranging from a single loss of 9p to multiple losses and gains (amplification/duplication) (Table 2, Figure 1A). Copy number neutral loss of heterozygosity (LOH) was observed in 3 SCCs. Overall, more LOH than gain was seen (14 vs. 9). Also large deletions of (almost) whole chromosome arms were more frequent than small deletions (<40 Mb; 12 vs. 8). In several tumors the detected aberrations were not consistently present in >70% of the tumor sample (i.e. not a full signal of loss or gain), suggesting tumor heterogeneity (Table 2). Four out of thirteen (30%) SCCs demonstrated loss of 9p21 where CDKN2A is located.
Analysis of AK samples showed chromosomal aberrations in 6 of 11 samples, ranging from a single up to multiple (11) aberrations in an affected sample (Table 2, Figure 1B). Again, more losses and copy number neutral LOH than gains were seen (21 vs. 3). Generally, the aberrations observed in AKs, were less extensive than those in SCCs. No correlation was detected among the AKs regarding the extent of chromosomal damage and the relative level of dysplasia. However, similar to SCCs, several AKs demonstrated genetic heterogeneity.
Notably, there was no clear correlation between the number and type of aberrations found in the AK and SCC from the same OTR. However the number of samples in each group was not large enough to perform any statistical analysis.

Cluster and single-gene analyses of expression profiles of SCCs and AKs
Fifteen SCC, fourteen AK and thirteen NS samples were included in the genome wide expression analysis (GWEA). Filtering for annotated and expressed genes resulted in 15,969 probes for further analyses. Subsequent hierarchical clustering, based on expression levels, resulted in two very distinct groups (Figure 2A), separating normal and tumor samples. The tumor samples were divided in three sub-clusters in one arm of the dendogram. Principal component analysis also demonstrated good separation of normal and tumor samples, already in the first component ( Figure 2B). In the second component, AKs could be separated from SCCs. No correlation was seen with the type of organ transplanted. However, since 13 of 15 patients were kidney recipients, no definite conclusions can be made about this matter. The same could be concluded for the immunosuppressive regimen (data not shown).
Statistical analysis of the expression data was performed using limma. The following comparisons between sample groups were made: SCC vs. NS, SCC vs. AK and AK vs. NS. For each comparison a FDR <1% was used, and the comparisons were made with several log 2 transformed fold change (log 2 FC) restrictions, namely log 2 FC > 0.5, log 2 FC > 1.0 and no log 2 FCrestriction (Table 3). With log 2 FC > 0.5, 2,087 probes, representing 2,009 genes, were differentially expressed between SCC and NS (Additional file 1). Among the upregulated genes in SCC were several genes known to be involved in keratinocyte differentiation and skin cancer, including keratins (e.g. KRT16, KRT17 and KRT6), matrix metalloproteinases (MMPs), S100 molecules and small proline-rich proteins (SPRRs). Most of these genes were highly overexpressed with a log 2 FC > 2.0. The cutaneous T-cell attracting chemokine, CCL27, was the most downregulated gene in SCCs compared to normal skin (log 2 FC = − 4.2).
With log 2 FC > 0.5, 977 probes (947 genes) were differentially expressed between AK and NS samples and expression of 571 probes (547 genes) was significantly different between SCCs and AKs (Additional file 1).
In total for the log 2 FC > 0.5 analysis, 180 probes were present in all three comparisons as visualized by a VennDiagram (Figure 3, Additional file 1). These probes represented 173 genes that were consistently up-or downregulated during the development of SCC from NS, via AK, and can be considered as genes involved in gradual skin cancer progression. Since SCC, AK and NS samples were available from each patient, we investigated whether the differentially expressed genes between sample groups differed per patient. To this end, a patient factor was included in the design of limma; however this did not change the results (data not shown).

Geneset enrichment analyses of expression of SCCs and AKs
In addition to single-gene analysis, affected molecular pathways were investigated employing gene set enrichment analysis (GSEA) using two different approaches. Firstly, differentially expressed probes between the sample groups from the limma analysis with log 2 FC > 0.5 (Additional file 1) were analyzed for gene set enrichment with DAVID Bioinformatic Resources. All annotated and expressed probes (n = 15,969) were used as background list. Results are shown in Additional file 2. Genes involved in epidermal development and keratinocytes were among the most enriched GO terms of the biological processes (BP) node in all comparisons. Also KEGG pathways were analyzed by DAVID. When comparing SCC with AK, genes involved in "focal adhesion" and  "pathways in cancer" were enriched among the DEGs (Additional file 2). In the second GSEA approach, the entire dataset was analyzed with the PGSEA package ( Figure 4). In this analysis expression profiles of AKs and SCCs were compared to those of NS. PGSEA analysis demonstrated that genes upregulated by RAS were enriched in SCC compared to NS, whereas this was not seen in AK. For both SCC and AK genes positively regulated by TNF, NFκB and fumarate-hydratase (FH, known to be involved in aggressive renal papillary cancers) [27] were enriched. Genes regulated by the oncogenes BRAF and SRC were not enriched in tumors compared to NS.
In addition to GSEA, the results from the limma analysis (log 2 FC > 0.5, FDR < 0.01, Additional file 1) were subjected to oPOSSUM analysis to investigate if the DEGs were controlled by specific transcription factor (TFs). The TFs whose transcription factor binding sites (TFBSs) were found enriched are listed in Table 4. In the list TFs of upregulated genes, the REL transcription factor class (incl. RELA and NFκB1) was overrepresented. These results correspond with the  PGSEA, where genes controlled by NFκB1 were upregulated in SCC and AK ( Figure 3). The forkhead (incl. FOXF2 and FOXD1) and ZH-finger (incl. SP1 and MZF1) TF classes were overrepresented in the TFBS of downregulated genes.

Comparison of genome-wide SNP analysis and GWEA
Expression levels of both AKs and SCCs were compared to those of NS based on their chromosomal location to identify any regional chromosomal expression bias using the reb package [28]. Since at the genomic level only a  samples were compared with gene expression profiles derived from normal skin (NS, n = 13) samples and analyzed using PGSEA for the gene lists that contain genes responsive to oncogenes or for indicated pathways. The genes that show increased expression to NS for each pathway are indicated with 'up'. List of genes that show decreased expression relative to control cells for each pathway are indicated with 'down' [28]. The resulting t-statistic for each gene list was plotted (−10 < t < 10), p < 0.005); red squares represent significant number of genes in the list with increased expression in tumor samples (AK or SCC) relative to NS; blue squares represent a significant number of genes in each list with decreased expression.
few, non-recurrent chromosomal alterations were found, we compared the results of the genome-wide SNP analyses with those of the regional expression level results per tumor. However, in only very few regions where at the SNP level alterations were seen the expression level of genes in this region were affected too (data not shown).

Validation of the GWEA with QPCR
To validate microarray results, 9 DEGs were selected based on their high log 2 FC, highly significant adjusted p-value, involvement in identified pathways from the GSEA and potential biological relevance for SCC development. Their expression was quantified by QPCR analysis in all patient samples. The expression of these 9 genes was normalized for the expression of 4 reference genes, which were selected based on their stable expression in the GWEA, using geNORM. Relative differences between sample groups were confirmed for all genes (Table 5, Figure 5), although fold changes were often higher in QPCR results compared to those from the GWEA (Table 5). This is probably due to the larger dynamic range of QPCR compared to array-based analysis.

Discussion
In the present study we performed genome-wide expression analysis on SCCs, AKs and NS from OTRs to identify common pathways involved in the cellular transformation from normal skin cells to AKs and eventually SCC. Data analysis disclosed large numbers of DEGs between NS, AKs and tumor samples and identified the activation of several known oncogenic pathways in SCC, some of which were also activated in AKs. Our unique isogenic approach in this genome-wide profiling study on a well-documented sample series of OTRs, would circumvent the influence of any possible patient-specific differences on the analysis, including type of immunosuppressive drug used, and would allow identification of genes and pathways specifically involved in progression of skin carcinogenesis. During the last years, several genome-wide expression studies on SCC and AKs have been conducted [21][22][23][24]. Most studies focused on genes differentially expressed in SCCs compared to normal skin, although some also included AK samples and cell lines. Since the various studies used different gene expression platforms, tumor selection criteria and isolation and labeling methods, it was impracticable to combine data from these studies with our data into one large dataset and perform a meta-analysis. However, we do see considerable overlap in results obtained, especially after pathway en gene set analyses.
Cluster analysis showed that samples derived from the same patient did not cluster together, indicating that gene expression differences between sample groups (NS, AK and SCC) were larger than differences between patients. So, inter-patient variation did not discernibly influence the results of the expression analysis in our study. Furthermore, gene expression profiles of three AKs demonstrating (mild-to-) severe dysplasia (AK_P-11, AK_P-18 and AK_P42) showed more similarity to that of SCCs than to the other AKs ( Figure 2). Even though the numbers were too small to make any statistically reliable statement, this result suggests that the progression of AK to SCC is a continuous process.
Limma analysis demonstrated that more than 40 % of the expressed and annotated probes were differentially expressed between SCC and NS; smaller numbers of DEGs were identified between AK and NS and between SCC and AK. When the analysis was restricted to probes with a log 2 FC > 0.5 or log 2 FC > 1.0 the number of differentially expressed probes reduced dramatically, indicating that it was possible to identify subtle but significant changes in expression level in many genes. However, the large number of DEGs hampered identification of genes crucial to SCC development. Moreover, at least for the comparison of SCC with normal skin, differential gene expression might also be caused by non-specific differences owing to increased proliferation and metabolic rate of fast-growing tumor cells compared to keratinocytes in the epidermis. Similar results were obtained previously, where many transcripts showed relatively small changes in expression level between SCC and sun-exposed skin [22]. Additional gene-targeted, approaches are needed to investigate which genes/pathways are crucial for AK formation and subsequent SCC development.
Next to single gene analysis, GSEA was performed to identify specific gene sets and pathways involved in SCC development. With the online resource tool DAVID both GO terms and KEGG pathways were analyzed. From the GO terms analysis, many gene sets involved in (epidermal) cell proliferation and differentiation were enriched in SCC and AK. These gene sets include genes encoding keratins, MMPs, S100 molecules and SPRRs, all (highly) upregulated in SCC and AK compared to normal skin. According to the KEGG analysis, several cancer pathways were enriched in SCC vs. AK or SCC vs. NS, and included 'pathways in cancer' and 'small cell lung cancer'. These cancer pathways comprise many genes downstream of known oncogenes, like RAS and MYC.
Upregulation of downstream targets of RAS was also identified by PGSEA when comparing the profiles of SCC to NS, indicating that this pathway is activated in SCC. NFκB1 and TNF pathways were activated in both AKs and SCCs, suggesting that activation of these pathways is an early event in the development of AK and SCC. NFκB1 is part of the REL TF family, which also includes REL and  RELA. Gene promoter analysis by oPOSSUM showed that upregulated genes in SCC and AK were enriched for TF binding sites for REL transcription factors. These TFs have shown to be important in both skin development and carcinogenesis [29,30]. The results from our oPOSSUM analysis concurred with those of a previous study [22]. Both PSGEA and DAVID demonstrated activated RAS signaling in SCC, but not in AK. Previous studies have suggested that RAS might be activated via amplification of the gene locus or by activating mutations (reviewed by Boukamp [15]). However, both copy number alterations (SNP analysis) and oncogenic hot-spot mutations in HRAS and KRAS (data not shown) were not found in our SCCs, suggesting other molecular events to be responsible for activating the RAS signaling. These might include overactivation of the epidermal growth factor receptor (EGFR) upstream of RAS, which is known to occur in a subset of cutaneous SCCs [31,32].
The accuracy of the GWEA was confirmed by validation experiments using QPCR. Differential expression could be confirmed for all nine genes tested. These genes were selected based on their log 2 FC and adjusted p-value, but also for their presumed biological relevance for the development to SCC. One of these genes, the keratinocytespecific chemokine CCL27 was confirmed to be downregulated progressively in AK and SCC, corresponding to previous findings that CCL27 was downregulated in keratinocyte-derived tumors [33]. As CCL27 is involved in mediating cutaneous homing of T-lymphocytes [34], its downregulation could signify possible immune evasion by cutaneous SCCs. Furthermore, our GWEA results are in agreement with previous reports that CCL27 is downregulated through activation of the RAS-MAPK-signaling pathway in human skin tumors [33]. Several DEGs identified by GWEA concerned genes involved in epidermal differentiation. Upregulation of intermediate filament keratins KRT6, KRT16 and KRT17, concurred with our previous work in which we found these proteins abundantly expressed in cutaneous SCCs of OTRS, in contrast to normal skin [35]. Also a recent study by Hudson et al. with samples from immunocompetent patients, showed that genes involved in epidermal differentiation were amongst the most DEGs between SCC and normal skin [22].
Next to the genes important in epidermal differentiation, genes involved in extracellular matrix (ECM) organization were differentially expressed in SCCs compared to NS. Among others, QPCR confirmed MMP1, MMP3, MMP9 and MMP10 to be highly upregulated in SCC compared to normal skin, sometimes already in AK. Although MMPs are normally induced temporarily in response to exogenous signals, they are known to play a role in tumor progression and metastasis. MMP1, MMP3 and MMP10 are virtually absent in normal skin, but often found upregulated in cutaneous SCC, although not in AK [36,37]. Moreover, MMP9 has been shown to be involved in the early development of cutaneous SCC in immunosuppressed individuals [38]. Even though it is known that the different MMPs have separate functions [37], MMPs often cooperate in parallel and are regulated by specific and non-specific inhibitors, such as serine proteases such as chymase, plasmin and kallikrein, to achieve targeted ECM degradation. We were able to identify several kallikrein proteins upregulated in SCC compared to normal skin. Upregulation of plasmin and chymase was not found, however these might be inhibited by serpin peptidase inhibitors, such as SERPINB3 and SERPINB4. These genes, which are known to produce circulating SCC-antigens 1 and 2 respectively, were highly upregulated in SCCs in our study and have previously been reported to promote survival of SCC cells in vitro [39].
The SCCs and AKs from the OTRs investigated in this study displayed very few genomic changes. The careful selection of tumor material and the results of the GWEA and QPCR analyses, make us believe that this lack of observed changes is not caused by interference of tumor material with that of normal cells. Nor do we have any indication that many alterations remained undetected because of the Illumina Genotyping Beadchip platform we used. This platform has been successfully used in various other genomic profiling studies [40,41]. Our results from the genome-wide SNP analyses contrast reports on previous studies in which SCCs displayed highly complex karyotypes, including many gains and losses of entire chromosome arms [15,42]. Also in AKs chromosomal aberrations have been reported to occur [20]. However, more recently Purdie et al. demonstrated that SCCs could be separated in genetically distinct a The Z-score determines whether a TFBS occurs more frequently in the set of co-expressed genes compared to pre-computed background set provided by the package. b The Fisher score (one-tailed Fisher exact) is calculated to determine the probability of non-random association between the co-expressed genes and the TFBS site of interest. c The number of genes in the differentially expressed gene set that could be aligned to the pre-computed background. subpopulations related to the differentiation status of the tumor [19]. Our results coincide with the observation in this study that well-differentiated SCCs possess fewer genomic aberrations compared to moderate and poorly differentiated tumors, albeit that the absolute number of aberrations per SCC in our study was lower. Although no relationship was detected by Purdie et al. between the patient's immune status and the number of aberrations observed, eight of nine tumor with 3 or less alterations originated from transplant patients. This observation and the results of our study suggest that perhaps the suppression of the immune system in OTRs may influence the development of SCCs in these patients. This phenomenon has been reported before by Rehman et al., who demonstrated that the LOH rate in SCCs from OTRs was less than half of that in SCCs from immunocompetent patients [43]. Moreover, there could also be a direct effect of the applied immunosuppressive drugs on the molecular pathogenesis of SCCs in OTRs, which has also been suggested before by others [44]. Larger compara-tive studies with SCCs from both OTRs and immunocompetent patients will be necessary to elucidate this matter further. The results from our genome-wide SNP analysis are unexpected in view of the many changes in mRNA expression levels found in both SCCs and AKs compared to normal skin. Next to a role for immunosuppression as seen in OTRs, the many differences in gene expression observed might be attributable to epigenetic alterations, such as DNA methylation or changes in miRNA expression. The fact that many genes were differentially expressed between tumors and normal skin may be caused by demethylation of these genes or modulation of their transcription factors (e.g. (de)methylation or (de) phosphorylation or interference from miRNA).

Conclusions
With this study we were able to demonstrate vast differences in gene expression profiles to exist between SCC, AK and NS from immunosuppressed OTRs. Moreover, we found that several pathways activated in SCCs were already activated in AKs, confirming the assumption that AKs are the precursor lesions of SCCs. These pathways include the NFκB1 and TNF pathways. RAS and MYC oncogenic pathways on the other hand appear to be specifically activated in SCC. Since the drastic changes in gene expression appeared unlinked to specific genomic gains or losses, the causal events driving SCC development require further investigation. The outcome of the current study also suggests that the early activation of NFκB1 and TNF pathways in the course of SCC development may offer opportunities for targeted preventive intervention in OTRs which may counter act "field cancerisation" by eradicating plaques formed by multiple AKs.

Patient material
Approval for the present studies was granted by the Leiden University Medical Center institutional review board. Patients were selected from the group of OTRs that are regularly seen at the dermatology clinic of the Leiden University Medical Center [45]. Patients with clinically suspected SCC were informed on the study and after informed consent was obtained, fresh frozen samples were obtained from the SCC. An AK on the forearm and unexposed normal skin (NS) from the buttock were taken and peripheral blood was drawn to serve as normal internal control in genome-wide SNP analysis. All clinical diagnoses were histologically confirmed. SCCs can be histologically categorized as "well", "moderate" or "poorly differentiated" tumors, based on the degree of keratinization and cellular atypia [46]. All SCCs included in this study were classified as well-differentiated. Most AKs showed mild dysplasia, with the exception of AK_P-11 and AK_P-18, both demonstrating mild to severe dysplasia and AK_P-42 with severe dysplasia. All samples were further handled in an anonymously coded fashion, following the medical ethical guidelines described in the Code 'Proper Secondary Use of Human Tissue' established by the Dutch Federation of Medical Sciences (http://www.federa.org).

RNA and DNA isolation
RNA and DNA were isolated from SCC and AK biopsy samples that contained at least 70% tumor cells, as determined by haematoxylin and eosin stained frozen sections. To ensure that the RNA and DNA represent the same tumor cells/biopsy, the various cut sections were alternately used for RNA and DNA isolation, respectively. From the sample of unexposed NS the epidermis was removed for further processing by cryosectioning parallel to the outer surface of the skin biopsy (cut to a depth where opaque whiting of the remaining dermal surface was observed).
RNA was extracted from frozen material using the RNeasy Fibrous Tissue kit (Qiagen, Hilden, Germany), which included proteinase K treatment (10 min at 55°C) of the lysed sample in RLT-buffer and on-column DNase treatment. DNA from the tumor and skin samples was isolated using the Genomic-tips 20/G kit (Qiagen). A salting-out procedure was used to isolate genomic DNA from the peripheral blood samples [47].
Both RNA and DNA were quantified using a Nanodrop (NanoDrop technologies, Wilmington, CA) and RNA was further evaluated for degradation with the Lab-on-a-chip assay on the 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA).

Genome-wide SNP analysis
Genomic DNA from AKs, SCCs and their corresponding normal controls (peripheral blood) were subjected to genome-wide SNP analysis. 400 ng of DNA was assayed with the Infinium II Sentrix HumanHap550v3 duo Genotyping BeadChip (Illumina, San Diego, CA, USA), containing over 550,000 unique tag SNP markers. Hybridizations were performed at the Leiden Genome Technology Center according to the manufacturer's instructions. Image analysis and quality control (call rate >98%) was performed in Illumina's BeadStudio Version 3.2 software using the genotyping module. Each BeadChip was selfnormalized using information contained within the array. Further analysis and visualization was performed with the Illumina Genome Viewer in Beadstudio. First, tumor data was paired with data from their corresponding normal control. Subsequently, examination of the log R ratio (signal intensity tumor/signal intensity normal control) and B allele frequency (BAF) was analyzed for the presence of copy number variations (CNVs) and/or loss of heterozygosity (LOH) in the tumor samples. By plotting the log R ratio across a chromosome, CNVs could be visualized by an increase or decrease in baseline signal. Secondly, the BAF was plotted across each chromosome to detect regions of LOH in tumor cells.
A region with LOH was defined as a stretch of heterozygous SNPs in normal DNA but homozygous in the matched tumor DNA.

Genome-wide expression analysis (GWEA)
Gene expression profiles were obtained using HumanWG-6 v2 Expression BeadChips (Illumina). In brief, 100 ng of total RNA was converted to cDNA and subsequently labeled cRNA using the Ambion Illumina TotalPrep RNA Amplification kit (Ambion, Austin TX, USA) according to manufacturer's instructions. The labeled cRNAs were hybridized overnight to the HumanWG-6 v2 BeadChips. After washing, the BeadChips were scanned using the Illumina BeadArray Reader to measure the fluorescence intensity of each probe.
Raw data was extracted from the BeadChip data files in Illumina's BeadStudio Version 3.2 software using the gene expression module. Background subtracted data was further analyzed in R-based Bioconductor package, lumi (version 1.12.4) [48]. In lumi, the data was transformed (variance-stabilizing transformation (VST)) [49] and normalized (robust spline normalization (RSN)) [48], resulting in log-transformed normalized data. Several quality control plots of both the normalized and unnormalized data were generated in the lumi package. The R-package illuminaHumanv2.db (version 1.4.1) was used for annotation. The data were purged of genes that did not meet the detection limit (expression-detection P-value >0.01) and/ or were not annotated. This selection resulted in 15,969 probes, representing 13,848 genes that were subjected to further analysis.
Hierarchical clustering and principal component analysis (PCA) were used to cluster the samples based on their expression levels, using the bioDist package (version 1.18.0).
The limma package (version 3.2.3) [50] was used to identify differentially expressed genes (DEGs) between SCC, AK and NS. A Benjamini and Hochberg False Discovery Rate (FDR) of 1% was used as cut-off to select significant DEGs and correct for multiple testing.
Gene set enrichment analysis (GSEA) for enriched Gene Ontology (GO) terms and KEGG pathways was performed with the significantly DEGs from the limma analysis using DAVID Bioinformatic Resources v6.7 (http://david.abcc. ncifcrf.gov) [51]. Secondly, GSEA on the entire data set was performed using the parametric gene set enrichment analysis (PGSEA) and regional expression bias (reb) packages (version 1.14.0) [28].
To identify activation of transcription factors (TFs) in AKs and SCCs, the DEGs from the limma analysis were investigated for over-represented TF binding sites (TFBSs) using the online analysis tool oPOSSUM [52]. Only the vertebrate taxonomic supergroup was selected and for each comparison the top 5 TFs were selected.
Both the gene expression data and genome-wide SNP analysis have been deposited in the NCBIs Gene Expression Omnibus [53] as a SuperSeries with accession number GSE32979 (http://www.ncbi.nlm.nih.gov/geo/query/ acc.cgi?acc=GSE32979).