Gene Expression Differences Associated with Human Papillomavirus Status in Head and Neck Squamous Cell Carcinoma

Materials and Methods

Patient selection and specimen collection. Thirty-six freshly frozen tumor samples were prospectively collected from patients undergoing surgery or biopsy for HNSCC at the University of North Carolina at Chapel Hill (21 patients) and Vanderbilt University (15 patients; see Table 1 and Supplementary Data). All tissues were snap-frozen in liquid nitrogen within 30 minutes of surgical resection or biopsy and kept at −80°C until further processing. All patients consented to participation in this study under protocols approved by the Institutional Review Boards at the two institutions. A previous gene expression study (12) included the 21 tumors from University of North Carolina reported here, but to allow comparison with the specimens from Vanderbilt and because the Agilent platform used in the prior study was discontinued, a completely new expression analysis was done using the Affymetrix platform (see below).

HPV detection and DNA sequencing. Tumor DNAs were tested for the presence of HPV DNA using a previously established PCR-based method (13). This method employs degenerate PCR primers (MY09 and MY11, WD72/76 and WD66/67/154) that are designed to represent highly conserved HPV L1 and E6 sequences present in all major types of HPV. In addition, all HPV-positive samples were also tested with a HPV16-specific PCR for E7 (primer A, 5′-GGACCGGTCGATGTATGTCT-3′ and primer B, 3′-TAAAACCATCCATTACATCCCG-5′). As a positive control for amplification, primers for β-globin are included with each sample (13). Optimal conditions for this combined PCR were determined using DNA from the cervical carcinoma cell line SiHa, which harbors on average two copies of HPV16 DNA per cell (13). Other positive control cell lines were CaSki (HPV16) and HeLa (HPV18). For each case, 200 ng of tumor DNA were tested for the presence of HPV DNA. PCR samples, which showed amplification products indicating the presence of HPV were purified using PCR purification columns (Qiagen, Valencia, CA) and subjected to bidirectional sequence analysis. In all of such cases, a positive identification of the HPV type could be made.

RNA isolation and DNA microarray analysis. Each tumor was examined by H&E staining to ensure presence of tumor and enriched by macrodissection to achieve a minimum of 70% tumor cells in each preparation. Total RNA was purified from frozen tumors using Qiagen RNeasy Mini kit according to the manufacturer's recommendations (Qiagen) using ∼10 to 20 mg of wet tissue from each sample. Fifty nanograms of the total RNA were amplified using NuGEN Ovation Biotin RNA Amplification and Labeling kit (NuGen, San Carlos, CA) according to the manufacturer's recommendations. The NuGEN Ovation amplification methodology uses an isothermal linear amplification using random hexamers. This technology provides sensitive and rapid whole-genome amplification without introducing a significant bias toward the 3′ end of the transcripts (14). Fifteen micrograms of biotin-labeled aRNA were fragmented, and the quality of the RNA was reconfirmed using the Agilent RNA 6000 Nano LabChip kit and Agilent 2100 bioanalyzer. The fragmented, biotin-labeled aRNA was combined with the hybridization mix and loaded on to the Affymetrix Human Genome U133 plus 2.0 GeneChip. After hybridization, the GeneChip was washed, stained with streptavidin/phycoerythrin conjugate and biotinylated antibody, and scanned according to the manufacturer's recommendations. The raw microarray data was normalized using Perfect Match software for further statistical analyses.

Gene expression data analysis. The genes that were differentially expressed in HPV⁺ and HPV⁻ tumors were selected based on Significance Analysis of Microarrays (15). The selected genes were verified for statistically significant prediction power using the class prediction model based upon the compound covariate method (16, 17). This class prediction model determined whether the patterns of gene expression can identify two classes of HPV⁺ versus HPV⁻ tumors. The accuracy of the classification rate using the selected genes was estimated using the leave-one-out cross-validation. The pattern among the statistically significant discriminator genes was investigated using hierarchical clustering algorithm (18, 19).

