Gene Variability between Perineural-Positive and Perineural-Negative Squamous Cell Skin Cancers

Background: Recurrent cutaneous squamous cell cancer (CSCC) is associated with poor outcomes with perineural invasion reported as a frequent finding in such lesions. Given the morbidity associated with late recurrence, identifying aggressive subtypes of CSCC at the time of primary excision is all the more necessary. This project sought to identify differentially expressed genes (DEGs) between perineural invasion-positive (PP) and perineural invasion-negative (PN) CSCC. Gene-based classification models for diagnosis of perineural invasion in CSCC were also developed. Method: Forty fresh-frozen surgical specimens of CSCC with presence or absence of histopathological perineural invasion were processed for RNA isolation and hybridization to Affymetrix-U219 DNA microarrays. Raw gene expression data were normalized by Robust Multi-array Averaging (RMA) and log2 transformed. DEGs were identified by empirical Bayes statistics using the Bioconductor limma package. BRB-ArrayTools software was used to develop gene expression-based sample classification models. Using leave-one-out cross-validation, the resulting accuracies of eight different classification algorithms were evaluated. Results: Twenty-one PP and 19 PN samples were analyzed. At a stringent limma p-value (p<0.001), 24 genes were differentially expressed between specimens. The cross-validated performance of the eight classification models exhibited a mean accuracy of 85-95%. Diagonal linear discriminant was most accurate at 95%, followed by Bayesian compound covariate at 94%. The poorest accuracy (85%) was observed for 1-Nearest neighbor and Support vector machines. For all eight methods, the sensitivities and specificities ranged from 79%-95%. Conclusion: Gene expression distinguishes between PP and PN CSCC. Classification models based on these gene patterns distinguish PP and PN cancers with strong statistical accuracy and may potentiate more timely and objective diagnosis of perineural invasion that could guide more comprehensive adjuvant therapies.


Introduction
Cutaneous squamous cell carcinoma accounts for approximately 20% of new cases of non-melanoma skin cancer diagnosed each year [1]. Some 2.5-14% of these cancers have perineural invasion and with simple excision [1], have been found to recur 46% of the time [2]. Many who present with clinical signs of perineural invasion, such as paresthesias or facial weakness, have a history of previous CSCC resection with reported negative margins [3,4]. These recurrences often require skull base or sub-cranial resections in order to obtain negative margins. The morbidity of such resections and the subsequent adjuvant radiotherapy is certainly much greater than a local excision at the primary site. Identifying lesions with perineural invasion at the time of primary excision could direct more extensive surgical management at that time or radiotherapy to limit recurrences and poor outcomes.
Molecular markers of perineural invasion have been investigated in multiple cancer types; however no concensus has been determined for CSCC of the head and neck. Markers in the head and neck, including MT1-MMP [5], p75 nerve growth factor receptor [6], tyrosine kinase receptor 4 [7][8][9], foxp3 [10], alpha heat shock protein [11], Ecad proteins [12], CD44 [13], and NGF and TrkA overexpression in adenoid cystic carcinoma [14] have all been identified as associated with perineural invasion, however there has been limited replicability. To date there have been no gene expression profiling studies aimed at profiling cancer tissue with perineural invasion to find specific genes implicated in the disease. This study aims to identify gene signatures linked to perineural invasion using microarray analysis in order to build an initial framework of genes with the hope of better honing these signatures to a more specific signature in future studies. This signature could potentially be used to create a diagnostic or screening assay to be used in the clinical setting for identifying samples that are high risk for perineural invasion.
Cancer Center of Wake Forest University) were searched for cases of CSCC of the head and neck after Institutional Review Board approval. All samples within the tumor bank had previously been placed into the bank after patient consent. Fresh frozen samples were utilized for microarray expression analysis. Additional samples were collected in a prospective fashion by identifying patients in our outpatient head and neck oncologic clinics, as well as a private dermatologic surgery clinic. These patients were consented for use of their tissue samples in this study.

Confirmation of tumor
All eligible samples underwent review by a faculty or fellow dermatopathologist in order to confirm the presence or absence of perineural invasion and to assess the adequacy of cancer tissue versus benign stromal tissue. Perineural invasion was determined histologically by using hematoxylin and eosin stain. Tumor samples were macrodissected (1 to 4 cores, approximately 2 mm 2 , were punched from each block) to obtain tissue with >80% cellularity of malignant epithelium. Tumors with sparse cancer cellularity were excluded from analysis.

RNA processing and quality analysis
Tumor tissues samples were disrupted under liquid nitrogen and homogenized according to the Qiagen QiaShredder protocol. Total RNA (approximately 1-10 ug per sample) was isolated from each tumor sample according to the Qiagen RNeasy Microarray Tissue Mini Kit protocol. Once adequate mass was confirmed using an Eppendorf BioPhotometer, the quality of the RNA was assessed using an Agilent 2100 Bioanalyzer. RNA suitable for microarray processing were selected using the following criteria: 1) RNA integrity number (RIN) >8.0; and 2) absorbance ratio (A260/A280) between 1.8 -2.1.

