Identification of differentially expressed genes in cutaneous squamous cell carcinoma by microarray expression profiling

Background Carcinogenesis is a multi-step process indicated by several genes up- or down-regulated during tumor progression. This study examined and identified differentially expressed genes in cutaneous squamous cell carcinoma (SCC). Results Three different biopsies of 5 immunosuppressed organ-transplanted recipients each normal skin (all were pooled), actinic keratosis (AK) (two were pooled), and invasive SCC and additionally 5 normal skin tissues from immunocompetent patients were analyzed. Thus, total RNA of 15 specimens were used for hybridization with Affymetrix HG-U133A microarray technology containing 22,283 genes. Data analyses were performed by prediction analysis of microarrays using nearest shrunken centroids with the threshold 3.5 and ANOVA analysis was independently performed in order to identify differentially expressed genes (p < 0.05). Verification of 13 up- or down-regulated genes was performed by quantitative real-time reverse transcription (RT)-PCR and genes were additionally confirmed by sequencing. Broad coherent patterns in normal skin vs. AK and SCC were observed for 118 genes. Conclusion The majority of identified differentially expressed genes in cutaneous SCC were previously not described.


Background
Nonmelanoma skin cancer (NMSC), encompassing both basal cell carcinoma (BCC) and squamous cell carcinoma (SCC), is the most common cancer in Caucasians, and over the last decade incidence has been increased dramatically worldwide [1]. Actinic keratosis (AK) is an early stage of SCC and approximately 10% of cases progress to SCC [2][3][4]. Ultraviolet radiation (UV) is the major risk fac-tor for this disease [5,6]. Carcinogenesis is a multi-step process indicated by a series of genes that are up-or downregulated during tumor progression. In normal vs. cancerous colon tissue only 1-1.5% of approximately 30,000 -50,000 functional genes per cell were differentially expressed resulting in 548 genes [7]. A comparable quantity was identified in breast cancer cells exhibiting 700 dysregulated genes [8]. Thus, the number of differentially expressed genes seem to be similar in different cancers. Many genes showing increased expression elevated in colon-cancer represent proteins that are considered to be involved in growth and proliferation while there were often attributed to differentiation in normal tissue. Analyzing skin-or head and neck-cell lines, genes that are associated with extracellular matrix production and apoptosis were disrupted in preneoplastic cells during SCC development, whereas genes that are involved in DNA repair or epidermal growth factors were altered at later stages [9]. Transformation of keratinocytes in response to UV-radiation was examined by microarray technology using normal human epidermal keratinocytes and SCC cell lines [10]. This study detected four clusters of differentially expressed genes in normal keratinocytes vs. skin cancer cells, which may play a role in the carcinogenic pathway. However, cell lines are often different from human tissues that has been demonstrated for both colon cancer [11] and ovarian cancer [12,13]. Thus, human cancer tissues are superior compared with cell lines analyzing differentially expressed genes. So far, only one study examined differentially expressed genes in NMSC tissue using nylon-filter DNA microarrays analyzing approximately 7,400 genes [14].
In the present study, we evaluated the different expression profile of 22,283 genes in normal skin biopsies vs. AK vs. cutaneous SCC. We focused on dysregulated genes best characterizing normal skin and NMSC (comprising both AK and SCC). Overall, 42 genes were up-regulated and 76 genes were down-regulated in skin cancer and the majority of differentially expressed genes were not described earlier.

