Characterization of Epigenomic Alterations in HPV16⁺ Head and Neck Squamous Cell Carcinomas

Head and neck cancer is the sixth most common cancer worldwide and accounts for 1% to 2% of cancer-related deaths (1). Over 90% of head and neck cancers are squamous cell carcinomas (HNSCC) and typically arise from the oral cavity, nasopharynx, oropharynx, hypopharynx, and larynx. The etiology of HNSCC includes tobacco exposure, alcohol consumption, and human papillomavirus (HPV) infection. The incidence of HPV-associated HNSCC has significantly increased over the last two decades (2). In addition to HNSCC, high-risk HPV infection, predominately HPV16 and HPV18, has been causally associated with many different cancer types including those of the cervix, penis, anus, vulva, and vagina. HPV16 is the most commonly implicated HPV type in HNSCC (3) and frequently involves the oropharynx (3, 4). Of clinical relevance, HPV⁺ oropharyngeal squamous cell cancers (OPSCC) expressing wild-type p53 and high levels of p16 have better prognosis and response to chemoradiation (2, 5, 6).

The pathophysiology of HPV-associated HNSCC involves integration of HPV virus into the host genome, upregulation of oncoproteins E6 and E7 (7), and degradation of tumor suppressors p53 and Rb, respectively (8). Additional molecular alterations occur at the epigenetic level and include DNA histone modification, chromatin remodeling, modulation of noncoding RNAs and, most notably, aberrant DNA methylation. It has been postulated that these methylation changes may have important contributions to neoplastic progression, treatment response phenotypes and prognosis (5, 9). Several studies have demonstrated that HPV⁺ HNSCCs and cervical cancers have significant overexpression of DNA methyltransferases and higher levels of CpG island hypermethylation (9–11). Indeed, HPV-positive OPSCCs have been shown to express distinct gene methylation profiles as compared with HPV-negative OPSCCs (6, 12).

HPV-associated HNSCC (and in particular, oropharyngeal SCC) demonstrates distinct clinical behavior and prognosis, with improved survival and lower risk of recurrence, as compared with their HPV-negative counterparts (13). Indeed, as a whole, HPV⁺ HNSCC demonstrate increased response to radiation-based therapy, and consequently, there has been substantial interest in treatment de-escalation strategies to optimize the therapeutic ratio for patients affected with this disease (14). Nonetheless, it is notable that approximately 20% of HPV-positive HNSCC have poor treatment response and prognosis, and the biological basis for this disparate behavior remains to be fully elucidated (15). Beyond treatment, it is important to note that despite the predicted continued increase in incidence particularly of HPV⁺ HNSCC, there are no established or standardized screening methods for this disease (16).

Identifying the spectrum of epigenetic changes across the genome may provide clues as to the molecular mechanisms that drive the degree of treatment response and, in turn, prognosis. Critical methylation events may contribute to the differences seen between strong and poor treatment responders in both patients with HPV⁺ and HPV⁻ HNSCC. This line of investigation may yield actionable targets or pathways that could be manipulated to enhance treatment response and prognosis. Furthermore, the ability to detect DNA methylation events in patient saliva or plasma may allow the leveraging of epigenetic biomarkers for the purposes of screening, early diagnosis, and surveillance.

Although several studies have examined the genome-wide methylation differences between HPV⁺ and HPV⁻ HNSCC (17–24), comprehensive analyses that include comparisons with normal head and neck tissues, which allow insight into carcinogenesis-associated epigenetic changes, are lacking. Moreover, the HPV type–specific impact on methylation patterns has not been fully evaluated. The objective of this study is to leverage genome-wide methylation array data from The Cancer Genome Atlas (TCGA) to characterize differential methylation alterations associated with both HPV16⁺, other non-HPV16 HPV⁺, and completely nonviral HNSCC in comparison with normal head and neck tissues.

Read counts available in the Supplementary Material from Cao and colleagues (25) and Tang and colleagues (26) were downloaded and merged. TCGA samples with more than 100 reads from a specific virus type were considered positive for that virus type. Samples with no reads from any virus type were labeled as nonviral tumor samples, whereas tumor samples with 1 to 100 reads were categorized as “low counts” and not used in our analysis.

