Tissue-Specific Gene Expression of Head and Neck Squamous Cell Carcinoma In Vivo by Complementary DNA Microarray Analysis

Objectives  To identify distinct gene expression profiles of human head and neck squamous cell carcinomas (HNSCCAs) using complementary DNA (cDNA) microarray analysis and to create a preliminary, comprehensive database of HNSCCA gene expression.

Patients and Methods  Nine patients with histologically confirmed HNSCCAs, staged according to the American Joint Committee on Cancer, were enrolled. The HNSCCA tumor tissue and normal mucosal tissue were harvested at the time of surgery. A cDNA library was constructed from the paired fresh-frozen human surgical specimens of HNSCCAs and nonmalignant epithelial tissues. Biotinylated RNA was transcribed from the cDNA library and hybridized to high-density microarrays containing approximately 12 000 human genes. Altered gene expression of HNSCCAs was identified by comparison to corresponding normal mucosal tissues after a bayesian statistical analysis of variance. Results were analyzed using the gene database of the National Institutes of Health. Hierarchical clustering of the genomic data sets was determined by similarity metrics based on Pearson correlation.

Results  Hierarchical clustering analysis revealed that the gene expression profiles obtained from the nonselected panel of 12 000 genes could distinguish the tumors from nonmalignant tissues. Gene expression changes were reproducibly observed in 227 genes representing previously identified chemokines, tumor suppressors, differentiation markers, matrix molecules, membrane receptors, and transcription factors that correlated with neoplasia, including 46 previously uncharacterized genes. Moreover, significant expression of the collagen type XI α1 gene and a novel gene was reproducibly observed in all 9 tumors, whereas these genes were virtually undetectable in their corresponding, adjacent nonmalignant tissues.

Conclusions  Complementary DNA microarray analysis of human HNSCCAs has produced a preliminary, comprehensive database of tumor-specific gene expression profiles and provided important insights into modeling gene expression changes implicated in carcinogenesis. A large-scale analysis of gene expression carries the future potential of identifying sensitive molecular markers for early tumor detection, prognosis, and novel targets for interceptive therapeutics.

AN OVERWHELMING body of evidence suggests that the phenotypic diversity of head and neck squamous cell carcinomas (HNSCCAs) is preceded or accompanied by corresponding genotypic changes. According to an established model by Kinzler and Vogelstein,^1,2 the development of neoplasia is a multistep process involving the accumulation of at least 7 genetic events for its progression into malignancy, primarily affecting various tumor suppressor genes and signaling tyrosine kinases.² Such modified genomes presumably lead to aberrations in the normal transcriptional program, resulting in the substantial deregulation of a battery of genes. Thus far, the Cancer Genome Anatomy Project index of tumor genes has already classified more than 40 000 genes that are directly or indirectly active in 1 or more cancers.^3,4

Nevertheless, systematic investigations into this complex molecular circuitry of carcinogenesis has been hampered by conventional techniques such as RNA Northern blot hybridization and ribonuclease protection assays that are both tedious and time-consuming and are typically performed by screening one gene at a time. Somewhat more sophisticated methods, such as differential display⁵ and Serial Analysis of Gene Expression,⁶ have been used to screen larger number of expressed genes. However, technical limitations render these techniques nonconducive to a comprehensive scale genomic survey.

Recently, technological advances have resulted in the emergence of an innovative technique based on the principle of complementary DNA (cDNA) library hybridization to high-density cartesian coordinate arrays containing ordered oligonucleotides of defined genes and expression sequence tags (ESTs), so-called microarrays or gene chips.⁷ Such an array platform permits simultaneous monitoring of thousands of gene expression events in a single experiment, providing the researcher with a new arsenal to analyze underlying pathomechanisms on a genome-wide scale. Although this approach has been applied to solid tumors of the head and neck, prior studies were limited by the paucity of genes loaded on the array (590-1190 genes),^8-10 had technical limitations due to the lower sensitivity of the nylon-based arrays,^9,10 or had statistical limitations due to sample size or inter-sample variation.^9-11 As such, comprehensive gene expression profiles that offer prognostic value or the identification of tumor marker genes that are exclusively expressed in HNSCCA tissues still remain elusive.

