Gene Expression Profiles Identify Epithelial-to-Mesenchymal Transition and Activation of Nuclear Factor-κB Signaling as Characteristics of a High-risk Head and Neck Squamous Cell Carcinoma

Introduction

Head and neck squamous cell carcinoma (HNSCC) is the fifth most common cancer in the United States, with the main risk factors being tobacco and alcohol use ( 1). Unfortunately, the majority of patients present with advanced disease stages that require aggressive therapy. Because of the debilitating nature of aggressive therapy and the limited ability to identify patients at high-risk for treatment failure, better biomarkers of prognosis as well as more effective and less toxic treatment for the high-risk patients are desperately needed.

With the development of DNA microarray technology, it is now possible to better predict survival using gene expression profile of the primary tumors independent of tumor-node-metastasis (TNM) staging ( 2– 6). Previously, we used gene expression profiling to identify a group of HNSCC patients with a high-risk of recurrence using a 582 intrinsic gene set ( 2). The patients in the high-risk group (labeled group 1 in ref. 2) had worse recurrence-free survival (RFS) compared with other patients (P = 0.02, log-rank test). However, the widespread clinical use of array-based gene expression profiles has been limited by the availability of fresh frozen tumor specimens and the relatively short follow-up information. The ability to assay RNA obtained from formalin-fixed tissue would be a great advancement that allows for analyses of large existing tissue collections with sufficient clinical information that could be correlated with clinical outcomes. We therefore determined feasibility of microarray gene expression analyses of formalin-fixed tissue and compared the generated predictive gene lists with an independent frozen tumor set generated in our previous study ( 2).

By combining the gene expression data sets from the formalin-fixed and previously analyzed frozen tumor tissues as training and testing sets, we identified a 75-gene list that is highly predictive of RFS as a reliable prognostic biomarker. In addition to the focus of biomarker discovery, we also wanted to understand the molecular characteristics of the high-risk tumors by studying the predictive genes individually; however, the gene list did not obviously implicate a single molecular pathway because the selected genes represented various cellular functions. To determine predominant biological process, we interrogated the expression data using Gene Set Enrichment Analysis (GSEA; ref. 7) and have shown that the gene sets involved in epithelial-to-mesenchymal transition (EMT), nuclear factor-κB (NF-κB) activation, and cell adhesion are significantly enriched in the high-risk tumors. Our data suggest that EMT, NF-κB, and cell adhesion pathways are important targets and provide impetus to test currently available anticancer agents targeting the pathways for the high-risk HNSCC patients.

Materials and Methods

Patient selection and specimen collection. Forty samples, including 33 formalin-fixed tumor blocks, 6 frozen tumors (matching 6 of the formalin-fixed tumors), and 1 formalin-fixed normal mucosal epithelium, from 29 patients with primary HNSCC of the oral cavity, oropharynx, and larynx who were treated between 1992 and 2005 were obtained following Institutional Review Board approval from Vanderbilt University Medical Center, Veterans Administration Medical Center (Nashville, TN) and The University of Texas M.D. Anderson Cancer Center (Houston, TX). We specifically included specimens that had been in storage for longer periods to study RNA integrity over time.

RNA isolation and DNA microarray analyses. H&E-stained slides for each tumor were examined, and areas with >70% tumor cellularity were chosen for macrodissection. RNA isolation and amplification were done using Arcturus Paradise kit (Arcturus, Mountain View, CA) as suggested by the manufacturer. The quality of the RNA was confirmed using Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA). The amplified RNA (aRNA) was labeled with ENZO BioArray High Yield RNA Transcript Labeling kit (Affymetrix, Santa Clara, CA), and biotin-labeled aRNA (15 μg) was fragmented. This aRNA was loaded on to the Affymetrix Human Genome X3P GeneChip and processed according to the manufacturer's recommendations. The raw microarray data were uploaded to the Vanderbilt Microarray Shared Resource (VMSR) database.