Relative expression in normal skin and NMSC
Of the 22,283 transcripts and expressed sequence tags (EST) investigated on each oligonucleotide microarray, 118 genes were detected differentially expressed in normal skin, AK, and SCC by Prediction Analysis of Microarrays (PAM) analysis excluding 81 genes (EST and genes with the description "consensus includes...") (see Methods). For each of the 6 normal skins, 4 AK, and 5 SCC (Table 1), the relative expression of each gene was examined. We have controlled the RNA quality of each human specimen analyzing the fragment sizes, and the ratio of 5'and 3'-ends using microarray test-chips with 24 human control genes of 6 randomly selected specimens and only non-degraded mRNA specimens were used in this study.
PAM analysis was used to identify genes to classify and to best characterize normal skin, AK or SCC. The rate of misclassification on the basis of individual cross validation plots was 0% (0.0) for normal skin, 25% (0.25) for AK, 20% (0.20) for SCC, and 13% (0.13) for all three classes. The first gene list contained 200 genes best characterizing normal skin (6), AK (4), and SCC (5). Under the exclusion of EST and genes with the description "consensus includes..." (n = 81) we identified 118 dysregulated genes ( Table 2).
Hierarchical clustering was performed with the identified 118 genes based on similarities of expression levels independently of the assigned class (normal skin, AK or SCC). Gene trees display gene similarities as a dendrogram, a tree-like structure made up of branches. This nested structure forces all genes to be related to a certain level, with larger branches representing the more distantly related genes ( Figure 1). The gene CHI3L1 was detected with two      The accession number, the symbol, the description of the genes, their function, their chromosome localization, the change fold, the raw signal (mean value), and the p values of the ANOVA analysis are shown. The accession numbers in bold represent the 42 genes identified by PAM and ANOVA (p < 0.05), which were significantly differentially expressed. The symbols marked in bold represent the genes verified by quantitative real-time RT-PCR. A) Up-regulated genes (42) in non-melanoma skin cancer. B) Down-regulated genes (76) in non -melanoma skin cancer. * significant expression difference verified by quantitative real-time RT-PCR. ** positive with two different affymetrix numbers (209395_at and 209396_s_at). No., numbers; T, Tumor (AK and SCC), N, normal skin; AK, actinic keratosis; SCC, squamous cell carcinoma; n.s., not significant. Cluster map analysis of 118 genes identified by PAM of 15 different specimens resulting in two classes (9 neoplastic skin lesions and 6 normal skin) Figure 1 Cluster map analysis of 118 genes identified by PAM of 15 different specimens resulting in two classes (9 neoplastic skin lesions and 6 normal skin). Prediction Analysis of Microarrays (PAM) using nearest shrunken centroids was performed with 22,283 genes, which were present on the microarray platform (Affymetrix) to identify genes best characterizing normal skin, actinic keratosis (AK), and squamous cell carcinoma (SCC). Hierarchical clustering was performed with 118 genes identified by PAM (CHI3L1 was detected with two independent Affymetrix probes and are included twice, marked with an asterisk). Thirteen genes verified by quantitative real-time RT-PCR are marked in bold. Each color patch represents the normalized expression level of one gene in each group, with a continuum of expression levels from dark blue (lowest) to dark red (highest). The minimal set of informative genes is given by HUGO Gene Nomenclature Committee (HGNC) symbols. Group (1-9) are non-melanoma skin cancer AK (1, 5, 6, and 8), and SCC (2-4, 7, and 9) showing different expression levels compared with six cases of normal skin (10-15). Numbers 6 and 10 were specimens with pooled RNAs. 15) 10 (N-02,08, 10,13,15) 11 (N-17) 12 (N-18) 13 (N-19)  different affymetrix numbers and are included twice in the cluster map ( Figure 1). The "Condition Tree" groups samples together based on similar expression profiles by standard correlation with the GeneSpring software 6.1 resulting in two classes. The specimens with the pooled RNA (normal and AK) grouped together with the nonpooled RNA in class 1 (normal skin) and class 2 (AK and SCC) (Figure 1). In class 1, all 6 normal skin specimens grouped together, and class 2 consisted of 4 AK and 5 SCC. Thus, statistical differences in the expression levels of such genes were not detected in carcinoma in situ (AK) vs. invasive cancer (SCC). Furthermore, the pooled normal skin specimens from 5 immunosuppressed patients grouped together with 5 non-pooled normal skin specimens from immunocompetent patients ( Figure 1). Thus, the expression levels of the selected genes were independent of systemic immunosuppression.
ANOVA analysis identified 364 genes including 7 EST, which were significantly differentially expressed between normal skin, AK, and cutaneous SCC (p < 0.05). Using p < 0.15 the gene list contained 2,197 genes including 42 EST. The overall agreement rate of the identified genes by ANOVA using p < 0.05 or p < 0.15 and PAM was 36% (42 of 118) or 78% (92 of 118), respectively. To identify potentially dysregulated genes between AK and SCC, we have performed ANOVA analysis in these two groups (p < 0.05), and no gene was significantly differentially expressed in AK compared with SCC. For further analysis we have used the 118 genes that have been identified by PAM and the 42 significantly differentially expressed genes identified by both methods are highlighted in Table  2.