In this study, we demonstrate the application of the powerful 12 000-gene, photolithographic microarray based on gene transcripts hybridized from 9 human surgical specimens of HNSCCA and corresponding nonmalignant epithelial tissues. The resulting data revealed intriguing insights into the biological pathways of upper aerodigestive tract malignancy. Hierarchical clustering analysis revealed that the gene expression profiles obtained from the raw panel of 12 000 genes alone were sufficient to distinguish most tumors from nonmalignant tissues. Moreover, significant expression patterns were observed with 2 genes that are reproducibly expressed in all 9 tumors, whereas these 2 genes were virtually undetectable in normal epithelial tissue. Using novel statistical software, we developed a model for genes involved in tumorigenesis with the hope that this approach may eventually identify specific therapeutic avenues for prevention, prognosis, and treatment of HNSCCAs.

Surgical biopsy specimens were immediately stabilized by freezing the tissue in liquid nitrogen in the operating room. For RNA extraction, tissues were cut to approximately 1 cm³, suspended in TRIzol reagent (tissue:TRIzol = 1:15, wt/vol) RNA extraction buffer, then immediately subjected to homogenization by Polytron handheld homogenizer (model PT1200C; Brinkmann Instruments, Westbury, NY) at a maximum setting of 30 seconds. Total RNA was purified according to published TRIzol extraction protocols (QIAGEN, Valencia, Calif). After the ethanol precipitation step in the TRIzol extraction procedure, an additional cleanup step using an isolation kit (RN-easy Total RNA; QIAGEN) was incorporated for better yields of labeled complementary RNA (cRNA) during the in vitro transcription-labeling reaction. The RNA concentration was determined by spectroscopy, and the quality was confirmed by agarose gel electrophoresis in 3-(N-morpholino)propanesulfonic acid (MOPS) buffer, according to well-established protocols.

Double-stranded cDNA was synthesized from total RNA using a custom synthesis kit (SuperScript II Double-Stranded cDNA Synthesis Kit; Life Technologies, Carlsbad, Calif). Following ethanol extraction, in vitro transcription reactions were performed with a labeling kit (BioArray HighYield RNA Transcript Labeling Kit; Enzo, New York, NY) according to the manufacturer's protocol. Purified, labeled cRNA was quantified by spectrophotometric analysis, qualitatively analyzed for size distribution by gel electrophoresis, and fragmented to 30 to 60 base fragments with Tris-acetate (pH 8.1, 40mM), potassium acetate (100mM), and magnesium acetate (30mM) in a 20-µL volume heated for 35 minutes to 94°C.

U95A human microarrays (GeneChip; Affymetrix, Santa Clara, Calif) were prehybridized with hybridization buffer (1.0M sodium chloride, 100mM 4-morpholinoethane sulfonic acid (MES) [pH 6.6], 20mM EDTA, and 0.01% Tween 20) at 45°C for 10 minutes with rotation. Fragmented cRNA (15 µg for standard GeneChip arrays, 5 µg for test chip), sonicated herring sperm DNA (0.1 mg/mL final), GeneChip Eukaryotic Hybridization Controls (Affymetrix), and acetylated bovine serum albumin (BSA) (0.5 mg/mL) were then added to the hybridization buffer and incubated with the arrays for 16 hours at 45°C after which the probe array was washed (Fluidics station; Affymetrix) and amplified with staining solution (100mM MES, 0.1M sodium chloride, 0.01% Tween 20, and 0.005% antifoam) containing 10 µg/mL of streptavidin-phycoerythrin (Molecular Probes, Eugene, Ore). Acetylated BSA (2 µg/µL final) was applied for 10 minutes at 25°C. After 2 washes, the signal was amplified with a second staining solution (100mM MES, 0.1M sodium chloride, 0.01% Tween 20, and 0.005% antifoam) containing 3 µg/mL of anti–streptavidin biotinylated goat antibody (Vector Laboratories, Burlingame, Calif), normal goat IgG (0.1 mg/mL final), and acetylated BSA (2 µg/µL final) for 10 minutes at 25°C. Additional streptavidin-phycoerythrin staining at 25°C and a final wash were performed.