Confirmation of expression data by real-time PCR. Total RNA from the formalin-fixed tumors was available for real-time PCR analysis. Total RNA (50 ng) was amplified using the NuGEN WT-Ovation RNA Amplification kit (NuGEN, San Carlos, CA; ref. 8). The amplified cDNA was cleaned using the Qiagen PCR purification kit (Qiagen, Valencia, CA). Amplified cDNA (25 ng) was used per reaction, and the probes were obtained from Applied Biosystems (Foster City, CA). Four genes were analyzed by real-time PCR: KRT14, ACTN1, COL5A1, and PLEC1 using Applied Biosystems Taqman FAM-labeled probes. The endogenous genes B2M, UBC, and Eukaryotic18S were used as internal calibration standards. The average of these three internal genes was used to normalize the real-time PCR results from the set of four genes. Analysis of each sample was done in triplicate on an Applied Biosystems 7900HT instrument (Applied Biosystems). Real-time PCR data were analyzed by the 2^−ΔΔCT method as described previously ( 9, 10).

Processing of the data set generated on Affymetrix X3P GeneChips. The raw data containing ∼61,000 probe sets (representing ∼47,000 transcripts) were normalized using GCRMA across the 40 CEL files and filtered for expressed genes by eliminating measurements with <5 in log 2 space and retaining genes that are present in >60% of the samples. K-nearest neighbor was used to impute the missing values. Principal component analysis was applied to detect any bias introduced by sample collection or processing. A bias in the data surfaced in the top two principal components. Segregation of the sample in these components seemed to be due to the age of the blocks; however, other factors, such as the time between the tissue collection and formalin fixation or the length of fixation, may also be involved, which are difficult to specify. The bias was eliminated by doing singular value decomposition to adjust for the bias ( 11, 12). Independent validation by both molecular classification and survival predictor was improved after this adjustment.