The raw Illumina450 IDAT files selected for HNSCCs were downloaded from TCGA and normalized using Control Probes with background correction as implemented in the minfi bioconductor package. Significant differentially methylated regions (DMR) were selected using the following criteria: FDR-corrected P value < 0.01, mean beta difference > 0.3 for HPV16 versus nonviral tumors. The mean beta difference threshold was increased to 0.4 for the normal versus HPV16 and normal versus nonviral comparisons, to account for the expected larger number of differences when comparing normal versus tumor tissue (27). In addition, probes with no significant neighbors, defined as being within 500 base pairs (28), were not considered to be significant.

All TCGA data were downloaded and processed as we have previous described (29). In short, the normalized and debatched RNA sequencing (RNA-seq) data file (EBPlusPlusAdjustPANCAN_IlluminaHiSeq_RNASeqV2.geneExp.tsv) were downloaded from Genomic Data Commons (https://gdc.cancer.gov/about-data/publications/pancanatlas). Samples listed in the Merged Sample Quality Annotation were removed due to potential quality issues. The expression values were log₂ transformed using log₂(x+1). For characterization of the tumor immune microenvironment, we used the same genes as previously reported (30).

The PanCancer Atlas DNA copy-number data (GISTIC data) were extracted from the Copy Number file (broad.mit.edu_PANCAN_Genome_Wide_SNP_6_whitelisted.seg http://api.gdc.cancer.gov/data/00a32f7a-c85f-4f86-850d-be53973cbc4d). Samples listed in the Merged Sample Quality Annotations spreadsheet (merged_ sample_quality_annotations.tsv) were excluded from further analysis. Based on GISTIC2.0 (31) DNA copy-number analysis, the following numeric values were assigned with the cBioPortal interpretation presented in parentheses: -2 = Deep Deletion (possibly a homozygous deletion); -1 = Shallow Deletion (possible heterozygous deletion); 0 = Diploid; +1 = Gain (a few additional copies, often broad); +2 = Amplification (more copies, often focal).

The raw Illumina450 IDAT files were downloaded from GEO (GSE68379; ref. 32) and processed in the same fashion as the TCGA data. Cell lines with missing TCGA “Tumor Type” were removed from the analyses.

The GSE38266 (33–35) dataset was downloaded from GEO and the processed data were used for validation analyses.

For principal component analysis (PCA) modeling, tumor samples with known viral type (n = 320) and normal samples (n = 50) were included. Of note, CpG probes with more than 20% missing values were removed. Each column was centered but not scaled before the PCA model was calculated.

The initial objective was for the selection of specific probes that were uniquely altered when comparing nonviral tumors with HPV16 tumors but not when comparing normal tissue with all tumor tissues. The converse goal was then to identify probes that differed in methylation between normal and all tumor tissues but were not altered in association with viral status. The following criteria were used to select the optimal probes for the partial least squares (PLS) model: Probes with greater than 20% missing values were removed. OF note, the TCGA methylation data were generated using the no longer available Illumina 450K array. To ensure future comparability of our results to the newer Illumina EPIC Chip, we removed probes that are not on this latter platform. HPV16-specific probes: For the HPV16⁺ PLS model, any probe that was differentially methylated in normal versus tumor with Δ β-value > 0.1 or P < 0.05 was removed, thus ensuring that only HPV16⁺-specific probes were selected. Subsequently, each probe was ranked using the geometric mean of the Δ β-value and the −log10(P value) for HPV16⁺ versus nonviral tumors. The top 100 significant probes (P < 0.0001 and Δ β-value > 0.2) with increased methylation in HPV16⁺ tumors compared with nonviral tumors and the top 100 probes with decreased methylation in HPV16⁺ tumors compared with nonviral tumors were selected for the HPV16⁺ PLS model. Tumor- versus normal-specific probes: For the tumor versus normal model, probes that were differentially methylated in HPV16 versus nonviral tumors with Δ β-value > 0.1 or P < 0.05 were removed, thus ensuring that only tumor versus normal-specific probes were selected. Following this selection, each probe was ranked using the geometric mean of the Δ β-value and the −log10(P value) for tumor versus normal. The top 100 significant probes (P < 0.0001 and Δ β-value > 0.2) with increased methylation in tumor compared with normal samples and the top 100 probes with decreased methylation in tumors compared with normal samples were selected for the final PLS model. A binary response variable was used with the designation of “1” for HPV16⁺ tumors and “0” for nonviral tumors in the HPV16 PLS model, whereas “1” was used for tumors and “0” for normal in the tumor versus normal PLS model. Samples with a predicted value above 0.5 for the HPV16 PLS model were considered to be HPV16-related and those less than or equal to 0.5 considered to be HPV16-unrelated. Cross-validation was used to estimate the optimal number of PLS components.

All data analyses and figures were performed/generated using MATLAB R2020a (The MathWorks, Inc.).

Using the criteria defined in the Materials and Methods section, we identified 59 HPV16⁺ tumors, 238 nonviral tumors, and 50 normal samples with available methylation data from the TCGA PANCAN HNSCC cohort, listed in Supplementary Table S1. Other virus types detected included: HPV33, n = 8; HHV5, n = 6; HPV35, n = 3; HHV1, n = 1; CYMV, n = 1; HBV, n = 1; HHV4, n = 1; and HPV56, n = 1. There were an additional 171 samples with low counts (<100) and 31 with missing information. All eligible head and neck cancers were included in this study. Of this group, a total of 56 of 297 cases could be distinguished as most likely oropharyngeal by TCGA-defined anatomical sites (e.g., base of tongue, tonsil, oropharynx; Supplementary Table S1).

A PCA model using Illumina 450 methylation array data from all 370 samples with known viral status and all probes with less than 20% missing values (n = 484,933) was performed. The score plot for the first two PCA components (PC1 and PC2), explaining 19.8% and 10.3% of the variation, respectively (Fig. 1A), revealed a clear separation of the normal samples (blue circles) from all tumor samples. This separation is seen in both PC1 and PC2, suggesting that the largest contribution of change in methylation levels arises from tumor versus normal differences. A plot of the third and fourth PCA components (Fig. 1B), accounting for 5.9% and 4.2% of the variation, demonstrated a separation of the HPV16⁺ tumor samples (red diamonds). HPV33⁺ (pink left triangles) and HPV35⁺ tumor (green circles) samples are noted to overlap with HPV16⁺ tumors, suggesting that these types all share a similar methylation pattern. That virus type segregations occur in the third and fourth PCA components, indicates that there are major differences in methylation patterns between tumor-associated HPV types.

Using our DMR selection criteria, we identified 130 significant probes (mapping to 57 genes) being differentially methylated between HPV16⁺ and nonviral tumors. Similar comparisons between normal tissue and nonviral tumors generated 96 significant probes, whereas normal samples versus HPV16⁺ tumors yielded 719 significant probes (Fig. 2A, probes listed in Supplementary Table S2A–S2C). The density scatter plot, Fig. 2B, demonstrates a linear relationship between Δ β-value for normal tissue samples versus HPV16⁺ tumors and normal tissue samples versus nonviral tumors. There is a clear vertical shift in the upper right quadrant indicating that there is a greater increase in the degree of methylation for HPV16⁺ tumors as compared with nonviral tumors.

Synaptonemal complex protein 2 (SYCP2, 20q) was among the most significant and representative of genes identified in this analysis that are hypomethylated in HPV16⁺ as compared with both nonviral tumors and normal samples. Figure 3A clearly demonstrates that HPV16⁺ samples (red) show a lower degree of methylation in SYCP2 compared with nonviral tumors (yellow) and normal samples (blue). Five of the 8 CpG-probes are significantly differentially methylated and are all located in the TSS1500-1st Exon region and may thus potentially impact the transcription of SYCP2. Each probe's correlation to corresponding RNA-seq gene expression levels is shown in Fig. 3B which demonstrates a strong negative correlation (r−0.4) for six of these probes. Methylation levels are also correlated between these six probes, as shown in the correlation heatmap (Fig. 3C). The negative correlation between the average methylation level of SYCP2 and gene expression was very high suggesting that expression of SYCP2 is regulated by methylation (Fig. 3D). There is a clear separation of the HPV16⁺ samples in both the methylation and gene expression spaces. By copy-number analysis for SYCP2, nonviral tumors demonstrate some variations falling into Shallow Deletion, Diploid, and Gain categories, but they do not seem to be related to overall expression levels (Fig. 3E). HPV16⁺ tumors samples show a progressive trend in association with SYCP2 expression from Shallow Deletions to Amplifications (Fig. 3F). These results indicate that the high expression of SYCP2 in HPV16 tumors is very likely to be epigenetically regulated by methylation and not by copy-number alteration events. Notably, out of 515 available HNSCC samples from the TCGA PanCancer Atlas HNSC dataset in cBioportal, SYCP2 was found to harbor very few mutations in HNSCC, with only 3 samples showing a truncating mutation and 4 samples harboring a missense mutation, all of unknown significance (36).

Msh homeobox 2 (MSX2, 5q) is an example of a significantly differentially hypermethylated gene in HPV16⁺ samples as compared with nonviral tumors and normal tissue samples (Fig. 4A). Hypermethylation was associated with lower expression of MSX2 in HPV16⁺ tumors. Differentially methylated probes are located in the TSS1500, TSS200, 5′UTR, 1st Exon and in the gene-body regions of MSX2 (Supplementary Fig. S1A–S1C). Expression was impacted by copy-number variation only in nonviral tumor samples (Supplementary Fig. S1D and S1E).

Helicase like transcription factor (HLTF, 2p) exhibited similar methylation characteristics to SYCP2, with hypomethylation in HPV16⁺ tumors and hypermethylation in nonviral tumors and normal tissues (Fig. 4B). All significant probes for HLTF were identified in the TSS1500 region and were negatively correlated with gene expression (Supplementary Fig. S2A–S2C). The expression of HLTF is also higher in HPV16⁺ tumors, indicating that the expression of HLTF is at least in part regulated by methylation. Copy-number variation influenced expression of HTLF in both nonviral and HPV16⁺ tumors (Supplementary Fig. S2D and S2E).

Paired like homeodomain 2 (PITX2, 4q) demonstrated hypermethylation in HPV16⁺ tumors with hypomethylation in normal tissues and most nonviral tumors (Fig. 4C). For PITX2, 70 CpG-probes were identified on the Illumina 450K chip of which 17 were used for calculating the average methylation level for PITX2. (Supplementary Fig. S3A–S3C). Average methylation levels for PITX2 correlated with gene expression which was not impacted by copy-number variations (Supplementary Fig. S3D and S3E).

GRAM domain containing 4 (GRAMD4, 22q) demonstrated hypomethylation in HPV16⁺ tumors and hypermethylation in nonviral tumors and normal tissues (Fig. 4D). Fifteen probes for GRAMD4 were negatively correlated with gene expression and were also correlated to each other (Supplementary Fig. S4A–S4C). GRAMD4 expression was affected by copy-number variations only in nonviral tumor samples (Supplementary Fig. S4D and S4E).

The HPV16⁺ PLS model (59 HPV16⁺ tumors and 238 nonviral tumors; all other virus types were excluded) was calculated using the 200 selected probes as described in the Materials and Methods section and listed in Supplementary Table S3. Cross-validation yielded two significant PLS components. The tumor versus normal PLS model (50 normal samples and 297 tumor samples) was also calculated using 200 selected probes, as listed in Supplementary Table S3. Two significant PLS components were identified by cross-validation. The results from the training of these two PLS models are shown in Fig. 5A. The HPV16⁺ PLS model score is plotted on the x axis, whereas the tumor versus normal PLS score is plotted on the y axis. In this two-dimensional representation, normal tissue, HPV16⁺ tumors, and nonviral tumors clearly separated into three quadrants.

In order to independently validate the PLS models, we applied the PLS model to samples not used in the training sets. We applied the HPV16⁺ PLS model to HNSCC samples from TCGA that were not used for derivation of the PLS model but with known viral status (Fig. 5B). Five of the eight HPV33 tumor samples (pink left triangles) were scored as “HPV16-related” with the remaining 3 closely approaching the 0.5 cut-off value. All three HPV35 tumor samples (green circles) also were scored as HPV16-related. This demonstrates that these specific HPV types are associated with the methylation of a similar set of genes, a concept also suggested by the PCA analysis (Fig. 1B). Notably, HPV16, HPV33, and HPV35 are closely related on the viral phylogenetic tree and fall within the alpha-9 species. All other virus types were scored as being unassociated with an HPV16-related profile.

Tumor T-cell infiltration has been associated with favorable prognosis (37), and HPV16⁺ tumors were associated with increased tumor infiltration by T cells (CD3E, CD4, and CD8B) and B cells (MS4A1/CD20) with superior effector functions as demonstrated by elevated CD247, Granzyme B, and Perforin and IFNγ (Fig. 5C and D). Increased tumor-infiltrating lymphocyte infiltration in HPV16⁺ tumors was associated with elevated levels of key chemokines that drive T-cell recruitment and trafficking (Fig. 5E) and a significant increase in STING expression (Fig. 5F) which is a crucial molecule implicated in tumor immunogenicity (38).

The HPV16⁺ PLS model was applied to the CESC (Cervical squamous cell carcinoma and endocervical adenocarcinoma) dataset from TCGA with similar annotation of HPV status (Fig. 6A). One hundred forty (96.6%) of 145 HPV16⁺ tumors were correctly scored as HPV16-related (red diamonds). Similar to the TCGA HNSCC data, all HPV35 tumors (n = 4) were scored as HPV16-related along with five out of six for both HPV33 (pink left triangles) and HPV31 (blue right triangles) tumors. These three HPV types (HPV35, HPV33, and HPV31) are all alpha-9 viruses. Other virus types that scored predominantly as HPV16-like included: alpha-5 HPV26 1 (1), alpha-6 HPV30 1 (1), and alpha-5 HPV69 1 (1). Interestingly, 25 (80.6%) of the 31 HPV18 tumors and 16 of 16 HPV45⁺ tumors were scored as HPV16-unrelated. This would indicate that methylation alterations associated with HPV18 and HPV45 (both alpha-7 species) are different from those seen with HPV16 and other alpha-9 counterparts (Fig. 6B). To further validate our HPV16⁺ PLS model, we applied it to the GSE38266 HNSCC dataset and demonstrate complete segregation between HPV⁻ and HPV⁺ samples (Fig. 6C). All HPV⁺ samples scored higher than 0.5, whereas two of the HPV⁻ samples were noted to score slightly above 0.5.

We also applied the HPV16⁺ PLS model to all remaining TCGA samples of all cancer types with available corresponding methylation profiling data (Fig. 6D). Most notably, 63 (40.4%) of 156 testicular germ cell tumors (TGCT) were classified as HPV16-related along with 15 (3.2%) of 475 lung adenocarcinomas (LUAD).

Our HPV16⁺ and tumor versus normal PLS models were applied to the Cancer Cell Line Encyclopedia (CCLE). All but one sample were scored correctly as “Tumor” (Fig. 6E), whereas 18 of the cell lines were classified as HPV16-related (Supplementary Table S4). These results suggest that the methylation changes in the selected probes are intrinsically associated with the tumor cell and not due to adjacent nontumor cells (e.g., immune, stromal) in the pathologic specimen. Seven of the HPV16-related cell lines were of cervical cancer origin, whereas six were derived from lung cancers.

Genome-wide alterations in methylation play a critical role in the development of numerous malignancies (39). Notably, it has also become increasingly apparent that these epigenetic events also contribute to the process of HPV-associated neoplastic progression, such as in a substantial number of head and neck malignancies (40). Although a number of studies have characterized methylation differences between HPV⁺ and HPV⁻ head and neck cancers (19–24), a comprehensive evaluation of differences that includes normal head and neck mucosa (allowing for an evaluation of potential epigenetic contributors to carcinogenesis) and HPV16⁺ (rather than all HPV types together) and nonviral HNSCC has been lacking.

Our analyses provide greater insight with respect to the contribution of methylation events in the process of general tumorigenesis as opposed to those that are uniquely associated with HPV-related carcinogenesis. In this study, we leveraged TCGA genome-wide methylation data for normal tissue and HNSCC to generate a PLS model that uniquely classifies: (i) tumor versus normal and (ii) HPV16⁺ tumors versus nonviral tumors. The PLS models provide insights into both shared and unique methylation alterations associated with normal tissue, HPV16⁺ and nonviral tumors in HNSCC. The HPV status and tumor versus normal PLS models were further validated in cervical squamous cell carcinoma (HPV⁺) from additional external datasets.

Our analysis also highlights variable patterns of epigenetic alterations between cancers associated with different HPV types. It is clear from both our PCA and PLS models that tumors infected with HPV16 have similar methylation patterns to tumors infected with HPV31, HPV33, and HPV35 (all alpha-9 species members) as compared with tumors infected with HPV18 and HPV45 (both alpha-7 species; ref. 41).

When comparing HPV16⁺ and nonviral HNSCC, we identified 130 significant probes mapping to 57 genes, of which 4 (SYCP2, MXS2, HLTF, and PITX2) have previously been reported to be associated with HNSCC.

SYCP2 is a testis-specific human gene with a protein product that is involved in chromosomal recombination during meiosis. We observed it to be differentially hypomethylated (with corresponding increase in gene expression) in HPV16⁺ versus nonviral HNSCC. This finding is concordant with several previous studies which have shown higher expression of SYCP2 in HPV-positive HNSCC as compared with premalignant/normal-matched tissue samples, and HPV-negative HNSCC (42, 43). SYCP2, along with another testis-specific gene testicular cell adhesion molecule 1 TCAM1, has been found to be upregulated in HPV-positive HNSCC by E6 and E7, the viral oncogenes of high-risk HPV (4, 42). In a recent study, Tripathi and colleagues found that higher levels of SYCP2 expression were associated with improved survival in HNSCC (43). SYCP2 hypomethylation has been observed among numerous other cancers and may represent a candidate epigenetically-mediated driver of malignancy (43, 44).

The MSX2 (MSH homeobox 2) protein is a transcription factor that plays a role in the regulation of cell proliferation, differentiation, and survival during development. Expression of this gene is low in most normal adult tissues but was found to be high in several human tumor cell lines (45). In colorectal cancer, MSX2 has been shown to have higher expression in tumor compared with normal colorectal tissue, and higher expression in tumor tissue was also negatively associated with prognosis (46). Studies of malignant melanoma and breast cancer tissues, however, demonstrated a positive association between increased MSX2 expression and prognosis (47). In HNSCC, MSX2 has been shown to be differentially hypermethylated in HPV-positive versus HPV-negative tumors (34). Our findings are consistent with these previous studies, with MSX2 being hypermethylated in HPV16⁺ tumors (with a corresponding reduction in gene expression) compared with nonviral HNSCC and normal tissues.

Supplementary Fig. S1A illustrates that not all probes for MSX2 were significantly differentially methylated, especially those in nonpromoter regions. Genes with significant differential methylation in nonpromoter regions, such as the gene body, often showed no gene expression differences in HPV16⁺ versus nonviral HPV tumors. This exemplifies the importance of using gene-level information when examining the location of differential methylation, as this may influence the selection of biologically significant epigenetically-mediated biomarkers.

HLTF (Helicase-like Transcription Factor), which belongs to the SWI/SNF family of proteins, is a DNA helicase involved with chromatin remodeling and DNA-damage repair. Hypermethylation of the HLTF gene promoter has been reported in more than 40% of colon cancers, as well as in gastric and esophageal cancers (48–50). In colon cancers, HLTF hypermethylation and decreased expression has been associated with increased tumor growth and genomic instability, as well as with poor prognosis (49, 50). In HNSCC, Capouillez and colleagues found that HLTF expression was associated with neoplastic progression (48). HLFT expression was also found to predict recurrence in hypopharyngeal (although not laryngeal) SCC (48). Our study appears to be the first to show that HLTF is differentially hypomethylated (with increased expression) in HPV16⁺ versus nonviral HNSCC and normal tissue.

PITX2 (Paired-like homeodomain transcription factor 2) is a transcription factor that plays an essential role in embryonic development, with involvement in tissue-specific cell proliferation, morphogenesis, and development of left-right asymmetry. High PITX2 expression has been found in colorectal cancers and other WNT/beta-catenin aberrant pathway tumors such as esophageal SCC (51). Low methylation/overexpression of PITX2 has been associated with tumor progression in thyroid and ovarian cancer, whereas hypermethylation has been found in acute myeloid leukemia and is associated with increased risk of disease in non–small cell lung cancer. Hypermethylation of PITX2 has been correlated with HPV-positive status, lower tumor grade and improved overall survival in HNSCC (51). Differential hypermethylation of PITX2 in HPV-positive versus HPV-negative HNSCC was demonstrated by Lechner and colleagues (34). Our findings that PITX2 is hypermethylated in HPV16⁺ HNSCC versus nonviral HNSCC and normal tissue are congruent with these previous studies.

GRAMD4 (Gram domain-containing protein 4) has not been previously reported to be associated with HNSCC. Also known as death-inducing protein, this gene is involved in apoptosis. We found that GRAMD4 is hypomethylated in HPV16⁺ HNSCC compared with nonviral HNSCC and normal tissue. This gene has previously been found to mimic the tumor-suppressive function of p53. However, it has also been found to have higher expression in hepatocellular cancer cell lines compared with normal liver cell lines (52), as well as in lung squamous cell carcinoma versus normal lung tissue (53). Our identification of this gene provides impetus for further investigation of its associations with HPV-associated HNSCC and its oncogenic mechanisms of action.

The effects of HPV in HNSCC including response to immunotherapy were recently reviewed by Lechien and colleagues (54). Our observations of signals of increased immune activation in this setting are concordant with previous reports of increased activation of T cells and chemokines (55). Decreased expression levels of CXCL10, CXCL9, GZMB, PRF1, IFNG, and CD8A have been linked to worse outcomes (56). The overexpression of both gene and protein levels of STING in HPV⁺ tumors, as compared with HPV⁻ tumors, is well established (57, 58).

The role of HPV16 in testicular germ cell tumorigenesis remains debatable. Studies investigating the association of HPV and testicular cancer have been sparse. Two meta-analyses by Garolla and colleagues and Heidegger I and colleagues (59, 60) reviewed HPV- and TGCT-related publications and their findings were inconclusive due to the small number of available studies. A total of five studies focused on investigating the association of HPV and TGCTs have been reported. An epidemiologic study by Strickler and colleagues demonstrated that only 5% of testicular cancers were positive for HPV16 (61). Garolla and colleagues demonstrated that in 155 patients with testicular cancer there was a significantly higher prevalence of HPV-infected semen as compared with controls (9.5% and 2.4%, respectively), however, the HPV status of these tumors were not directly investigated (62). Animal studies by Kondoh and colleagues using transgenic mice carrying HPV16 E6 and E7 open-reading frames demonstrated a high incidence of TGCT that were invariably HPV16-positive (63). In contrast to these findings, in a study of in testicular seminomas (n = 61) and normal tissues (n = 23), HPV16 and HPV18 genotypes were not detected (64). Concordantly, Rajpert-De Meyts and colleagues demonstrated that HPV was not detected in a cohort of 22 TGCTs (65). reviewed HPV- and TGCT-related publications and their findings were inconclusive due to the small number of available studies. A total of five studies focused on investigating the association of HPV and TGCTs have been reported. An epidemiologic study by Strickler and colleagues demonstrated that only 5% of testicular cancers were positive for HPV16 (61). Garolla and colleagues demonstrated that in 155 patients with testicular cancer, there was a significantly higher prevalence of HPV-infected semen as compared with controls (9.5% and 2.4%, respectively); however, the HPV status of these tumors was not directly investigated (62). Animal studies by Kondoh and colleagues using transgenic mice carrying HPV16 E6 and E7 open-reading frames demonstrated a high incidence of TGCT that was invariably HPV16-positive (63). In contrast to these findings, in a study of in testicular seminomas (n = 61) and normal tissues (n = 23), HPV16 and HPV18 genotypes were not detected (64). Concordantly, Rajpert-De Meyts and colleagues demonstrated that HPV was not detected in a cohort of 22 TGCTs (65). By application of our PLS model to the TGCT specimens from TCGA, we demonstrated that 40% of TGCTs displayed a clear HPV16-related methylation profile (Fig. 6C). The role of HPV in testicular cancer tumorigenesis remains controversial due to the limited number of studies. Further investigations are clearly warranted to better elucidate this association.

Similarly, although to a lesser degree, we observed that handful of lung adenocarcinoma samples (3.2%) were classified as HPV16-related by our PLS model. Notably Cao and colleagues examined RNA-seq data from TCGA and did not detect any HPV16 or other HPV types in any lung adenocarcinoma samples (25). To examine whether smoking status could be a potential confounder in our TCGA HNSCC cohort, we did examine HPV16 PLS scoring in smokers versus nonsmokers and demonstrated no differences. As such, it remains unclear as to the significance of an HPV16-like methylation profile among lung adenocarcinomas.

As previously noted, it is well established that HPV⁺ HNSCC have more favorable clinical outcomes and response to treatment as compared with HPV⁻ counterparts; however, there remains a small but significant subset of approximately 20% of HPV⁺ cancers that appear more treatment-resistant and with a worse prognosis (15). Our observations suggest that differences in the HPV-driven epigenotype may provide some underlying molecular explanations to the differences in clinical behavior. Furthermore, there is still, nonetheless, heterogeneity among HPV⁺ HNSCC (which are frequently defined by p16-positivity rather than genotyping). Variable HPV type–related methylomic changes may at least, in part, contribute to such differences. Nakagawa and colleagues have reported a high-DNA methylation subtype in HPV⁺ oropharyngeal SCC that was associated with good prognosis (20). Our findings support the notion that HPV infection results in type-specific biologically significant methylomic changes that may account for observed differences in response to treatment and prognosis in HNSCC.

The emerging importance of DNA methylation changes in the development of HNSCC has also led the recognition of their potential as biomarkers for a number of translational clinical applications. From a screening standpoint, there has been significant interest in developing novel/efficient methods of detecting premalignant head and neck or HNSCC at the earliest of stages. There has been evidence to suggest that strategies such as the identification of molecular alterations in oral gargle specimens may hold promise. Indeed, methylation in selected genes have been detectable in such specimens from patients with HPV⁺ HNSCC (66). More recently, the emergence of reliable circulating tumor DNA detection techniques has also allowed for the identification of methylation changes in patient blood. As an example, Misawa and colleagues demonstrated that a panel of 3 methylated genes were detectable in the plasma of patients with HPV⁺ HNSCC (67). As noted, previous studies have suggested that there are important differences in methylation changes between HPV⁻ and HPV⁺ HNSCC that could allow for etiology-specific detection panels for various platforms. However, our findings are relevant in highlighting that such changes actually vary depending on not only HPV infection in general but also, the underlying HPV type and have implications for more accurate detection panels. In addition to applications in screening, specific DNA methylation signatures (e.g., HPV type–specific) may serve as more refined biomarkers for treatment strategies such as de-escalation and patient prognosis. Broad demethylating treatment agents such as 5azaDC have been long used in hematologic malignancies and are the subject of HNSCC-related clinical trials. However, further understanding of epigenetically mediated carcinogenic elements and pathways may yield novel targeted strategies for therapy and/or enhancement of treatment sensitivity.

Translational development of such findings will require further both mechanistic/biologic and broader validating investigations related to the differential methylation in HPV⁻ and HPV⁺ type–specific HNSCC which may, in turn, lead to the identification of novel strategies for prevention, early diagnosis, surveillance, and personalized treatment strategies.

A. Berglund, E.M. Siegel, and R.M. Putney report grants from NCI/NIH and other support from University of Tennessee Health Science Center during the conduct of the study. S.A. Eschrich reports grants from NCI during the conduct of the study as well as personal fees and other support from Cvergenx, Inc. and personal fees from Cancer Medicine Journal outside the submitted work; in addition, S.A. Eschrich has a patent for 9,846,762 issued and licensed to Cvergenx, a patent for 8,660,801 issued and licensed to Cvergenx, a patent for 8,655,598 issued and licensed to Cvergenx, and a patent for 7,879,545 issued and licensed to Cvergenx. S. Kim reports grants from Bristol Myers Squibb, AstraZeneca, and NIH/NCI during the conduct of the study. D.N. Hayes reports serving on an advisory board for the pharmaceutical company Merck and being a founder, shareholder, and holder of intellectual property licensed to the company GeneCentric Therapeutics. D. Shibata reports grants from NIH/NCI during the conduct of the study. No disclosures were reported by the other authors.

A. Berglund: Conceptualization, data curation, formal analysis, visualization, methodology, writing–original draft. C. Muenyi: Investigation, writing–original draft. E.M. Siegel: Conceptualization, funding acquisition, writing–review and editing. A. Ajidahun: Project administration, writing–review and editing. S.A. Eschrich: Investigation. D. Wong: Writing–review and editing. L.E. Hendrick: Investigation. R.M. Putney: Data curation, writing–review and editing. S. Kim: Writing–review and editing. D.N. Hayes: Writing–review and editing. D. Shibata: Conceptualization, funding acquisition, writing–review and editing.

This project has been supported in part by grant R21-CA185921 (DS/ES) from the NCI. The results generated are, in part, based upon data available from the TCGA Research Network: https://www.cancer.gov/tcga. This work has been supported, in part, by the Biostatistics and Bioinformatics Shared Resource at the H. Lee Moffitt Cancer Center & Research Institute, an NCI-designated Comprehensive Cancer Center (P30-CA076292).

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.