The resulting arrays were scanned twice in the GeneChip system confocal scanner (Hewlett-Packard GeneArray Scanner; Hewlett Packard Co, Palo Alto, Calif). To obtain average difference intensities, data analysis was first performed using computer software (GeneChip Expression Analysis Software version 3.3; Affymetrix). The Cyber-T software package (available from the Internet at http://www.igb.uci.edu/servers/cybert/) was then used to determine average intensity values of genes across 9 independent patient sets and to compare values for control and tumor specimens. This list was further restricted based on whether genes were consistently present or consistently absent in the 9 cases. Using a bayesian statistical analysis, significant changes were ranked by the assignment of P values.¹² Genes displaying intensities below 30 were empirically considered absent. Fold changes in gene expression were calculated by comparing average intensities in normal tissue with that of HNSCCA lesions. For genes that were entirely absent in one subclass of specimens but present in the others, relative fold changes were assigned by comparing intensities to a normalized value determined by the Affymetrix software. Hierarchical clustering of the genomic data sets was determined by similarity metrics based on the Pearson correlation equation.^13,14 Implementation of the clustering algorithms and dendrogram graphics was performed with the combined Cluster and TreeView¹⁵ software package under uncentered correlation and average linkage clustering settings.

The U95A human gene microarray is composed of oligonucleotide probe sets representing approximately 12 000 full-length, partially sequenced EST genes, including transcripts representing constitutively expressed housekeeping genes for internal control and normalization among samples. Of the sequences represented, approximately 5500 genes (45%) were detected in cRNA produced from snap frozen biopsy specimens of both HNSCCA and control epithelial tissues (Table 1), indicating abundant representation.¹¹ The amplified cRNA generated no nonspecific or unusual hybridization patterns and had no appreciable degradation with respect to its 5′ region. The transcripts that represent the internal control genes were detected in abundance, suggesting that the integrity of the tissue RNA was of high quality (data not shown).

In contrast to the clonal nature of tissue culture experiments, investigations of gross biopsy specimens in vivo are often confounded by factors such as patient-to-patient variation, the proximity of "normal" control tissue to disease sites, and overall heterogeneity of cell populations, all of which may affect gene expression profiles. To measure the extent of heterogeneity and physiologic variation among samples, pairwise comparisons of total gene expression levels of all 18 biopsy specimens were performed based on the Pearson correlation coefficient (R² = [n(Σxy) − (Σx)(Σy)]/{[n(Σx)² − (Σx)²][n(Σy)² − (Σy)²]}). Scatterplots of normal tissue vs normal tissue, tumor tissue vs tumor tissue, and normal tissue vs tumor tissue displayed overall correlation coefficients of 0.76, 0.68, and 0.64, respectively (Figure 1). Although these values are modest by clonal standards, the regression curves and R² values consistently indicate that gene expression variability is greatest between normal and tumor tissues, followed closely by the overall variability observed within the tumor group, whereas the normal tissues independently displayed the least variability among samples. This intriguing result suggests that the gene expression patterns between normal tissues and tumors are most divergent and that the magnitude of those differences may be sufficient to uncover tumor-related genes from "biological noise" and the inherent genetic variability among patient-to-patient samples.

A natural basis for organizing gene expression data is to group genes with similar patterns of expression and to decipher common distinguishing features of the cancer tissue. The average linkage hierarchical clustering algorithm was performed to group genes on the basis of similarity in the pattern with which their expression varied throughout all samples. The same unbiased approach was used to group the 18 experimental tissues (from 9 tumors and 9 controls) on the basis of overall similarity in their gene expression profiles. We focused first on the nonexcluded panel of the entire 12 000 gene set. The input data were based on the average difference intensity gene expression value as determined by Affymetrix software, and the relationships among samples were summarized in a dendrogram (Figure 2). The results reveal that most normal tissues (7/9) group together and, likewise, most tumor tissues (6/9) cosegregated into a separate cluster, suggesting different gene expression patterns between the 2 major groups. Strikingly, the 2 distinct tumors that clustered furthermost from the common tumor group represented the only 2 cases with a prior history of HNSCCA among this series (Table 2). Although previous microarray studies have reported gene-based tissue segregation that distinguishes various pathophysiologic or histopathologic phenotypes from control tissue,^16-21 its discriminatory power was based on a preselected pool of differentially expressed genes. Our study suggests that the analysis of an unbiased panel of 12 000 genes alone may be sufficient for the detection of the major subtypes of HNSCCA tissue.

A common approach to identifying significant gene expression changes is based on selection by fold change relative to the reference tissue. In our study, a pairwise analysis of all 9 tumors with corresponding control tissues revealed an average of 1790 genes that displayed a greater than 3-fold expression change (data not shown). However, many of these observations were changes with high interexperimental variance and likely represented artifacts. An alternative approach is based on applying the t test within a bayesian statistical framework and calculating P values rather than fold changes to rank significant expression changes. This method was preferred, since it more conservatively identifies genes with expression changes and has been shown to reduce false-positive results in microarray studies with limited number of replicates.^12,22 In short, the bayesian algorithm generates rank order lists of genes that are most reproducibly different among groups.

After implementation of the combined t test–bayesian statistical analysis, significant changes in messenger RNA levels were detected in 615 genes (P<.01). Further increasing the stringency criteria in the selection of genes to those with P<.001 generated a register of 227 genes (Table 3). Moreover, 30 genes were identified with P<.00001 (Table 3). The molecular taxonomy of the 227 gene data set revealed a substantially enriched ratio of genes encoding extracellular proteins (25.1%), followed by cytosolic genes (22.5%), membrane-associated genes (17.6%), nuclear genes (13.7%), and mitochondrial genes (1.3%). In addition, we identified more than 30 previously characterized oncogenes (Table 3) and 45 HNSCCA-associated sequences (19.8%) that represent novel genes or ESTs.

According to the holistic model of cancer by Hanahan and Weinberg,²³ a malignant cell has to essentially acquire 6 biological alterations to dictate pathogenesis, namely self-sufficiency in proliferative growth signals, insensitivity to growth inhibitory signals, evasion of apoptosis, limitless replicative potential, sustained angiogenesis, and the induction of invasion. In support of this model, our study shows that a considerable number of the 227 HNSCCA-related genes can directly or indirectly feed into any number of these pathways (Table 3). Two prominent features of our data were the high representation of extracellular genes, most of which encode matrix proteins (56%), and a global down-regulation of extracellular genes involved in inflammation. Although the pathophysiologic role of such extracellular genes in HNSCCA is largely unknown, it is intriguing to speculate that this cluster of genes may represent novel molecular signatures specific for squamous cancers.^24-28

To test the efficacy of the bayesian statistical framework, hierarchical clustering was performed based on 1239 differentially expressed genes gated at P<.05 (Figure 3). Rather than merely recapitulating the affirmative trends seen in Figure 2, the combined algorithm resulted in a near-perfect clustering of tumor (8/9) and normal (8/9) tissues. The reproducibility between the 2 experiments was clearly evident. In both experiments, the tumors farthest from the common tumor group represented cases with prior histories of HNSCCA. Altogether, these improved results under bayesian analysis suggest that a highly precise tumor classification system can be generated by genes selected from a pool of less than 25% of the original sequences sampled.

Studies demonstrate that gene detection by microarrays correlate with standard methods such as Northern blot hybridizations, reverse transcriptase–polymerase chain reaction, and real-time quantitative polymerase chain reaction.^11,29 These quantitative and qualitative aspects of microarrays can be exploited to screen for putative molecular markers of HNSCCA. An ideal candidate for such a marker would be a transcript that is both robust in its expression and specific among all HNSCCA and normal tissue. In our study, we identified 2 genes that fulfilled these criteria. In both genes, significant expression was reproducibly observed in all 9 tumors, whereas these genes were virtually undetectable in their corresponding, adjacent nonmalignant tissues (Table 4). The identity of the first gene was the human collagen type XI α1 gene (COL11A1), which encodes an essential component of the interstitial extracellular matrix. Consistent with our study, this gene was detected in 27 (79%) of 34 colorectal cancer tissues, whereas it was not present in normal colon tissue.³⁰ Its role in malignancy, however, still remains unclear. The second gene (accession No. D79998) was a cDNA clone identified in the human myeloblast cell line KG-1.³¹ Interestingly, this gene is unrelated to any previously reported genes. Its 3.6-kilobase sequence reveals only partial homology to human chromosome 17 and incomplete homology with the remainder of the current human genome database. Therefore, we named this putative gene HSCA-1 (head and neck squamous carcinoma–associated gene), since this gene is highly likely to represent a novel gene. In the future, we will design specific oligonucleotide primer sets for HSCA-1 and COL11A1 to validate these findings by reverse transcriptase–polymerase chain reaction on future specimens of HNSCCAs.

The use of microarray analysis in determining gene expression in cancer has been applied in the past for several different cancers, such as colon,¹⁷ breast,^18,20,21 B-cell lymphoma,¹⁹ and, more recently, lung³² and head and neck cancers.^33-35 Through the use of this technique, Villaret et al³³ identified 13 genes that are overexpressed in HNSCCA when compared with normal tissue, of which 4 of these genes were novel. Using cluster analysis, Belbin et al³⁴ found 337 genes that distinguish 2 subgroups of HNSCCAs. They also identified 4 EST segments that clustered to only one subgroup.

Through cDNA microarray analysis and using the combined t test–bayesian statistical analysis, we were able to identify 227 genes, including 46 EST segments that were either overexpressed or underexpressed in HNSCCA when compared with normal epithelial tissue of the same patient. In a similar study, Al Moustafa et al³⁵ compared HNSCCA samples with normal tissue from the same patient and identified 213 statistically significant genes, 91 overexpressed and 122 underexpressed. Two genes of note that were underexpressed in HNSCCA were Claudin-7 and Connexin 31.1. In our study, we also find Caludin-7 to be underexpressed. Yet, we identify COL11A1 and the novel HSCA-1 gene as being overexpressed solely in HNSCCAs.

The gene expression profile of tumors can be performed using either RNA extracted from whole tumor tissue or that from isolated tumor cells. Pure population of tumor cells and normal cells can be isolated using the technique of laser capture microdissection.^11,36 This strategy is essential to identifying unique tumor markers, which may serve as potential therapeutic targets of cancer cells.³⁷ However, because the clinical behavior of a particular tumor depends on the characteristics of both tumor cells and stromal tissue, such as density of angiogenesis and lymphocyte infiltration, gene expression profile of whole tumor tissue is required to elucidate prognostic markers and biologic variations.

The ability to obtain quantitative information from multiple transcriptional profiles represents an exceptionally powerful means to explore cancer biology and to generate gene databases relevant to HNSCCA. Despite the challenges of heterogeneity in patient tissues, genes expected to play a role in tumorigenesis were selected, organized by their predicted protein localization, and clustered based on similarity. Such a model represents an exciting venture for further investigation, yet it is critical to be aware of the limitations imposed by microarray studies. The observed gene expression may reflect either increased transcription or decreased mRNA stability by the reference tissue. Certain transcripts may not necessarily encode bioactive proteins. Moreover, microarray experiments are vulnerable to problems of statistical inferences because many genes will show significant changes in gene expression purely by chance alone. Therefore, to generate biologically meaningful data, it was necessary to use additional statistical algorithms such as the bayesian analysis.¹²

These caveats notwithstanding, DNA microanalysis was an effective and efficient screening method to identify global gene expression profile. Several previously known genes associated with cancer were identified, supporting our statistical and experimental approach. In addition, 2 genes exclusively expressed in HNSCCA were identified: an isoform of the collagen gene observed in colon cancers and a putative novel gene, each of these molecules representing potential targets for further investigation. Taken together, quantitative analysis of gene expression by cDNA microarray analysis of human HNSCCA has produced a preliminary, comprehensive database of tumor-specific gene expression profiles and provides important insights into modeling gene expression changes implicated in carcinogenesis. A large-scale analysis of gene expression carries the future potential of identifying sensitive molecular markers for early tumor detection and novel targets for interceptive therapeutics.

Corresponding author and reprints: Fang-An Chen, MD, PhD, Department of Otolaryngology, New York University School of Medicine, Suite 7U, Skirball Building, 530 First Ave, New York, NY 10016 (e-mail: fangan.chen@med.nyu.edu).

Accepted for publication January 15, 2003.

This study was supported in part by the George E. Hall Research Fund. Dr Sok was partially supported by a Master Student Training Program grant from the National Institutes of Health, Bethesda, Md.

This study was presented at the annual meeting of the American Head and Neck Society Meeting, Boca Raton, Fla, May 11-13, 2002.

We gratefully acknowledge the assistance of J. Denis Heck, PhD, and Kim Nguyen, PhD, with the microarrays.