Molecular classification by intrinsic gene set analysis. The intrinsic gene set was determined as described previously ( 2). Briefly, eight pairs of tumors (six formalin-fixed and frozen pairs and two formalin-fixed pairs) from the same patients were used to generate a set of genes that are highly variable across the tumors but consistent within the pairs. Thus, only genes that are intrinsic to the tumors are obtained. Using this intrinsic gene set, all 40 samples were grouped using hierarchical clustering for molecular classification. To test the validity of the intrinsic gene set, previously published 60 frozen HNSCC samples analyzed on Agilent cDNA microarray were used as an independent test set. UniGene IDs (build #184) were assigned to each gene that was present on both X3P GeneChips and Agilent cDNA microarrays. Genes with identical UniGene identifiers were collapsed by taking the median value within a sample. The data points were quantile normalized across the two platforms. The UniGene identifiers corresponding to the intrinsic gene set from the X3P GeneChip data set were extracted. The data for each gene were median centered and analyzed with hierarchical clustering ( 13). Finally, the clustering was visualized by TreeView.

Determination of the predictive gene list for high-risk patients. To identify the minimum number of genes that can predict the high-risk tumors, 28 formalin-fixed samples with recurrence and survival data were trained using supervised principal components analysis for the identification of high-risk versus low-risk groups ( 14). The predictive high-risk gene list was tested on the 60 frozen tumor data set after mapping the genes from two platforms as described above. In addition, using high- and low-risk as a classifier, the deregulated pathways that are differentially enriched in the two groups were determined using GSEA, which allows interrogating the gene expression data by 1,325 a priori–defined sets of genes ( 7).

Statistical analyses for RFS. RFS time was determined as the time from diagnosis to recurrence, to disease-specific death, or to the last follow-up date. Univariate survival analysis was done using log-rank test using the SAS/STAT software package (SAS Institute, Research Triangle Park, NC) and plotted as Kaplan-Meier curves.

Results

Patient characteristics. Forty samples from 29 patients with HNSCC of oral cavity, oropharynx, and larynx were collected between 1992 and 2005 and analyzed for gene expression using DNA microarrays. Matched frozen tumor specimens were available for six of the formalin-fixed tumors. We did replicate or triplicate analysis on three of the formalin-fixed tumors to study reproducibility of the assay. We also analyzed a formalin-fixed normal oropharyngeal epithelium specimen. Detailed patient characteristics are described in Table 1 . Twenty-eight of 29 patients had known recurrence and survival data with the median follow-up of 33 months. Fourteen of 28 patients had recurrence after completion of the curative treatment with the average time to recurrence of 13.5 months.

RNA isolation from formalin-fixed tissue and DNA microarray hybridization. The total RNA yield from a single, macrodissected, 7-μm section of formalin-fixed paraffin-embedded tissue with a tumor surface area of 5 to 10 mm² ranged from 67 to 2,100 ng. Using 50 to 100 ng input total RNA for the linear amplification, the yield of aRNA ranged from 10.2 to 50 μg. Approximately 10% of the RNA isolations failed, and isolation had to be repeated; however, once enough RNA was isolated, all amplification reaction produced RNA of sufficient quality for hybridization to the DNA microarray. Compared with frozen tissues, there were less present calls on the hybridized DNA microarray from formalin-fixed tissues (∼8,000 versus 25,000 probes with present calls) determined by MAS5.0 output. A subtle bias, possibly due to more RNA degradation for older tissue blocks, was detected by Principal Component Analysis in the array data from samples with different age. However, the number of present calls on each array was not correlated with the age of the blocks. In addition, the bias was easily corrected using standard statistical methods that have been developed and used to correct batch effects of spotted arrays ( 11, 12, 14). The Pearson correlation coefficients between the matched frozen and formalin-fixed tumors ranged from 0.61 to 0.96 with average of 0.78. Because the samples were macrodissected for tumor cellularity >70%, the 78% concordance is within the expected experimental variation. The Pearson correlation coefficient from formalin-fixed samples repeated twice using same RNA isolation but two separate amplification was 0.96 (sample VU6018), whereas repeated experiments from two separate RNA isolations of the same tumor were 0.82 and 0.90 (samples MD147E5 and VU300663). These results are comparable with the Pearson correlation values obtained when studying RNAs obtained from fresh frozen tissues. Data from this study were deposited in the NIH Gene Expression Omnibus database under accession number GSE2837.

RNA expression analyses by real-time PCR. Microarray expression results for a subset of tumors were confirmed by separate real-time PCR analyses of KRT14, ACTN1, COL5A1, and PLEC1. Twenty RNA samples from formalin-fixed tumors were analyzed, and 17 samples gave reliable data on all three real-time PCR control genes, including B2M, UBC, and Eukaryotic 18S. We observed measurable expression of KRT14 in 15, ACTN1 in 10, COL5A1 in 6, and PLEC1 in 9 of 17 samples, albeit with high Ct values indicating the difficulty with which these signals were obtained. On average, KRT14 expression was increased 261-fold in high-risk tumors compared with low-risk tumors [95% confidence interval (95% CI), 12.1-5404.7; P = 0.0018] by real-time PCR, whereas 2.2-fold increase by DNA microarray (95% CI, 1.2-4.2; P = 0.018). ACTN1, COL5A1, and PLEC1 could not be statistically compared due to small numbers of samples that yielded reliable measurements. The low numbers of samples that had measurable levels of RNA by real-time PCR are probably due to the location and experimental condition of the commercially available probes, which are not optimized for RNA from formalin-fixed tumors. Affymetrix arrays use 11 different sequences from different locations on the gene (each 25 bases long) for each probe, and the expression value of the “probe set” or gene is estimated from hybridization to all 11 probes. In the Affymetrix X3P GeneChip, the probes are biased toward the 3′-end of the transcript. Therefore, the issue of probe location for short RNA fragments from formalin fixation is probably less critical, and the expression data are reliably measured, albeit less present calls compared with frozen tumors.

Molecular classification of HNSCC by intrinsic gene set. The intrinsic 950-gene set (1,101 Affymetrix probes) was generated from 8-paired tumor samples as described previously ( 2, 3) and used to classify the 40 samples on X3P GeneChips based on their gene expression levels. All six frozen and formalin-fixed tumor pairs and three formalin-fixed repeated samples were grouped immediately adjacent to each other in pairs at the terminal branches on hierarchical clustering ( Fig. 1 ; ref. 13). To test the validity of the gene set, molecular classification of previously published 60 frozen primary HNSCC tumors was repeated using 349 genes that were present on both cDNA microarray and X3P GeneChip probe sets using common UniGene identifiers ( 2). Despite the limitation of using only 349 genes (37% of the newly generated intrinsic gene set), the frozen tumors that previously mapped to the high-risk group were identified with >80% concordance and a κ statistic of 0.70, indicating good agreement between the two gene expression assessments. We repeated this analysis using the intrinsic gene set from the frozen tumor data set and applied this to the formalin-fixed tumor data set with similar results (data not shown). These findings indicate that the gene expression data generated from formalin-fixed tumors is informative, and the data yields comparable results against frozen tumors and vice versa.

• Download figure
• Open in new tab
• Download powerpoint

Figure 1.

Hierarchical clustering of 40 samples analyzed on X3P GeneChips, including 33 formalin-fixed tumors, 6 matching frozen tumors, and 1 formalin-fixed mucosal epithelium, using the 950 intrinsic gene set for molecular classification. The frozen samples are labeled as a suffix to a sample name with “-Frozen.” The repeated samples are labeled as a suffix to the same name with “-2” or “-3”. The tumors with multiple samples are labeled with blue bar.

Determination of genes predictive of the high-risk group. The intrinsic gene set analysis uses unsupervised statistical analysis. To determine whether more accurate classification could be obtained using supervised analysis, supervised principal components analysis was used to predict each sample individually using the recurrence event and time as a supervising variable. A 75-gene list was determined to be the most predictive of recurrence ( Fig. 2 ). Of the 75 predictive genes, 42 genes had known functions ( Table 2 ). Twenty-six of 42 genes had higher expression, and the remaining 16 genes had lower expression in the high-risk group. The same tumors that were used for training were classified based on the expression of these 75 genes and analyzed for correlation with RFS. Results were plotted as Kaplan-Meier plots and analyzed by log-rank test ( Fig. 3A ). Patients with tumors classified as high-risk had a statistically significantly worse RFS compared with patients assigned to the low-risk group (median follow-up, 33 months; median survival high-risk group 11 months, low-risk group 65 months; P = 0.002, log-rank test). The frozen tumor set was used as an independent validation set. Only 28 genes of 75 predictive genes overlapped between the two platforms. However, the difference of RFS between the high- and low-risk groups in the frozen tumors based on the 28 genes was still statistically significant (median follow-up, 28 months; median survival high-risk group 27 months, low-risk group 29 months; P = 0.033, log-rank test; Fig. 3B). Our results strongly suggest that there are significant biological differences between the low- and high-risk tumors based on the findings that only partially overlapping gene sets can distinguish the groups.

• Download figure
• Open in new tab
• Download powerpoint

Figure 2.

Hierarchical clustering of 28 tumors from formalin-fixed tissue with recurrence and survival data using the 75-gene list for the prediction of high-risk versus low-risk for poor RFS. The supervised principal components analysis was used to predict the risk for the individual samples. Red, high-risk tumors; blue, low-risk tumors.

• Download figure
• Open in new tab
• Download powerpoint

Figure 3.

Kaplan-Meier plot of RFS based on the risk prediction by the 75-gene predictive gene list for the 28 patients with formalin-fixed tumors (A) and for the 60 patients with frozen tumors (B).

Molecular characteristics of the high-risk tumors. The predictive genes for recurrence were involved in various cellular functions ( Table 2) making it difficult to determine any unifying pathways. Therefore, we applied GSEA to further interrogate molecular pathways or gene sets that distinguished high-risk from low-risk tumors. Eight gene sets were determined to be significantly different with a false discovery rate of <25%, which is the accepted statistical cutoff value of GSEA ( 7). Four of the eight gene sets that were enriched in the high-risk group were defined by genes involved in EMT ( 15), NF-κB signaling ( 16), and cellular adhesion (two gene sets procured from independent sources; refs. 7, 17). Four of the eight gene sets that were enriched in the low-risk group were defined by genes involved in mRNA splicing, ¹² cell cycle regulation (two gene sets procured from independent sources; refs. 7, 17), and mRNA processing. ¹² The enriched gene sets and their individual gene lists are also detailed in Supplementary Table S1.

To confirm the GSEA result suggesting the deregulation of NF-κB signaling in the high-risk group, we obtained the identities of 277 unique genes (305 clones on murine cDNA array) that were known to be modulated by NF-κB from a published study on squamous cell carcinoma generated from a keratinocyte cell line ( 18). UniGene mapping showed that 99 of the 277 unique genes were present on the Human Agilent microarray and designated as the NF-κB signature (Supplementary Table S2). Using this NF-κB signature, the 60 frozen HNSCC tumors were grouped using hierarchical clustering. As expected, most of the high-risk tumors clustered separately from other tumors again (P < 0.0001, Fisher's exact test), meaning that gene expression profiles indicating aberrant NF-κB signaling could distinguish high-risk from low-risk tumors ( Fig. 4 ).

• Download figure
• Open in new tab
• Download powerpoint

Figure 4.

A, hierarchical clustering of 60 frozen tumors from Chung et al. ( 2) using the 99-gene NF-κB signature known to be modulated by NF-κB (extracted from Loercher et al.; ref. 18). Left, 17 of the 21 samples in the high-risk group from the original study (red) with a distinct NF-κB signature, suggesting deregulated NF-κB signaling (P < 0.0001, Fisher's exact test). B, genes that had higher expression in the high-risk group among the 99-gene NF-κB signature.

Discussion

This report is the first to show that analyses of formalin-fixed tissue for global gene expression are feasible for class prediction analyses. This was accomplished with only 50 to 100 ng total RNA isolated from formalin-fixed tumors and yielded >8,000 data points per sample. The feasibility of obtaining DNA microarray data from archived formalin-fixed tissue specimens collected 13 years before analyses was also shown. Because formalin-fixed tissue is routinely stored in clinical pathology laboratories, gene expression array analyses of formalin-fixed tissue will allow correlative studies from large well-defined patient populations, including cohorts from completed clinical trials with multiyear clinical follow-up. In addition, because of institutional expertise, cost, or convenience, investigators use a wide variety of microarray platforms making comparison of results difficult. Our data suggest that variability introduced from the use of two different platforms (Affymetrix X3P GeneChip and Agilent Human cDNA microarray) as well as different sample preparation (formalin-fixed and frozen tumors) could be, at least in part, overcome by correction of the systematic bias using a well-established statistical methods, such as principal components analysis and distance-weighed discrimination ( 11, 12, 19), allowing the combination of the array data for statistical analyses. The ability to combine new and previously published data will address some of the issues with cost of array experiments.

More importantly, our results indicate that clinically and biologically meaningful information can be obtained as shown through reproducibility of the molecular classification and recurrence prediction between formalin-fixed and frozen tumors. By combining gene expression data from the formalin-fixed and frozen tumors, we have determined a 75-gene list that is highly predictive of recurrence. Among the genes that had higher expression in the high-risk of recurrence group, keratin 14, matrix metalloproteinase (MMP) 2, stratifin, and galectin 1 were shown previously to associate with poor prognosis. Keratin 14 was identified as a marker of high-risk patients and validated by immunohistochemistry in our previous study ( 2). MMPs are a family of endopeptidases involved in the breakdown of extracellular matrix (ECM), and MMP-2 degrades type IV collagen, which is a major structural component of basement membrane ( 20). In a study by Ruokolainen et al. ( 21), high expression of MMP-2 in HNSCC tumors by immunohistochemistry was associated with shorter survival as well as lymph node and distant metastases. Stratifin, also known as 14-3-3σ protein, is required for the transition of stem cells to transit-amplifying cells in keratinocytes, and the down-regulation of the protein leads to immortalization of primary human keratinocytes ( 22). Stratifin is thought to be an epithelial cell type–specific protein, and its role in p53-mediated G₂-M cell cycle arrest by inhibiting the activation of cyclin-dependent kinase 1 in colorectal cancer is well established ( 23, 24). Higher expression of stratifin in the high-risk tumors could be the result of p53 mutation or deregulation, which is present in 42% to 91% HNSCC ( 25– 27). Galectin 1 is one of the members of the β-galactoside-binding proteins and thought to modulate cell to matrix interaction and cell proliferation ( 28). Galectin 1 was shown recently to be a hypoxia-induced protein, and HNSCC tumors with positive immunohistochemistry staining in the ECM for galectin 1 had significantly worse overall survival ( 29).

The ultimate goal of the clinical use of gene expression data is to improve clinician's ability to stratify treatment based on tumor aggressiveness. Once high-risk patients are identified using the predictive 75-gene list, appropriate treatment will still have to be determined for this group of patients. The traditional approach has been treatment intensification by increasing drug dose, by combining treatment modalities, or by increasing treatment frequency. This approach has resulted in some improvement in disease outcome; however, it has also increased severe treatment-related toxicities and mortality. Understanding of differences in the underlying biology of the high-risk tumors may identify therapeutic targets that will have increased efficacy and decreased toxicities. When we tried to understand the biological contribution of each of the 75 genes to the disease recurrence, 42 genes had known molecular function, and the unifying pathways were difficult to delineate. Therefore, we applied GSEA to gain insight into the biological difference of the high-risk tumors using the entire 8,000 genes that passed the stringent filtering criteria in addition to the independent evaluation of each 42 genes. The most significantly enriched gene sets in the high-risk tumors were on pathways related to EMT, NF-κB activation, and cell adhesion deregulation, which may be a requirement of the EMT process.

The transformation of EMT phenotype has associated with changes in the morphology to spindle-shaped and motile fibroblastoid phenotype and in loss of tight- and adherens-junction proteins, which allows the tumor cells to pass through the basement membrane ( 15, 30, 31). EMT has been associated with late stage of tumor progression and metastasis ( 32). Presence of EMT in HNSCC has been seen in our previous molecular classification study as one of the poor prognosis indicators ( 2). Recently, EMT has gained significant attention clinically in non–small cell lung cancer due to the association with resistance to epidermal growth factor receptor tyrosine kinase inhibitors, such as erlotinib ( 33, 34). In addition, with the well-known involvement of Src kinase in EMT and the development of Src kinase inhibitors, such as AZ0530 and dasatinib, in clinical trials, the understanding of EMT has increased clinical significance and importance of EMT as a process that can be targeted by novel drugs. Activation of NF-κB is associated with risk for distant metastases and poor clinical outcome in different cancer types, including HNSCC ( 35– 39). Therefore, we interrogated our data using an independent NF-κB signature generated by Loercher at el. ( 18) in addition to the gene set procured in GSEA ( 7). Examples of the genes that were up-regulated within the NF-κB signature were MYC, PTEN, and HIF1-α, which are important regulators of growth and angiogenesis. It is known that activation of NF-κB by tumor necrosis factor can be augmented by activation of AKT ( 40). Therefore, increased PTEN expression observed in the high-risk tumors may represent a cellular compensatory response to AKT activation in these tumors.

Although our data suggest that formalin-fixed paraffin-embedded tumors can be used to generate informative gene expression data by DNA microarrays, we advocate continued collection of frozen tissue using standardized protocol in future prospective studies because of the versatility of frozen specimens for many molecular assays and because our data suggest that fresh frozen tumors yield better quality of data. However, the use of formalin-fixed tissues from archived samples for gene expression carries the tremendous advantage of increasing the number of available specimens annotated with complete clinical data. The ability to combine various data sets from multiple laboratories will make a significant effect in translational research. Our study also provides additional evidence for the power of gene expression analysis to distinguish patients who are at high risk of recurrence. This constitutes independent validation of the results obtained in our earlier study ( 2). Comparison of gene expression from high- and low-risk tumors revealed that high-risk tumors have deregulated genes involved in EMT as well as those involved in NF-κB signaling. These pathways should be targeted as novel treatment options, including currently available agents, such as Src kinase inhibitors, and NF-κB or proteasome inhibitors.