Expression microarray analysis
Forty samples [21 perineural-positives (PP) and 19 perineuralnegatives (PN)] that met tissue and RNA quality criteria were profiled on the Affymetrix GeneAtlas U219 human genome array strip in the Comprehensive Cancer Center's Cancer Genomics Shared Resource. For each sample, 250 ng of total RNA was used as a starting template for the reverse transcription and labeling reactions. Exogenous PolyA controls were spiked into each sample to monitor quality of amplification reaction. The amplified and biotinylated cRNA targets were hybridized to the microarrays, stained, washed and scanned according to standard Affymetrix protocols. Log intensity distributions and pair-wise correlations between arrays were examined to assess quality of each microarray hybridization. The resultant raw data (CEL files) were normalized using the RMA (Robust Multi-array Average) algorithm [15] as implemented in the Affymetrix GeneAtlas Instrument Control Software and log2 transformed for downstream statistical analyses.

Analytical plan
Differentially expressed genes (DEGs) between PN and PP samples were identified by using random variance model for univariate tests implemented in BRB-Array Tools [16]. At two different significant thresholds of p values 0.001 and 0.005, we obtained 24 and 130 probe sets, respectively. In addition, we used Significant Analysis of Microarray (SAM) with target proportion of false discoveries of 0.05 to obtain 134 probe sets. The above selected sets of genes were used in pathway enrichment analyses using the NIAID's DAVID microarray resource [17] and construction of classification models.
BRB-Array Tools software [16] was used to develop gene expression-based sample classification models. Using leave-oneout cross-validation, the resulting accuracies of eight different classification algorithms were evaluated. They are (1) The Prediction Analysis for Microarrays (PAM), using the shrunken centroid algorithm [18]; (2) The Compound Covariate Predictor, using a weighted linear combination of log-intensities for genes that are univariately significant at the specified level [19]; (3) The Diagonal Linear Discriminant Analysis (DLDA), using linear discriminant analysis, but ignoring correlations among the genes to avoid overfitting the data [20]; (4-5) The 1 or 3 Nearest Neighbor Predictor, using Euclidean distance as the distance metric to predict the class of test samples; (6) Nearest Centroid Prediction, predicting a test sample belonging to a class corresponding to the nearest centroid; (7) Support Vector Machine (SVM), with linear kernel functions, to separate the data subject to penalty costs on the number of specimens misclassified [21]; (8) The Bayesian compound covariate predictor, using weighted average of the log expression values of the selected genes, with the weights being the t statistics of differential expression in that training set [22]. For each outcome of the eight methods, accuracy=[(A+B)/n], sensitivity=[A/(A+C))] and specificity=[B/(D+B)] were calculated. A represents the number of perineural samples predicted as perineural, B presents the number of non-perineural samples predicted as nonperineural, C represents the number of perineural samples predicted as non-perineural, D represents the number of non-perineural samples predicted as perineural, and N represents the total number of samples.

Results
Twenty-one PP and 19 PN samples were analyzed and deemed adequate for analysis. At a stringent p-value threshold of p<0.001, 24 genes were differentially expressed between PP and PN specimens (Table 1). At a p-value threshold of p<0.005, 130 genes were differentially expressed between PP and PN specimens. Gene expression-based classification models were developed using both p value thresholds, p<0.001 and p<0.005. At the p-value threshold of p<0.001, the crossvalidated performance of the eight classification models exhibited a mean accuracy of 85-95%. DLDA was most accurate at 95%, followed by Bayesian compound covariate at 94%. The poorest accuracy (85%) was observed for 1-Nearest neighbor and Support vector machines. For all eight methods, the sensitivities and specificities ranged from 79%-95% (Table 2).
At the p-value threshold of p<0.005, the cross-validated performance of the eight classification models exhibited a mean accuracy of 82-92%. PAM was most accurate at 92%, followed by the Compound covariate predictor, DLDA, and Bayesian compound covariate predictor 90%. The poorest accuracy (82%) was observed for 1-nearest neighbor, 3-nearest neighbors, and Support vector machines. For all eight methods, the sensitivities and specificities ranged from 63%-90%. We performed gene ontology/pathway enrichment analysis on selected DEGs at both p-value thresholds of p<0.001 and p<0.005, however, no significant enrichments were observed.