Verification of selected genes by quantitative real-time RT-PCR
To verify the different expression levels of mRNA measured by microarray technology, we selected 13 up-or down-regulated genes with low through high change folds from the list of differentially expressed genes. These included 9 genes with a higher (change folds by microarray analysis 1.31 ->10) and 4 with a lower expression level (change folds by microarray analysis 1.43 -2.50) in neoplastic skin lesions vs. normal skin ( Table 2). Gene specific intron-flanking primers were designed for 9 upregulated genes (RAB31, MAP4K4, IL-1RN, NMI, IL-4R, GRN, TNC, MMP1, and CDH1) and 4 down-regulated genes (ERCC1, APR-3, CGI-39, and NKEFB) ( Table 3). Unspecific PCR products were not obtained for all genes shown by electrophoresis of the PCR amplicons in agarose gels. Gene-specificity of all 13 genes was confirmed by sequencing of the PCR product of each gene. The results of the real-time RT-PCR were consistent with the microarray data ( Figure 2). All genes showed the predicted expression level either higher or lower in normal skin vs. Gene Forward (5'-3') Reverse (5'-3')   APR-3  GGT TCT GAT TTC GTC CCT GA  CAG CAT TAG CTC TCG TGT CG  CDH1  TGA AGG TGA CAG AGC CTC TGG AT  TGG GTG AAT TCG GGC TTG TT  CGI-39  CGT CAA AGG TGA AGC AGG AC  ATT ATG CTC CAG TGC CCG TA  ERCC1  GGG AAT TTG GCG ACG TAA TTC  GCG GAG GCT GAG GAA CAG  GRN  CAG TGG GAA GTA TGG CTG CT  TTA GTG AGG AGG TCC GTG GT   IL-1RN  GGA AGA TGT GCC TGT CCT GT  CGC TTG TCC TGC TTT CTG TT  IL4R CAC

