- Research article
- Open Access
Molecular profiling of cutaneous squamous cell carcinomas and actinic keratoses from organ transplant recipients
BMC Cancer volume 13, Article number: 58 (2013)
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).
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.
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.
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.
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 . Moreover, some studies have suggested that SCCs from OTRs display a more aggressive behavior .
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 tumor-tolerance 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 . However, the oncogenic potential of HPV in cutaneous SCCs remains unclear . In a recent transcriptome sequencing study no evidence of active viral gene expression was found in large group of SCCs .
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 ). 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 . These clusters are thought to precede tumor formation  and are more frequent in OTR compared to immunocompetent patients . However, not every SCC contains a mutated TP53 gene .
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. ). 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 .
At the chromosomal level, several studies have shown that cutaneous SCCs can display complex karyotypes with large numbers of allelic imbalances [17–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–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.
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). Side-matched 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 log2 transformed fold change (log2FC) restrictions, namely log2FC > 0.5, log2FC > 1.0 and no log2FC-restriction (Table 3). With log2FC > 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 log2FC > 2.0. The cutaneous T-cell attracting chemokine, CCL27, was the most downregulated gene in SCCs compared to normal skin (log2FC = − 4.2).
With log2FC > 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 log2FC > 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 log2FC > 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)  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 (log2FC > 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 . Since at the genomic level only a 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 log2FC, 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.
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–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 log2FC > 0.5 or log2FC > 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 . 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 .
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 ). 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 log2FC and adjusted p-value, but also for their presumed biological relevance for the development to SCC. One of these genes, the keratinocyte-specific chemokine CCL27 was confirmed to be downregulated progressively in AK and SCC, corresponding to previous findings that CCL27 was downregulated in keratinocyte-derived tumors . As CCL27 is involved in mediating cutaneous homing of T-lymphocytes , 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 .
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 . 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 .
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 . Even though it is known that the different MMPs have separate functions , 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 .
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 . However, more recently Purdie et al. demonstrated that SCCs could be separated in genetically distinct subpopulations related to the differentiation status of the tumor . 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 . 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 . 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).
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.
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 . 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 . 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 .
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 self-normalized 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) . In lumi, the data was transformed (variance-stabilizing transformation (VST))  and normalized (robust spline normalization (RSN)) , 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)  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) . 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) .
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 . 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  as a SuperSeries with accession number GSE32979 (http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE32979).
Quantitative RT-PCR (QPCR)
For 10 NS, 13 AK and 15 SCC samples sufficient amounts of total RNA were available for cDNA synthesis using the iScript™ cDNA Synthesis Kit (Bio-Rad, Veenendaal, the Netherlands) according to manufacturer’s instructions. Primers (Invitrogen, Breda, the Netherlands) were designed for nine genes of interest and four reference genes (Table 6). The reference genes were selected based on their stable expression in all samples in the GWEA. QPCR reactions were performed using the SYBR Green Supermix (Bio-Rad) on the CFX384™ real-time PCR detection system (Bio-Rad). The PCR-program included initialization (6 min at 95°C), 45 cycles of denaturation (15 s at 95°C), annealing (30 s at 60°C) and elongation (30 s at 72°C), final elongation for 1 min at 72°C and a DNA melting curve (55°C to 95°C through 0.2°C increments every 10s). All samples were tested in duplicate. Specificity of PCR products was evaluated by size in agarose gel electrophoresis followed by DNA sequence analysis. Serial dilutions of cDNA from spontaneously immortalized keratinocytes (own laboratory) and universal human reference RNA (Stratagene, Santa Clara, CA USA) were included to determine PCR-efficiencies.
The stable expression of the reference genes was validated by analyzing their QPCR data using the geNORM method (M < 1.0, CV < 0.5) in the freely available qBase software [54, 55]. QPCR data of the genes of interest were normalized based on the expression of the reference genes using the normalization factor from geNORM. Relative gene expression between the different sample groups was statistically analyzed in Excel, using a two-sided Student’s T-test assuming unequal variance. P-value ≤ 0.05 was considered significant.
LH carried out the experiments, performed statistical analyses and drafted the manuscript. SC carried out the QPCR experiments and drafted the manuscript; JNBB and HCW selected patients and collected the samples; LH, JNBB, FRG, RW, LM, CPT and HV designed the study; all authors saw and approved the final manuscript.
Cornelis P Tensen and Harry Vrieling shared senior authorship.
B allele frequency
Basal cell carcinoma
Differentially expressed gene
Differentially expressed probe
False discovery rate
Genome-wide expression analysis
Human papilloma virus
Kyoto encyclopedia of genes and genomes
log2 fold change
Loss of heterozygosity
Organ transplant recipient
Principal component analysis
(parametric) geneset enrichment analysis
Reversed transcriptase polymerase chain reaction
Squamous cell carcinoma
Tumor suppressor gene
Transcription factor binding-site
Bouwes Bavinck JN, Hardie DR, Green A, Cutmore S, MacNaught A, O’Sullivan B, Siskind V, van ver Woude FJ, Hardie IR: The risk of skin cancer in renal transplant recipients in Queensland, Australia. A follow-up study. Transplantation. 1996, 61: 715-721. 10.1097/00007890-199603150-00008.
Bordea C, Wojnarowska F, Millard PR, Doll H, Welsh K, Morris PJ: Skin cancers in renal-transplant recipients occur more frequently than previously recognized in a temperate climate. Transplantation. 2004, 77: 574-579. 10.1097/01.TP.0000108491.62935.DF.
Fortina AB, Piaserico S, Caforio AL, Abeni D, Alaibac M, Angelini A, Iliceto S, Peserico A: Immunosuppressive level and other risk factors for basal cell carcinoma and squamous cell carcinoma in heart transplant recipients. Arch Dermatol. 2004, 140: 1079-1085. 10.1001/archderm.140.9.1079.
Glover MT, Deeks JJ, Raftery MJ, Cunningham J, Leigh IM: Immunosuppression and risk of non-melanoma skin cancer in renal transplant recipients. Lancet. 1997, 349: 398-
Brash DE, Rudolph JA, Simon JA, Lin A, McKenna GJ, Baden HP, Halperin AJ, Ponten J: A role for sunlight in skin cancer: UV-induced p53 mutations in squamous cell carcinoma. Proc Natl Acad Sci USA. 1991, 88: 10124-10128. 10.1073/pnas.88.22.10124.
Ullrich SE: Mechanisms underlying UV-induced immune suppression. Mutat Res. 2005, 571: 185-205. 10.1016/j.mrfmmm.2004.06.059.
Proby CM, Harwood CA, Neale RE, Green AC, Euvrard S, Naldi L, Tessari G, Feltkamp MC, de Koning MN, Quint WG, et al: A case–control study of betapapillomavirus infection and cutaneous squamous cell carcinoma in organ transplant recipients. Am J Transplant. 2011, 11: 1498-1508. 10.1111/j.1600-6143.2011.03589.x.
Nindl I, Gottschling M, Stockfleth E: Human papillomaviruses and non-melanoma skin cancer: Basic virology and clinical manifestations. Dis Markers. 2007, 23: 247-259.
Arron ST, Ruby JG, Dybbro E, Ganem D, Derisi JL: Transcriptome sequencing demonstrates that human papillomavirus is not active in cutaneous squamous cell carcinoma. J Invest Dermatol. 2011, 131: 1745-1753. 10.1038/jid.2011.91.
Ko CJ: Actinic keratosis: facts and controversies. Clin Dermatol. 2010, 28: 249-253. 10.1016/j.clindermatol.2009.06.009.
McGregor JM, Berkhout RJ, Rozycka M, ter Schegget J, Rozycka M, Bouwes Bavinck JN, Crook T: p53 mutations implicate sunlight in post-transplant skin cancer irrespective of human papillomavirus status. Oncogene. 1997, 15: 1737-1740. 10.1038/sj.onc.1201339.
Jonason AS, Kunala S, Price GJ, Restifo RJ, Spinelli HM, Persing JA, Leffell DJ, Tarone RE, Brash DE: Frequent clones of p53-mutated keratinocytes in normal human skin. Proc Natl Acad Sci USA. 1996, 93: 14025-14029. 10.1073/pnas.93.24.14025.
Rebel H, Mosnier LO, Berg RJ, Westerman-de Vries A, van Steeg H, van Kranen HJ, de de Gruijl HJ: Early p53-positive foci as indicators of tumor risk in ultraviolet-exposed hairless mice: kinetics of induction, effects of DNA repair deficiency, and p53 heterozygosity. Cancer Res. 2001, 61: 977-983.
de Graaf YG, Rebel H, Elghalbzouri A, Cramers P, Nellen RG, Willemze R, Bouwes Bavinck JN, de Gruijl FR: More epidermal p53 patches adjacent to skin carcinomas in renal transplant recipients than in immunocompetent patients: the role of azathioprine. Exp Dermatol. 2008, 17: 349-355. 10.1111/j.1600-0625.2007.00651.x.
Boukamp P: Non-melanoma skin cancer: what drives tumor development and progression?. Carcinogenesis. 2005, 26: 1657-1667. 10.1093/carcin/bgi123.
Tsang DK, Crowe DL: The mitogen activated protein kinase pathway is required for proliferation but not invasion of human squamous cell carcinoma lines. Int J Oncol. 1999, 15: 519-523.
Ashton KJ, Weinstein SR, Maguire DJ, Griffiths LR: Chromosomal aberrations in squamous cell carcinoma and solar keratoses revealed by comparative genomic hybridization. Arch Dermatol. 2003, 139: 876-882. 10.1001/archderm.139.7.876.
Purdie KJ, Lambert SR, Teh MT, Chaplin T, Molloy G, Raghavan M, Kelsell DP, Leigh IM, Harwood CA, Proby CM, et al: Allelic imbalances and microdeletions affecting the PTPRD gene in cutaneous squamous cell carcinomas detected using single nucleotide polymorphism microarray analysis. Genes Chromosomes Cancer. 2007, 46: 661-669. 10.1002/gcc.20447.
Purdie KJ, Harwood CA, Gulati A, Chaplin T, Lambert SR, Cerio R, Kelly GP, Cazier JB, Young BD, Leigh IM, et al: Single nucleotide polymorphism array analysis defines a specific genetic fingerprint for well-differentiated cutaneous SCCs. J Invest Dermatol. 2009, 129: 1562-1568. 10.1038/jid.2008.408.
Rehman I, Takata M, Wu YY, Rees JL: Genetic change in actinic keratoses. Oncogene. 1996, 12: 2483-2490.
van Haren R, Feldman D, Sinha AA: Systematic comparison of nonmelanoma skin cancer microarray datasets reveals lack of consensus genes. Br J Dermatol. 2009, 161: 1278-1287. 10.1111/j.1365-2133.2009.09338.x.
Hudson LG, Gale JM, Padilla RS, Pickett G, Alexander BE, Wang J, Kusewitt DF: Microarray analysis of cutaneous squamous cell carcinomas reveals enhanced expression of epidermal differentiation complex genes. Mol Carcinog. 2010, 49: 619-629. 10.1002/mc.20636.
Padilla RS, Sebastian S, Jiang Z, Nindl I, Larson R: Gene expression patterns of normal human skin, actinic keratosis, and squamous cell carcinoma: a spectrum of disease progression. Arch Dermatol. 2010, 146: 288-293. 10.1001/archdermatol.2009.378.
Ra SH, Li X, Binder S: Molecular discrimination of cutaneous squamous cell carcinoma from actinic keratosis and normal skin. Mod Pathol. 2011, 24: 963-973.
Hofbauer GF, Bouwes Bavinck JN, Euvrard S: Organ transplantation and skin cancer: basic problems and new perspectives. Exp Dermatol. 2010, 19: 473-482. 10.1111/j.1600-0625.2010.01086.x.
Braakhuis BJ, Tabor MP, Kummer JA, Leemans CR, Brakenhoff RH: A genetic explanation of Slaughter’s concept of field cancerization: evidence and clinical implications. Cancer Res. 2003, 63: 1727-1730.
Sudarshan S, Linehan WM, Neckers L: HIF and fumarate hydratase in renal cancer. Br J Cancer. 2007, 96: 403-407. 10.1038/sj.bjc.6603547.
Furge KA, Chen J, Koeman J, Swiatek P, Dykema K, Lucin K, Kahnoski R, Yang XJ, Teh BT: Detection of DNA copy number changes and oncogenic signaling abnormalities from gene expression data reveals MYC activation in high-grade papillary renal cell carcinoma. Cancer Res. 2007, 67: 3171-3176. 10.1158/0008-5472.CAN-06-4571.
Gugasyan R, Voss A, Varigos G, Thomas T, Grumont RJ, Kaur P, Grigoriadis G, Gerondakis S: The transcription factors c-rel and RelA control epidermal development and homeostasis in embryonic and adult skin via distinct mechanisms. Mol Cell Biol. 2004, 24: 5733-5745. 10.1128/MCB.24.13.5733-5745.2004.
Loercher A, Lee TL, Ricker JL, Howard A, Geoghegen J, Chen Z, Sunwoo JB, Sitcheran R, Chuang EY, Mitchell JB, et al: Nuclear factor-kappaB is an important modulator of the altered gene expression profile and malignant phenotype in squamous cell carcinoma. Cancer Res. 2004, 64: 6511-6523. 10.1158/0008-5472.CAN-04-0852.
Rittie L, Kansra S, Stoll SW, Li Y, Gudjonsson JE, Shao Y, Michael LE, Fisher GJ, Johnson TM, Elder JT: Differential ErbB1 signaling in squamous cell versus basal cell carcinoma of the skin. Am J Pathol. 2007, 170: 2089-2099. 10.2353/ajpath.2007.060537.
Commandeur S, van Drongelen V, de Gruijl FR, Ghalbzouri AE: Epidermal growth factor receptor activation and inhibition in 3D in vitro models of normal skin and human cutaneous squamous cell carcinoma. Cancer Sci. 2012, 103: 2120-2126. 10.1111/cas.12026.
Pivarcsi A, Muller A, Hippe A, Rieker J, van Lierop A, Steinhoff M, Seeliger S, Kubitza R, Pippirs U, Meller S, et al: Tumor immune escape by the loss of homeostatic chemokine expression. Proc Natl Acad Sci USA. 2007, 104: 19055-19060. 10.1073/pnas.0705673104.
Morales J, Homey B, Vicari AP, Hudak S, Oldham E, Hedrick J, Orozco R, Copeland NG, Jenkins NA, McEvoy LM, et al: CTACK, a skin-associated chemokine that preferentially attracts skin-homing memory T cells. Proc Natl Acad Sci USA. 1999, 96: 14470-14475. 10.1073/pnas.96.25.14470.
Commandeur S, de Gruijl FR, Willemze R, Tensen CP, El Ghalbzouri A: An in vitro three-dimensional model of primary human cutaneous squamous cell carcinoma. Exp Dermatol. 2009, 18: 849-856. 10.1111/j.1600-0625.2009.00856.x.
Son KD, Kim TJ, Lee YS, Park GS, Han KT, Lim JS, Kang CS: Comparative analysis of immunohistochemical markers with invasiveness and histologic differentiation in squamous cell carcinoma and basal cell carcinoma of the skin. J Surg Oncol. 2008, 97: 615-620. 10.1002/jso.21006.
Tsukifuji R, Tagawa K, Hatamochi A, Shinkai H: Expression of matrix metalloproteinase-1, -2 and −3 in squamous cell carcinoma and actinic keratosis. Br J Cancer. 1999, 80: 1087-1091. 10.1038/sj.bjc.6690468.
Chebassier N, Leroy S, Tenaud I, Knol AC, Dreno B: Overexpression of MMP-2 and MMP-9 in squamous cell carcinomas of immunosuppressed patients. Arch Dermatol Res. 2002, 294: 124-126. 10.1007/s00403-002-0299-x.
Ahmed ST, Darnell JE: Serpin B3/B4, activated by STAT3, promote survival of squamous carcinoma cells. Biochem Biophys Res Commun. 2009, 378: 821-825. 10.1016/j.bbrc.2008.11.147.
Peiffer DA, Le JM, Steemers FJ, Chang W, Jenniges T, Garcia F, Haden K, Li J, Shaw CA, Belmont J, et al: High-resolution genomic profiling of chromosomal aberrations using Infinium whole-genome genotyping. Genome Res. 2006, 16: 1136-1148. 10.1101/gr.5402306.
Yachida S, Jones S, Bozic I, Antal T, Leary R, Fu B, Kamiyama M, Hruban RH, Eshleman JR, Nowak MA, et al: Distant metastasis occurs late during the genetic evolution of pancreatic cancer. Nature. 2010, 467: 1114-1117. 10.1038/nature09515.
Ashton KJ, Carless MA, Griffiths LR: Cytogenetic alterations in nonmelanoma skin cancer: a review. Genes Chromosomes Cancer. 2005, 43: 239-248. 10.1002/gcc.20183.
Rehman I, Quinn AG, Takata M, Taylor AE, Rees JL: Low frequency of allelic loss in skin tumours from immunosuppressed individuals. Br J Cancer. 1997, 76: 757-759. 10.1038/bjc.1997.457.
O’Donovan P, Perrett CM, Zhang X, Montaner B, Xu YZ, Harwood CA, McGregor JM, Walker SL, Hanaoka F, Karran P: Azathioprine and UVA light generate mutagenic oxidative DNA damage. Science. 2005, 309: 1871-1874. 10.1126/science.1114233.
Wisgerhof HC, van der Geest LG, de Fijter JW, Haasnoot GW, Claas FH, le Cessie S, Willemze R, Bouwes Bavinck JN: Incidence of cancer in kidney-transplant recipients: a long-term cohort study in a single center. Cancer Epidemiol. 2011, 35: 105-111. 10.1016/j.canep.2010.07.002.
Breuninger H, Black B, Rassner G: Microstaging of squamous cell carcinomas. Am J Clin Pathol. 1990, 94: 624-627.
Miller SA, Dykes DD, Polesky HF: A simple salting out procedure for extracting DNA from human nucleated cells. Nucleic Acids Res. 1988, 16: 1215-10.1093/nar/16.3.1215.
Du P, Kibbe WA, Lin SM: lumi: a pipeline for processing Illumina microarray. Bioinformatics. 2008, 24: 1547-1548. 10.1093/bioinformatics/btn224.
Lin SM, Du P, Huber W, Kibbe WA: Model-based variance-stabilizing transformation for Illumina microarray data. Nucleic Acids Res. 2008, 36: e11-
Smyth GK: Linear models and empirical bayes methods for assessing differential expression in microarray experiments. Stat Appl Genet Mol Biol. 2004, 3: Article3-
Huang DW, Sherman BT, Lempicki RA: Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc. 2009, 4: 44-57.
Ho Sui SJ, Mortimer JR, Arenillas DJ, Brumm J, Walsh CJ, Kennedy BP, Wasserman WW: oPOSSUM: identification of over-represented transcription factor binding sites in co-expressed genes. Nucleic Acids Res. 2005, 33: 3154-3164. 10.1093/nar/gki624.
Edgar R, Domrachev M, Lash AE: Gene expression omnibus: NCBI gene expression and hybridization array data repository. Nucleic Acids Res. 2002, 30: 207-210. 10.1093/nar/30.1.207.
Vandesompele J, de Preter K, Pattyn F, Poppe B, van Roy N, de Paepe A, Speleman F: Accurate normalization of real-time quantitative RT-PCR data by geometric averaging of multiple internal control genes. Genome Biol. 2002, 3: RESEARCH0034-
Hellemans J, Mortier G, de Paepe A, Speleman F, Vandesompele J: qBase relative quantification framework and software for management and automated analysis of real-time quantitative PCR data. Genome Biol. 2007, 8: R19-10.1186/gb-2007-8-2-r19.
The pre-publication history for this paper can be accessed here:http://www.biomedcentral.com/1471-2407/13/58/prepub
The authors would like to thank the patients for the cooperation in this study, Wim Zoutman and Coby Out for their technical assistance, Branco Misovic and Jelle Goeman for their help with the statistical analysis of the GWEA, the Leiden Genome Technology Center for their technical assistance and the pathology department is for the histopathological evaluation of the tumor samples. The work was supported by a scholarship of the LUMC for LH. The authors have no conflict of interest.
The authors declare that they have no competing interests.
Electronic supplementary material
Additional file 1: Differentially expressed probes between SCC and NS, between AK and NS, between SCC and AK, and the progression list with FDR < 0.01 and log2FC > = 0.5, including the VST transformed RSN normalized data. (XLS 3 MB)
Result DAVID analysis for the differentially expressed probes between the sample groups from the limma analysis with log
Additional file 2: 2 FC > 0.5 and FDR < 0.01. (XLS 201 kb) (XLS 202 KB)
About this article
Cite this article
Hameetman, L., Commandeur, S., Bavinck, J.N.B. et al. Molecular profiling of cutaneous squamous cell carcinomas and actinic keratoses from organ transplant recipients. BMC Cancer 13, 58 (2013). https://doi.org/10.1186/1471-2407-13-58
- Cutaneous squamous cell carcinoma
- Non-melanoma skin cancer
- Actinic keratosis
- Organ transplant recipient
- Gene expression
- Genomic profiling