Chromosome mapping of expression data. DIGMAP was done as described before (11). Briefly, chromosomal locus information for Affymetrix probes was retrieved in a batch mode from our local Gene Annotation Project database. Genes exhibiting significant differential expression were then identified by T test (significance threshold P < 0.01) using the MEV software package (TMEV2.1; ref. 20). The T scores were log-transformed reciprocal Ps [log₁₀(P⁻¹)]. The output files from these statistical analyses were processed by the DIGMAP Viewer and differential flag regions mapping programs. We implemented a sliding window method using perl scripts to compute total T scores per million base pairs (Mbp) from neighboring genes. In this study, we determined that a window size of 10 genes was optimal for visualizing differential flag regions without loss of sensitivity and low noise levels. The sliding windows overlap each other by one gene locus (i.e., step size = 1) to cover the entire chromosome, and the gene expression profiles are displayed as a moving average per Mbp. Normalized T scores for a window size 10 were calculated by summing 10 T scores (T) from within the window, then divided by the window length in actual genomic distance. Criteria used for identifying differential flag regions were a normalized T score >2 SD from the mean of total normalized T score (in this case, the cutoff value is 4.9) for all sliding windows.

Confirmation of expression data by real-time PCR. Total RNA from seven of the eight HPV⁺ tumors was available for real-time PCR (RT-PCR) analysis, whereas an equal number of seven RNAs were chosen from the HPV⁻ tumors for comparisons of expression levels. Fifty nanograms of total RNA were amplified using the NuGEN WT-Ovation RNA Amplification kit (NuGEN; ref. 14). The amplified cDNA was cleaned using the Qiagen PCR purification kit (Qiagen). Five genes among the 91 statistically significant genes were analyzed by RT-PCR: TAF7L, CDKN2A, SYCP2, RFC4, and NAP1L2 using Applied Biosystems Taqman FAM labeled probes (Applied Biosystems, Foster City, CA). An additional RT-PCR assay was done to test for HPV16-E6 expression in seven of the eight HPV⁺ tumors. The endogenous genes 18S, PPIA, and GUSB were used as internal calibration standards. The average of these three internal genes was used to normalize the RT-PCR results from the set of five genes and HPV16-E6. Twenty-five nanograms of amplified cDNA were used per reaction, and the probes were obtained from Applied Biosystems. Analysis of each sample was done in quadruplicate on an Applied Biosystems 7900HT instrument (Applied Biosystems).

Statistical analyses. Descriptive statistics were generated and tested with Fisher's Exact and Wilcoxon rank sum tests using the SAS/STAT statistical analysis package (SAS Institute, Research Triangle Park, NC). RT-PCR data were analyzed by the 2^−ΔΔCT method as described previously (21). Briefly, the average C_t was calculated for the four replicate analyses of the three control genes (18S, PPIA, and GUSB), and this value was subtracted from the average C_t calculated from the four replicate analyses for the genes of interest. Expression differences were compared using these normalized ΔC_t values between the HPV⁺ and HPV⁻ tumors, and the observed differences were tested using Student's t test. Two-tailed Ps < 0.05 were considered statistically significant.

Results

HPV detection in HNSCC tumors. A total of 36 DNA samples obtained from HNSCC specimens, representing all subclasses except nasopharynx (Table 1), were subjected to HPV analysis using PCR amplification of E6, E7, and L1 sequences. An example of an analysis of HPV E6 is shown in Fig. 1. Seven tumors were positive for both E6 and L1 PCR analyses, whereas one tumor was positive for L1 and not for either E6 or E7 (970108), and one tumor was positive for E6 and E7 but not for L1 (300171). Based on these results, the eight tumors that were positive for E6 and E7 were classified as HPV⁺, whereas the one specimen that was only positive for L1 was classified as HPV⁻. DNA sequence analysis of the E6 PCR products revealed that all tumors harbored type 16 HPV. RNA expression of HPV16-E6 was confirmed by RT-PCR in seven of the eight HPV⁺ tumors that were available for this analysis, including the L1-negative sample (300171), which also showed E6 expression. Patients with HPV⁺ tumors were on average younger than those with HPV⁻ tumors (median age, 49 versus 58 years), although this difference did not reach statistical significance (P = 0.08, Wilcoxon rank sum test; Table 1). HPV was significantly overrepresented in tumors originating from the oropharynx (seven of eight HPV⁺ tumors), whereas HPV was also observed in one of the nine tumors originating from the larynx. None of the 15 oral cavity tumors analyzed harbored HPV. No significant differences were observed with respect to race, tobacco use, alcohol use, clinical and pathologic stage, or tumor differentiation.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 1.