Discussion
Perineural invasion of the head and neck has been identified as a marker of poor outcomes, decreased survival, increase locoregional recurrence and shorter time to recurrence [4]. In the head and neck, cancer cells may spread along the entire cranial nerve network in a retrograde or anterograde fashion in a relatively low resistance manner. Once patients become symptomatic, many times with facial paresthesias or facial paresis, the tumor has often spread into larger named nerves or to the skull base, making surgical options more morbid. Though new excision techniques, such as Mohs micrographic surgery, allow serial frozen sectioning to track margins, these techniques generally treat smaller and less advanced skin cancers with more subtle perineural invasion [1]. Artifact from frozen sectioning has contributed to the underreporting of perineural invasion with this technique, despite its well-documented ability to clear microscopic disease [23]. Therefore, diagnostic or screening methods utilizing gene markers of perineural invasion or more aggressive subtypes could decrease underreporting,  regardless of excision technique, and identify high-risk cancers at an earlier stage.
A discussion of high-risk CSCC would not be complete without mentioning other histologic and clinical factors that have been identified as influencing clinical aggressiveness of CSCC. Fitzpatrick skin types II and III, immunosuppression as is often encountered with organ transplantation, cigarette use, and prior history of nonmelanoma skin cancer are clinical factors that have been found to predict more aggressive CSCC [24]. Lymphovascular invasion, poorly differentiated histology, diameter of 2 cm or greater, depth beyond the dermis, and location on ear or lip are tumor-specific factors that have been associated with increased risk of recurrence [25]. Despite our goal to specifically distinguish between PP and PN samples, it is important to remember that the physiology behind recurrence and poor outcomes in this patient population is multi-factorial.
To our knowledge, expression profiling studies seeking to identify genes involved in perineural invasion have not been reported. However, markers of aggressive subtypes and progression of squamous cell carcinoma have been documented. Kivisaari et al used tissue microarrays and immunohistochemistry in order to find MMP-7 up-regulation in epidermolysis bullosa-associated CSCC, known to be an aggressive subtype, suggesting a potential therapeutic effect of epidermal growth factor receptor antagonists in treatment of advanced CSCC [26]. Farshchian et al found elevated expression of Serpin Peptidase Inhibitor Clade A Member 1 (Serpin A1) in more aggressive subtypes of CSCC [27]. Nathan et al investigated markers of the MEK/ERK pathway in metastasis in patients undergoing elective parotidectomy with CSCC. They found increased expression of pS6 in CSCC with parotid metastasis compared to those without parotid metastasis. Though this work is useful in helping to differentiate between aggressive subtypes, we wanted to look at specifically perineural invasion and the gene markers associated with this process.
This study is unique in that we used DNA microarrays to examine the differential expression patterns of over 20,000 genes in PP and PN tissues. By testing classification models whose gene features were selected based on multiple significance thresholds, we found that the threshold of p<0.001 (24 gene set) demonstrated the best predictive performance overall, likely owing to the statistical power of using the top 24 statistically significant genes versus the top 130. Of these 24 genes, none could be identified in the literature as having previously been associated with perineural disease in CSCC. By using leave-oneout cross validation to control for over-fitting, we were able to develop classification models with high accuracy with a range of 85-95%. In developing these classification models we found that these 24 genes have a robust ability to distinguish PP and PN tissues regardless of the classification algorithm used, though some perform better than others.
Of the 24 genes comprising our classifiers, no underlying biological processes or pathways could be discovered by gene ontology enrichment analysis. This may be due, in part, to the fact that a number of these genes lack functional characterization, while the remaining functionally characterized genes represent relatively diverse biological functions with few direct implications in cancer. One intriguing observation, however, that could shed light on the pathological differences between PP and PN, is the >2-fold reduction of TXNIP expression observed in PP samples. TXNIP is a recently identified tumor suppressor gene [28] that negatively regulates energy metabolism by inhibiting uptake and utilization of glucose [29,30]. In this capacity, down-regulation of TXNIP expression could conceivably contribute to the more aggressive clinical behavior of perineural invaison through augmentation of energy metabolism, though whether TXNIP inhibition enhances perineural invasion and subsequent recurrence is unknown and warrants further study.
A weakness of this study was the small sample size analyzed. As the goal of this project was to create pilot data that could serve as groundwork for future studies, we estimated that forty samples would be adequate to identify a number of lead genes capable of discriminating PP and PN. We believe the findings presented here provide a strong rationale for subsequent and larger population-based studies aimed at honing and validating the prognostic performance of the 24-gene classifier towards identifying PP disease at the time of primary excision.

Conclusion
Our work shows that gene expression patterns can distinguish PP and PN CSCC. Internally cross-validated classification models based on these gene patterns distinguish PP and PN cancers with high sensitivity and specificity. Gene-based classification models may potentiate more timely and objective diagnosis of perineural invasion and may have utility in guiding more comprehensive adjuvant therapies.

Author Contributions
ACM designed the study, acquired samples, was primarily responsible for the study, and participated in manuscript preparation and final approval. LM designed portions of the study, in particular the methods, participated in manuscript preparation and final approval. JC designed portions of the study, in particular the data analysis, participated in manuscript preparation and final approval. JDB participated in the conception of the study, assisted with acquiring funding, acquired samples, participated in manuscript preparation and final approval. AC acquired samples, personally performed all tissue work, and participated in manuscript preparation and final approval. All authors have read and approved the final manuscript.