Discussion
We have examined the expression levels of 22,283 genes in human biopsies of normal skin and cutaneous squamous cell carcinoma (AK, and SCC) by microarray technology. One hundred and eighteen genes were differentially expressed in normal skin vs. skin cancer and fulfilled the criteria used for PAM based cluster map analysis.
Expression profiling using oligonucleotide microarrays is a useful tool to identify tumors, to distinguish different tumor entities, and to differentiate between progressing and non-progressing neoplastic lesions [7,[15][16][17][18]. In this study, we used mRNA from skin biopsies without microdissection resulting in high RNA amounts and subsequently no amplification of the RNA transcripts were required. On the other hand, we cannot avoid a mixture of dysplastic and non-dysplastic cells in our specimens. A mixture of normal epithelial cells and tumor cells are most likely present in cancerous lesions (AK) but are unlikely in normal skin specimens. If the tumor specimen contained normal and dysplastic cells, an increased gene expression in cancer vs. normal skin or vice versa was not detected by microarray analysis. Thus, the number of differentially expressed genes detected in our study represent a subset of all differentially expressed genes in skin cancer and genes showing only low differences are most likely to be unidentified.
The number of differentially expressed genes in normal tissue vs. colon cancer and breast cancer was 548 and 700, respectively [7,8]. In our study, we detected 118 genes excluding EST best characterizing normal skin and epithelial skin cancer. These genes represent approximately 20% of all genes expected to be differentially expressed in skin cancer. So far, only one study examined the expression profile in human biopsies of NMSC and skin cancer cell lines by microarray analyzing approximately 7,400 genes [14]. Although there was only a minimal overlap between human tissue and cell lines, five genes were differentially expressed both in vivo and in vitro, namely fibronectin 1, annexin A5, glyceraldehyde 3-phosphate dehydrogenase, zinkfinger protein 254, and huntingtin-associated protein interacting protein.
Of these genes the calcium and phospholipidbinding protein annexin 5 was over-expressed (ratio 2.1) and annexin 1 showed a slightly over-expression in cutaneous SCC. In our study using another approach annexin 1 was over-expressed in AK and cutaneous SCC (change fold 1.74) showing that 1 of 5 genes (20%) differentially expressed in the study of Dooley and colleagues [14] could be confirmed. Lamin A and C showed a significant higher expression in AK and SCC compared with normal skin analyzed by immunohistochemistry [19], and these genes were also up-regulated in our study. Furthermore, enzymes of the mitochondrial chain namely cytochrome c oxidase, cytochrome b, and NADH dehydrogenase were the majority of down-regulated genes. In prostatic intra-epithelial neoplasia, a high mutation rate of NADH subunits of the respiratory chain complex I was observed similar to lung and head and neck cancer [20]. Delsite and colleagues [21] examined breast cancer cell lines and suggested that a lack of mitochondrial genes leads to increased oxidative stress, reduced DNA repair, and genetic instability. Furthermore, mitochondrial dysfunction leads to an increased production of reactive oxygene species (ROS), inhibition of apoptosis, activation of oncogenes, and inactivation of tumor suppressor genes, and thus is involved during carcinogenesis [22,23]. In our study, the majority of these enzymes was down-regulated in NMSC, indicating that mitochondrial dysfunction is possible associated with cutaneous SCC.
The reliability of the identified 118 genes by microarray technology was verified and confirmed by real-time RT-PCR analyzing 13 genes, and the following discussion is based on these genes. NMSC (CDH1, MAP4K4, IL-1RN, IL-4R, NMI, GRN, RAB31, TNC, and MMP1) CDH1 (E-Cadherin) is a representative of the classic cadherin family and a calcium-dependent cell-cell adhesion glycoprotein, mutations are correlated with a variety of cancers and loss of function may lead to cancer progression and metastasis [24]. In our study we observed an over-expression of CDH1 in skin cancer vs. normal skin, although a lower expression rate was expected. This may be due to a wrongly identified gene, a loss of function due to mutations, a post-transcriptional regulation of this gene or another mechanism in skin cancer vs. other cancers.

genes up-regulated in
MAP4K4 is a representative of the serine/threonine protein kinase family activating MAPK8/c-Jun N-terminal kinase (JNK) [25]. JNK signal transduction pathway participates in the proliferation, differentiation, and apoptosis of osteoblasts and is functionally operative in the malignant transformation of osteoblasts and the subsequent development and progression of human osteosarcomas [26]. IL -1RN (interleukin 1 receptor antagonist), IL-4R (interleukin 4 receptor), NMI (N-Myc and Stat interactor), and GRN (granulin) are genes involved in cell communications. IL-1RN is a representative of the interleukin 1 cytokine family inhibiting the activities of IL-1A and IL-1B, and modulates a variety of interleukin 1 related immune and inflammatory responses [27]. The homozygous genotype IL-1RN*2/2 of the IL-RN gene was strongly associated with early-stage gastric cancer [28]. IL-4R develop allergic reactions, modulate the function of monocytes and macrophages and has been shown overexpressed in a variety of human cancer cells in vitro and in vivo like melanoma, breast, ovarian, renal, and head and neck [29]. NMI interacts with the oncogenes C-myc and Nmyc and other transcription factors containing a ZIP, HLH, or HLH-Zip motif first isolated and characterized by Bao and Zervos [30]. In addition, NMI interacts with all Stats except of Stat2 and augments Stat-mediated transcription in response to cytokines IL-2 and IFN-gamma [31]. A novel pathogenic mechanism of the transcription factor complex NMI, BRCA1 and c-Myc is the activation of telomerase, which is a key enzyme in carcinogenesis [32]. The growth factor GRN stimulates progression and metastasis of breast cancer and is involved in a variety of cancers such as clear cell renal carcinoma, invasive ovarian carcinoma and glioblastoma [33]. In our study all 5 genes showing different functions in tumorigenesis were also over-expressed in skin cancer.
Rab31 represent a family of monomeric GTP-binding proteins and belongs to the Ras family [34,35]. Ras is a protooncogene, which is evolutionary conserved and is involved in various cancers [36], but the precise role of Ras, especially Ha-ras in NMSC is unknown. TNC is an extracellular matrix protein with anti-adhesive effects, and involved in tissue interactions during fetal development and oncogenesis. TNC was associated with breast and lung cancer [37] and was over-expressed in vulvar intraepithelial neoplasia. Matrix metalloproteinases (MMP) are involved in extracellular matrix degradation and cancer invasion [38,39]. Tsukifuji and colleagues [40] reported an over-expression of MMP-1, MMP-2, and MMP-3 in skin cancer (16 Ak, 6 AK with SCC, and 15 SCC). We detected an increased expression rate of MMP-1 and MMP-9 in skin cancer and both genes showed the highest expression rate in AK/SCC indicated by the change-folds of MMP-1 (<10), and MMP-9 (4.70) by microarray analysis. All three genes are considered to be involved in a variety of cancers, they were over-expressed in our study and may play a role in the cancerogenesis of NMSC, and thus are interesting candidates for further studies.

genes down-regulated in NMSC (ERCC1, APR-3, CGI-39, and NKEFB)
ERCC1 has a high homology with the yeast excision repair protein RAD10 [41], is reduced in testis neoplasms [42] and ovarian cancer cell lines [43]. APR-3 is considered to be involved in apoptosis and was identified using subtractive hybridization strategy in order to clone apoptosisrelated genes [44]. NKEFB encodes a representative of the peroxiredoxin (Prx) family of antioxidant enzymes and may play a role in cancer development [45]. Prx II was strongly expressed in mature endothelial cells of benign vascular tumors, whereas it was weakly or not expressed in immature endothelial cells in malignant tumors of Kaposi's sarcoma and angiosarcoma [46]. In our study these genes were also down-regulated in skin cancer, and thus were consistent with the expected expression level observed in other carcinoma.

Conclusion
In conclusion, we identified 42 genes up-regulated and 76 genes down-regulated in cutaneous squamous cell carcinoma (AK and SCC) vs. normal skin, which represent approximately 20% of the genes differentially expressed in skin cancer. The majority of genes which known functions in other cancers was consistent with our results of differentially expressed genes in NMSC. These 118 genes either individually or more likely together or a subset of these genes may prove useful for diagnostic approaches.

Patients
Biopsies were obtained from 5 organ-transplanted (TX) recipients (3 kidney, 1 heart, and 1 liver, 58-73 years, median age 66 years) each normal skin, AK, and SCC. The time since transplantation ranged from 2 through 23 years (median 11 years), and no rejection was observed. All patients had multiple NMSC, such as AK, SCC and/or basal cell carcinoma, and lesions were mainly located on sun-exposed areas. The specimens from TX recipients of 5 normal skin and 2 AK specimens were pooled due to the low RNA amount of the individual specimens that was not sufficient for further microarray analyses. Furthermore, we have included 5 normal skin specimens from age-matched non-immunosuppressed individuals (17-74 years, median age 61 years). Thus, we have analyzed 6 normal skin, 4 AK, and 5 SCC specimens by microarray technology (Table 1). All clinical specimens were collected under standardized conditions by the same clinician (TF). From each organ-transplanted patient, punch biopsies (diameter 4 mm) of normal tissue, AK, and SCC were collected. Half of the tissue was transferred to liquid nitrogen within 2 minutes of resection and stored at -70°C until RNA isolation was performed. The other half of each biopsy was fixed in formalin, embedded in paraffin and sections were stained with hematoxylin and eosin for histological evaluation. All clinical diagnoses, normal skin, AK, and SCC were confirmed by histology.
The same 15 RNA specimens (or representative subsets) were used for quantitative real-time reverse transcription (RT)-PCR of 13 selected genes for verification ( Figure 2

RNA isolation and microarray hybridization
Total RNA was isolated using a modified RNeasy Micro Kit protocol (Qiagen, Hilden, Germany). The modification included the homogenization of the frozen tissue in 300 μl of buffer RLT (Qiagen) with 20 ng Glycogen (Roche, Mannheim, Germany) using a rotar-stator homogenizer "Ultra Turrax T25" (Janke & Kunkel, Staufen, Germany). The homogenized tissue was digested with 10 μl Proteinase K (10 mg ml -1 ) (Roth, Karlsruhe, Germany) at 55°C for 15 min. Subsequently the sample was digested with DNase I (Invitrogen, Karlsruhe, Germany). Quantification of isolated RNA was performed using UV-spectroscopy and the quality was determined both by A 260 /A 280 ratio and Agilent bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). Five microgram total RNA was used for cDNA synthesis with 5 pmol μl -1 T7-oligo(dT) 24 primer and was performed at 43°C for 90 minutes with the "Superscript First-Strand Synthesis-System" for RT-PCR (Invitrogen). Second-strand synthesis was performed with complete cDNA. The cDNA solution was incubated at 16°C for 2 hours followed by an incubation step for 20 min with 6 U T4-DNA polymerase at 16°C and the reaction was stopped using 10 μl of 0.5 M EDTA. The double stranded cDNA was purified by phenol/chloroform, ethanol precipitated and the pellet was resuspended in 12 μl of DEPC water. Labeled cRNA was generated from the cDNA sample by an in vitro transcription reaction that was supplemented with biotin-11-CTP and biotin-16-UTP (Enzo Diagnostics, Farmingdale, NY, USA) according to the manufacturer. The cRNA was quantified by A 260 , and the quality was determined using the labchip bioanalyzer (Agilent). Only cRNA specimens with a high quality were selected for further analyses. Fragmented cRNA (15 μg) was used to prepare 300 μl hybridization cocktail (100 mM MES, 1 M NaCl, 20 mM EDTA, 0.01% Tween-20) containing 0.1 mg ml -1 of herring sperm DNA, and 0.5 mg ml -1 acetylated bovine serum albumine. Control cRNA was used in order to compare hybridization efficiencies between arrays and to standardize the quantification of measured transcript levels and was included as component of the 'Eukaryotic Hybridization Control kit' (Affymetrix, Santa Clara, CA, USA). The cocktails were heated to 95°C for 5 minutes, equilibrated at 45°C for 5 minutes, and clarified by centrifugation. The cocktail was hybridized to HG U133A arrays (Affymetrix) at 45°C for 16 hours. The arrays were washed and stained with a streptavidin-conjugated fluor using the GeneChip fluidics station protocol EukGE-WS2 (Affymetrix) according to the manufacturer's instructions. Arrays were scanned with an argon-ion laser confocal scanner (Hewlett-Packard, Santa Clara, CA) with detection at 570 nm. Data were extracted using Microarray Suite version 5.0 (Affymetrix) and linearly scaled to achieve an average intensity of 2,500 per gene. Text files were exported to determine the intensity of each interrogating oligonucleotide perfect match probe cells or mismatch probe cells. In addition, the ratios of 5'-and 3'-ends of mRNA were analyzed of six randomly selected specimens (two of each group) using microarray test-chips (Test3 Array) containing 24 human housekeeping/maintenance genes (Affymetrix) and RNA degradation was not observed.

Bioinformatic analysis
The Data Mining Tool 3.0 (Affymetrix) and GeneSpring software package 6.1 (Silicon Genetics, Redwood City, CA, USA) were used for different replicates and statistical analyses were performed in order to compare between cancer stages. For each hybridization, the intensities were normalized in three steps, (1) data transformation, (2) per chip, and (3)  All processings from raw data, normalized raw data and pdetection values of the microarray experiments to the final tables and figures and a description of the method are provided as supplemental material with the series number GSE2503 [47]. Thus, the entire process of analysis is completely transparent and the description of the methodology used is according to the MIAME standard (minimum information about a microarray experiment).
We used PAM for classification of tumors and identifications of genes that were significantly different expressed between three groups. PAM is a statistical technique for class prediction from gene expression data using nearest shrunken centroids. The technique has advantages in accuracy, especially when more than two classes are considered to be examined [48,49] as it is required for this study. PAM ranks genes using a panelized t-statistic and uses soft-thresholding to identify a gene set for classification. Data analysis was performed with 22,283 genes of all 15 specimens depending on their class (normal skin, AK, or SCC). The number of genes used was controlled by a thresholding parameter, which was determined with a 10fold cross-validation. We used the imputation engine method with the k-nearest neighbor (n = 10), and the threshold 3.5 was chosen to minimize the overall error rate. This cross validation also allows a judgment of the classification quality. For the detailed mathematic procedure, see Tibshirani et al. [48].
In addition, we have independently applied the ANOVA model using two different p-values (p < 0.05 and p < 0.15) to identify dysregulated genes between three groups (normal skin, AK, and SCC) and two groups (AK and SCC) to focus on differences between these two groups. Multiple testing corrections were performed by the false discovery rate of Benjamini and Hochberg for all analyses.
Hierarchical clustering was performed with the genes best characterizing normal skin and NMSC identified by PAM. Genes of all 15 specimens with different expression profiles were grouped by standard correlation with Gene-Spring software package 6.1. Hierarchical clustering of the genes was based on similarities of expression levels.

Quantitative real-time RT-PCR
Real-time RT-PCR with the LightCycler system (Roche) was used as an independent method to validate the microarray expression data and to assess quantitative gene expression. Thirteen genes were selected including 9 upregulated and 4 down-regulated genes in NMSC with low, moderate, and high change folds. In addition, 3 of the 13 genes (MMP1, RAB31, and TNC) were verified with an increased number of different specimens of immunosuppressed and immunocompetent patients (22 normal skin, 11 AK, and 15 SCC) (see Methods section 'Patients'). RT was performed with the "Superscript First-Strand Synthesis-System" (Invitrogen) using oligo-dT as described by the manufacturer. The concentration of cDNA was quantified with "OliGreen ssDNA Quantitation Kit" (Molecular Probes, Leiden, Netherlands). Specific PCR primers for the target genes were designed using the Primer3 software program [50], and synthesized by Metabion (Planegg-Martinsried, Germany). The primers of each gene were located in different exons to exclude DNA contamination. Amplification mix (20 μl) contained 20 ng of cDNA, 500 nM of each primer, 2 μl LightCycler FastStart Reaction Mix Syber Green I (Roche), 3 mM MgCl 2 and sterile double distilled water. The concentration of MgCl 2 varied depending on each specific primer pair between 3-5 mM. PCR reaction was initiated with 10 min denaturation at 95°C followed by 40 cycles (95°C for 10 sec, 60°C for 5 sec, and 72°C for 10 sec). Fluorescence detection was performed immediately at the end of each annealing step and the purity of each amplification product was confirmed by generating melting curves. All specific RT-PCR products of target genes were purified by gel extraction and confirmed by sequencing with gene specific primers (Table 3) using the DNA sequencing kit and the ABI PRISM 310 Genetic Analyzer (Applied Biosystems, Foster City, USA). A negative control without reverse transcriptase was included in each PCR experiment. The expression of RPS9 was used to control equal RNA loading and to normalize relative expression data for all other genes analyzed. The copy ratio of each analyzed cDNA was determined as the mean of two experiments. The U-Test of Wilcoxon, Mann, and Whitney was applied for estimation of differentially expressed transcripts identified by real-time RT-PCR. A pvalue < 0.05 was considered significant for alpha.