HPV analysis in HNSCC samples by amplification of E6 and β-globin. Positive control cell lines (lanes 1-3); negative controls, no DNA added (lanes 22 and 23). SiHa (lane 1), HeLa (lane 2), CaSki (lane 3), size standard (lanes 4 and 21), DNA samples derived from HNSCC tissue (lanes 5-20). Patient samples (lanes 5, 8, 12, and 18) show the shorter amplification product of HPV E6. The β-globin amplification product was visualized in all samples, indicating the integrity of the starting DNA.

Gene expression differences with HPV status in HNSCC tumors. To identify the genes that were differentially expressed between the eight HPV⁺ and 28 HPV⁻ tumors, statistical analysis using HPV status as the supervising variable was done (15). Among the ∼25,000 genes on the DNA microarray, 91 differentially expressed genes were highly statistically significant with a false-discovery rate (FDR) of <0.3% in classifying HPV⁺ versus HPV⁻ tumors (Fig. 2; Supplementary Data). Thus, <1 of the 91 genes (<1.1%) in this list is expected to make the cutoff to be classified as differentially expressed due to a chance event. The accuracy of the 91-gene set to predict HPV status by leave-one-out cross-validation was 100%. Of the 91 genes, 89 genes had higher average expression levels in the HPV⁺ tumors, whereas only two genes were on average expressed at lower levels in the HPV⁺ tumors (Table 2; Supplementary Data). One of the most significant differentially expressed genes was cyclin-dependent kinase inhibitor 2A (CDKN2A), which encodes the p16^INK4A tumor suppressor protein. Other genes with higher expression in the HPV⁺ group were cell cycle regulators other than p16^INK4A (p18 and CDC7), transcription factors (TAF7L, RFC4, RPA2, and TFDP2), the cell adhesion molecule (TCAM1), and several sequences defined by expressed sequence tags only (Table 2; Supplementary Data). The genes that had lower expression in the HPV⁺ group were NAP1L2 and KIRREL. Data from this study were deposited in the NIH Gene Expression Omnibus database under accession no. GSE3292.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 2.

Cluster diagram of 91 genes that are differentially expressed between HPV⁺ and HPV⁻ HNSCC tumors. HPV⁺ tumors form a separate cluster (right).

Gene expression cluster on chromosome 3q. HNSCCs harbor a large number of cytogenetic changes, but it is not known if any of these are associated with HPV status of the tumor and differential expression of the genes. To identify potential cytogenetic gains or losses using gene expression differences, we aligned the known chromosomal location of all genes to their expression levels using DIGMAP (11). Clusters of genes across the full genome were then tested for expression differences between the HPV⁺ and HPV⁻ groups. A graphical representation of this analysis is shown in Fig. 3. The most significant cluster of genes with increased expression in HPV⁺ tumors was observed on the long arm of chromosome 3. The region ranged from 3q24-ter, with several subclusters of genes with increased expression levels. We then ranked all genes in this cluster according to their expression difference between HPV⁺ and HPV⁻ tumors (Table 3).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 3.

Analysis of gene expression in HPV-associated HNSCC across the human genome. A, normalized T score, indicating areas of significantly higher number of genes with high levels of expression in HPV⁺ versus HPV⁻ tumors. B, detailed plot of chromosome 3, indicating the areas on the long arm (3q24-ter) with highly significant scores.

RNA expression analysis by RT-PCR. Microarray expression results for a subset of differentially expressed genes and HPV DNA results were confirmed by separate RT-PCR analyses. Seven of the eight HPV⁺ tumors were available for RT-PCR analysis, and all tested positive for E6 expression as measured by RT-PCR, whereas none of seven selected HPV⁻ tumors expressed detectable levels of HPV16-E6 (see Supplementary Fig. 1). Furthermore, expression levels of CDKN2, CDKN2A, TAF7L, SYCP2, and RFC4, which were increased in HPV⁺ tumors with FDR < 0.3%, as shown by microarray analyses, were also significantly increased in the RT-PCR analysis (see Supplementary Fig. 1). On average, TAF7L expression was increased 220-fold in HPV⁺ tumors [P = 0.0001; 95% confidence interval (95% CI), 28-1,744], SYCP2 increased 56-fold (P < 0.0001; 95% CI, 16-199), RFC4 increased 3.4-fold (P = 0.019; 95% CI, 1.3-9), and CDKN2 (p16^INK4A) increased 35-fold (P = 0.006; 95% CI, 3-358). Of these genes, RFC4 is located on the chromosome 3 region identified by DIGMAP (3q27.3). Decreased expression was found for NAP1L2, although the results did not reach statistical significance (7-fold decreased, P = 0.285; 95% CI, 0.2-286).

Discussion

New data on the effects of HPV infection as a possible etiologic agent in HNSCC has emerged in recent years, indicating that HPV⁺ tumors may represent a distinct subset that can be identified by specific molecular features (8, 9). However, a comprehensive analysis of the gene expression differences between HPV⁺ and HPV⁻ HNSCCs has not been reported. In this article, we determined global gene expression in eight HPV⁺ and 28 HPV⁻ tumors and compared the two groups to identify differences in gene expression patterns.

Patient characteristics. In terms of demographics, our study agreed with previously published results showing that HPV is mostly associated with tumors from the oropharynx, and that patients with HPV⁺ tumors are generally younger than those with HPV⁻ tumors (7, 22). In our patient population, seven of the eight HPV⁺ tumors originated in the oropharynx, and patients with HPV⁺ tumors were significantly younger than those with HPV⁻ tumors. We were unable to determine any (negative) association with smoking because of the relatively small size of our sample and the fact that all patients had some history of tobacco exposure. From these observations, we conclude that our patient cohort is representative of the general HNSCC patient population.

Gene expression analysis. The supervised analyses of the gene expression data showed 91 genes that were differentially expressed in HPV⁺ and HPV⁻ tumors. Because all but one of the HPV⁺ HNSCCs analyzed arose in the oropharynx, it is impossible to rule out confounding effects due to tumor site. However, in our previous molecular classification study of 60 HNSCCs, we did not find a specific gene expression pattern based on the anatomic subsite (12). Further studies examining larger numbers of HPV⁺ HNSCCs from multiple subsites will be needed to clarify this issue. However, it is also reassuring that the current gene list based on difference in HPV status identifies previously identified HPV-specific features in uterine cervical carcinomas, such as increased expression of CDKN2A (p16^INK4A), RFC4, MCM2, LIG1, and TFDP2, and decreased expression of NAP1L2 (23). We also found overrepresentation of highly expressed genes on chromosome 3q24-qter. Increased copy numbers of this chromosomal region were previously reported in association with positive HPV status (5). For these reasons, we believe that the observed gene expression differences are based on the presence or absence of HPV rather than on any hypothetical differences between oropharyngeal carcinomas and other upper airway squamous cell carcinomas.

Among the most prominent genes up-regulated in HPV⁺ tumors, several are involved in transcription and cell cycle regulation. One such cell cycle inhibitor is CDKN2A, which encodes the p16^INK4A tumor suppressor protein. p16^INK4A is a cyclin-dependent kinase inhibitor in the Rb tumor suppressor pathway. Increased expression of p16^INK4A may potentially reflect loss of a negative feedback loop associated with inactivation of the retinoblastoma tumor suppressor protein (pRB) by HPV E7. Increased levels of p16^INK4A are strongly correlated with HPV infection in cervical and in head and neck carcinomas (8, 24). Overexpression of p16^INK4A is highly correlated with the presence of HPV and has been used as a surrogate marker for HPV (24).

One of the most significantly up-regulated genes was TAF7L, a TATA box binding protein–associated factor, an X-linked gene that is normally only expressed in testis and spermatogonia (25). TAF7L plays a role in regulation of the transcription factor IID during male germ cell differentiation, where it replaces TAF7 in the transcription factor IID protein complex (26). This complex is required for transcription by promoters targeted by RNA polymerase II. TAF7 has been implicated in activator protein transcription regulation as a possible cofactor that binds to c-Jun (27). Expression of the 32-kDa subunit of replication protein A (RPA), encoded by RPA2, was also increased in HPV⁺ tumors. RPA is a tetrameric protein complex that binds ssDNA and facilitates DNA damage sensing (28). The HPV E1 protein, an ATP-dependent viral DNA helicase, binds and recruits RPA to sites of DNA replication, through binding of the 70-kDa subunit of RPA (RPA1; ref. 29).

In the microarray analysis, there were only two genes significantly down-regulated in HPV⁺ tumors: NAP1L2 and KIRREL (Table 2). However, RT-PCR analysis of NAP1L2 did not reach statistical significance, presumably because only two of the seven HPV⁻ tumors showed high NAP1L2 expression levels compared with the HPV⁺ tumors (see Supplementary Fig. 1). NAP1L2 encodes a member of the nucleosome assembly protein (NAP) family with high expression levels in brain and a possible role in neuronal cell cycle regulation (30). KIRREL (NEPH1) is a nephrin-related member of the immunoglobulin superfamily involved in cell-cell interactions related to selective ultrafiltration by podocytes in the kidney (31). However, expression is not limited to the kidney, and KIRREL mRNA is also found in the colon, pancreas, heart, and spleen (32). It is currently unclear what role a down-regulation of expression of these two genes would have in HPV-mediated carcinogenesis.

DIGMAP. Assignment of expression analyses with known chromosomal location of genes using DIGMAP revealed a cluster on chromosome 3q24-qter with a large number of genes with higher levels of expression in HPV⁺ versus HPV⁻ tumors. Interestingly, a prior study using comparative genome hybridization on 25 primary tonsillar carcinomas suggested of a gain on chromosome 3q24-ter in HPV⁺ tumors compared with HPV⁻ tumors (5). The observed expression difference in our series is highly significant, which indicates that either the expression levels are greatly increased in the tumors with chromosome gains on 3q, or that gene transcription is elevated in the tumors with unaltered DNA copy numbers. The 3q region identified by DIGMAP includes a known amplicon at 3q26 that has been described in HNSCC and lung carcinomas (33–35). Candidate oncogenes within this amplicon are PIK3CA (35, 36), SCRO (37), the RNA component of telomerase (TERC; ref. 38), ZNF639 (39), and p63 (40). Of these genes, SCRO and TERC are not represented on the Affymetrix Human 133U Plus 2.0 GeneChip. Amplification and overexpression of PIK3CA has been reported for HNSCC (35) and cervical carcinomas (41), although PIK3CA amplification has not been associated with HPV infection (42) and did not reach statistical significance in the present study.

The other differentially expressed genes within chromosomal band 3q26 (MYNN, NDUFB5, and ECT2) have not previously been associated with amplification of this region. ECT2 belongs to the dbl family of oncogenes that encode large, growth-regulatory molecules with guanine nucleotide exchange factor activity of towards Rho-family GTPases (43). Deregulation and mislocation of ECT2 has been linked to malignant transformation in mouse 3T3 cells. Other genes in the 3q24-ter region identified by DIGMAP do not localize to the 3q26 amplicon, which implies that the difference in gene expression profiles between HPV⁺ and HPV⁻ tumors can not be fully explained by just amplification at 3q26. Furthermore, the differentially expressed genes in this chromosomal region were mainly involved in DNA replication and transcription, such as replication factor C4 (RFC4), TFDP2, POLR2H, SHOX2, and SFRS10. RFC4 encodes the 37-kDa subunit that is part of the heteropentamer RFC (also called activator 1), an accessory protein required for the elongation of primed DNA templates by DNA polymerase δ and ε (44). Proliferating cell nuclear antigen is also part of the elongation complex. Proliferating cell nuclear antigen is a target of retinoblastoma signaling, but RFC does not respond to this pathway and is not responsive to E2F signaling (45).

Through global gene expression analyses of HNSCC tumors, we have identified differences in gene expression profile of HPV⁺ and HPV⁻ tumors and gained insight into the biological effect of HPV infection in HNSCC. The most prominent genes up-regulated in HPV⁺ tumors were involved in transcription and cell cycle regulation. These genes are interesting from a clinical targeting standpoint because some have already been identified for potential therapeutic benefit and novel agents have been developed that are currently in clinical trials. Examples are transcription inhibitors, such as histone deacetylase inhibitors (SAHA, MS-275, etc.; ref. 46), and cell cycle inhibitors, such as CDK and Chk kinase inhibitors (flavopiridol, UCN-01, 17-AAG, etc.; ref. 47). Combined with previous analyses of chromosomal imbalance, analyses such as these may help to further identify potential therapeutic targets in altered chromosomal regions. Additionally, gene expression changes driven by HPV status independent of chromosomal abnormalities may also identify cellular mediators of HPV oncogenesis.