The virome of HPV-positive tonsil squamous cell carcinoma and neck metastasis.

Oropharyngeal squamous cell carcinoma (OPSCC) represents the most common HPV-related malignancy in the United States with increasing incidence. There is heterogeneity between the behavior and response to treatment of HPV-positive oropharyngeal squamous cell carcinoma that may be linked to the tumor virome. In this prospective study, a pan-pathogen microarray (PathoChip) was used to determine the virome of early stage, p16-positive OPSCC and neck metastasis treated with transoral robotic surgery (TORS) and neck dissection. The virome findings of primary tumors and neck lymph nodes were correlated with clinical data to determine if specific organisms were associated with clinical outcomes. A total of 114 patients were enrolled in the study. Double-stranded DNA viruses, specifically Papillomaviridae, showed the highest hybridization signal (viral copies) across all viral families in the primary and positive lymph node samples. High hybridization signals were also detected for signatures of Baculoviridae, Reoviridae, Siphoviridae, Myoviridae, and Polydnaviridae in most of the cancer specimens, including the lymph nodes without cancer present. Across all HPV signatures, HPV16 and 18 had the highest average hybridization signal index and prevalence. To our knowledge, this is the first study that has identified the viral signatures of OPSCC tumors. This will serve as a foundation for future research investigating the role of the virome in OPSCC. Further investigation into the OPSCC microbiome and its variations may allow for improved appreciation of the impact of microbial dysbiosis on risk stratification, oncologic outcomes, and treatment response which has been shown in other cancer sites. Copyright:

INTRODUCTION

Head and neck squamous cell carcinomas (HNSCCs), including oropharyngeal squamous cell carcinoma (OPSCC), are classically associated with risk factors such as tobacco and alcohol use [1]. However, human papillomavirus (HPV) has emerged as one of the primary causes of carcinogenesis for the majority of OPSCCs in North America and Europe [2]. Important clinical differences exist between HPV-positive and HPV-negative OPSCC which ultimately impact treatment decisions. HPV-positive OPSCCs have a higher incidence amongst younger patients with higher performance status, lower tobacco consumption, and higher socioeconomic status [3, 4]. Compared to HPV-negative OPSCCs, HPV-positive OPSCCs have an improved overall survival (OS) and disease-free survival (DFS) [5, 6], which is at least partially related to the increased sensitivity to chemotherapy and radiation therapy [7]. The survival advantage of HPV-positivity persists even after adjusting for confounders [8], suggesting a difference in the tumor biology between HPV-positive and -negative OPSCCs.

An improved understanding of prognostic factors and tumor biology is necessary to improve clinical outcomes and maximize quality of life in long-term survivors of HPV-positive OPSCC. Various studies have demonstrated that there are unique microbiomes associated with certain cancer types [9-14], including in oral and oropharyngeal squamous cell carcinomas [15]. These specific microbiomes may create a microenvironment conducive to oncogenesis, or it may be that the microenvironment in malignancy may promote specific microbiomes. Regardless, the interplay between the host and microbiome impacts the immune system on a local and systemic level that contributes to oncologic outcomes [16, 17]. Our previous work on the mixed group of p16-negative and -positive oral and oropharyngeal squamous cell carcinomas demonstrated that there were distinct viral and other microbial signature patterns associated with this type of malignancy [15]. Amongst the viruses, HPV, specifically HPV16 was the most detected molecular signature which is consistent with the fact that approximately 95% of HPV-positive OPSCCs are caused by HPV16 [18]. Additional viral signatures detected in our study included other HPVs, Herpesviridae, Poxviridae, Retroviridae and Polyomaviridae, and host integration hotspots were identified for HPV16 and JC polyomavirus [15].

In the present study, we utilized a pan-pathogen microarray (PathoChip [19]) to determine the viral microbiome of p16-positive, T1 or T2 tonsil squamous cell carcinomas with no positive nodes or a single <6 cm lymph node metastasis treated with transoral robotic surgery (TORS) and neck dissection. The virome findings were then correlated with clinical data to determine prognostic and diagnostic markers.

RESULTS

Clinical characteristics and outcomes

After excluding inadequate tissue samples, a total of 114 patients were enrolled in the study. There were 23, 22, 23, 22, and 24 samples in the NPR (“negative-node primary”), NLN (“negative lymph node”), PPR (“positive-node primary”), PLN (“positive lymph node”), and control groups, respectively (Figure 1). Table 1 lists the demographic characteristics for patients included in the study. The difference in the proportions of patients who smoked from the negative-node SCC, positive-node SCC, and control groups was significant (71%, 74%, and 36%, p = 0.031). The average age for patients from the negative-node SCC, positive-node SCC, and control groups was significantly different (58.4, 57.0, and 46.4 years, p < 0.001).

Figure 1

Study enrollment.

Shown is the design of the study, including the three cohorts of patients and five cohorts of microbiome analyses. All patients with tonsil SCC underwent TORS radical tonsillectomy with elective neck dissection. Control patients underwent tonsillectomy alone as part of treatment for obstructive sleep apnea. SCC = squamous cell carcinoma; TORS = transoral robotic surgery; OSA = obstructive sleep apnea; ND = neck dissection.

Table 1

Demographic characteristics for all patients

	Tonsil SCC, negative nodes (n = 23)	Tonsil SCC, positive node (n = 23)	Control (n = 24)	P-value
Age (mean (SD), years)	58.4 ± 9.8	57.0 ± 9.8	46.4 ± 10.3	<0.001
Sex (no. male, %)	19 (79%)	20 (87%)	19 (76%)	0.903
Race (no., %)
White	22 (92%)	21 (91%)	23 (96%)	0.779
Black	1 (4%)	2 (9%)	1 (4%)	0.779
Hispanic	1 (4%)	0 (0%)	1 (4%)	0.7
Smokera (no., %)	17 (71%)	17 (74%)	9 (36%)	0.031
Diabetes (no., %)	2 (8%)	3 (13%)	5 (25%)	0.398
Hypertension (no., %)	6 (25%)	8 (35%)	8 (32%)	0.727

SCC = squamous cell carcinoma.

aSmokers included current and former smokers.

Disease and treatment characteristics for patients with SCC are listed in Table 2. There were no statistically significant differences between pathologic factors between the negative-node and positive-node groups. Significantly more patients from the positive-node cohort received post-operative radiation therapy (70% vs 33%, p = 0.01). There were no statistically significant differences in the 2-year OS or DFS between the negative-node and positive-node groups.

Table 2

Disease characteristics for patients with SCC

	Tonsil SCC, negative nodes (n = 23)	Tonsil SCC, positive node (n = 23)	P-value
Primary size (mean (SD), cm)	2.37 ± 0.87	2.53 ± 0.77	0.51
Positive margins (%)	4%a	0%	0.32
LVI (%)	54%	74%	0.16
PNI (%)	0%	0%	0.60
Node size (mean (SD), cm)	n/a	3.42 ± 1.15	n/a
ECS (%)	n/a	0%	n/a
Adjuvant treatment			n/a
Chemotherapy	8%b	9%	0.97
Radiation therapy	33%	70%	0.01
Outcomes
2-year overall survivalc (%)	90%	95%	0.58
2-year disease free survivalc (%)	81%	90%	0.97

SCC = squamous cell carcinoma; LVI = Lymphovascular invasion; PNI = Perineural invasion; ECS = Extracapsular extension.

apositive margin treated with chemoradiation including cetuximab.

bincludes 1 patient who received cetuximab.

cpatients with insufficient follow up were excluded from analysis.

Viral signatures in p16-positive tonsil squamous cell carcinoma

Various DNA and RNA viral signatures were associated with the tonsil cancer groups (Figure 2). Double-stranded DNA viruses had the highest ratio amongst the viral signatures for all cohorts (p < 0.001), followed by single-stranded RNA (Figure 2A). There were 101 virus families with signatures. A hybridization signal index of > 5 was used to select virus families with relatively higher intensity signatures, of which there were 43. High hybridization signals were detected from Siphoviridae, Myoviridae, Reoviridae, Polydnaviridae, and Baculoviridae in most of the cancer specimens, including NLN (Figure 2B, Supplementary Table 1). Viral sequences belonging to the double-stranded DNA virus Papillomaviridae showed the highest hybridization signal across all viral families in the PPR, PLN, and NPR samples (Figure 2B, Supplementary Table 1). HPV detection by PathoChip was 95.7%, 100%, 100%, 86.4%, and 20.8% of PPR, NPR, PLN, NLN, and control specimens, respectively.

Figure 2

Viral signatures for tonsil cancer cohorts and controls by PathoChip screen.

(A) Proportion of different viral signatures detected significantly in the tonsil cancer samples represented in bar graphs divided by viral types. (B) High level signatures of viral families for cancer cohorts, with the average hybridization signal as a bar graph and the prevalence as dots. (C) High level HPV signatures divided by subtype with the total hybridization signal as a bar graph and the prevalence as dots are shown in red and blue, respectively. NPR (“negative-node primary”), PPR (“positive-node primary”), NLN (“negative lymph node”), PLN (“positive lymph node”), HPV (human papillomavirus), *=p < 0.01, **= p < 0.001, ***= p < 0.0001.

Fourteen HPV subtype signatures were detected with high hybridization signals amongst all cohorts, with the two high-risk subtypes HPV16 and HPV18, detected (Figure 2C). Across all HPV signatures, HPV16 had the highest average hybridization signal intensity and prevalence (Figure 2C). HPV16 was present in 73.9%, 95.7%, 95.5%, 0.00%, and 16.7% of PPR, NPR, PLN, NLN, and control specimens, respectively. Additionally, HPV18 maintained a high prevalence amongst the cohorts and was present in 65.2%, 69.6%, 95.5%, 45.5%, and 0.00% of PPR, NPR, PLN, NLN, and control specimens, respectively. Low-risk HPV6b had a prevalence of 34.8%, 34.8%, 31.8%, 0.00%, and 0.00% in PPR, NPR, PLN, NLN, and control specimens, respectively.

Due to the low number of deaths in each cohort, no significant probes were identified on Kaplan–Meier analysis of OS and DFS. Specifically, we could only perform survival analysis in the NPR group because of the limited number of deaths in the other cohorts. HPV16 probes had the highest hybridization signal indexes, prevalence, and number of probes within the NPR group. Comparisons of overall survival for high and low hybridization signal indexes for HPV16 probes in the NPR group were not significant (p = 0.082). There were no significant differences seen on evaluation of pathologic factors (tumor stage (T1 versus T2), PNI, LVI and tumors with a positive lymph node (NPR versus PPR)) in different hierarchical clusters of cancer samples for both HPV16 and HPV18 probes (Figure 3).

Figure 3

Hierarchical clustering of tonsil squamous cell carcinoma cohorts based on viral signature detection pattern with associated pathological features (perineural invasion (PNI), lymphovascular invasion (LVI), and tumor stage (T stage)).

Hierarchical clustering for HPV16 (A) and HPV18 (B) viral probes are represented as heat maps for each cohort. Clustering was performed by R program using Euclidean distance, complete linkage and non-adjusted values. Clustering of the samples using NBClust software [Calinski and Harabasz index, Euclidean distance, complete linkage]. Chi-square test was applied and showed no significant differences of proportions of tumor stage (T1 versus T2), PNI, and LVI in different hierarchical clusters of cancer samples. NPR (“negative-node primary”), PPR (“positive-node primary”), NLN (“negative lymph node”), PLN (“positive lymph node”), CTRL (“control”), HPV (human papilloma virus), PNI (perineural invasion), LVI (lymphovascular invasion), T stage (tumor stage).

DISCUSSION

We previously described the microbiome of oral and oropharyngeal squamous cell carcinomas [15]. In the current study, we evaluated the virome for a more homogeneous group of p16-positive, T1 or T2 tonsil squamous cell carcinomas with no positive lymph nodes or a single positive level 2 lymph node <6 cm, analyzing both the primary tumor and the lymph node specimens. To our knowledge, this is the first report of the virome of p16-positive tonsil squamous cell carcinoma. Our motivation for evaluating both pathologically positive and negative lymph node specimens in addition to primary tumors was to improve our understanding of the microenvironment necessary for or induced by regional metastasis. Furthermore, we attempted to determine clinical correlations with specific viral signatures.

In accordance with the NCCN treatment guidelines at the time of treatment for p16-positive T1-2, N0-N1 [20], patients in the study were treated with TORS primary resection and ipsilateral neck dissection followed by radiation or chemotherapy based on multidisciplinary review of pathology. We found no significant survival difference between N0 and N1 patients, which is consistent with prior studies [21, 22]. Haughey et al. demonstrated that increasing pathological N-classification was not associated with worse OS; in fact, N0 cases demonstrated worse OS than N1 in their analysis [21]. Based on prior studies of tobacco and HPV-associated OPSCC [6], the high prevalence of smoking (>70%) in our cohorts would be expected to negatively impact the clinical outcomes.

HPV infection is typically acquired through sexual contact, and the prevalence of oral HPV infection increases with the number of lifetime sexual partners and cigarettes smoked per day [23]. Prior larger studies have demonstrated that healthy controls have a prevalence of oral HPV infection of around 5% which is lower than our control group with 21%, possibly due to the higher rates of smoking in our cohort [23-25]. HPV16 is the most commonly detected carcinogenic HPV in oral samples of healthy patients [24], with prevalence around 1% [23-25]. We demonstrated a higher control prevalence of HPV16 (17%), but a low prevalence of the carcinogenic HPV18 (0%).

Though typically cleared by the immune system, in some individuals, the viral DNA integrates into the host genome [4] and can lead to OPSCC [2]. Persistent HPV16 DNA after treatment of HPV-positive OPSCC has been associated with worse survival [26]. Immunohistochemistry (IHC) for p16, a protein that increases with inhibition of pRb by E7, is currently the most widely utilized surrogate HPV biomarker [4, 27, 28], and is the method utilized at our institution for determining HPV-positivity in OPSCC. Interestingly, several studies have demonstrated some discordance between p16-positivity by IHC and the presence of HPV DNA via PCR [29-31]. Our results indicate a high concordance between p16-positivity and the presence of HPV signatures detected by PathoChip (100% prevalence in NPR and PLN and ~96% in PPR). Our previous work on primary oral and oropharyngeal squamous cell carcinoma demonstrated an average HPV16 prevalence of 98% [15], which was higher than the prevalence for the primary oropharyngeal tumor samples in the current study (~85%).

In contrast to our prior analysis of the mixed group of p16-negative and -positive oral and oropharyngeal squamous cell carcinomas [15], the current study contains a more homogeneous group of p16-positive OPSCC. Our previous study on p16-negative and -positive oral and oropharyngeal squamous cell carcinomas demonstrated high hybridization signals for Herpesviridae, Poxviridae, Retroviridae, Reoviridae, Orthomyxoviridae, and Polyomaviridae [15]. Similarly, the current study identified high hybridization signals for Baculoviridae, Reoviridae, Siphoviridae, Myoviridae, Polydnaviridae. Approximately 15% of human cancers are known to be caused by viruses [32], the majority of which are double-stranded DNA viruses [33]. The viruses with the highest hybridization signals in the cancer cohorts were predominantly double-stranded DNA viruses, with some exceptions including the double-stranded RNA virus Reoviridae. We suggest that the viruses identified in this study may represent signatures for p16-positive OPSCC.

The finding that Reoviridae had high hybridization signal in all cancer cohorts is particularly interesting because this virus is known to have strong specificity for cancer cells, particularly those with increased Ras activity [34-36]. Furthermore, Reovirus is known to have potent anti-tumor activity in many tumor models including HNSCC [37-39]. It is entirely possible that co-infection of Reovirus is increased in HPV-positive OPSCC compared to HPV-negative OPSCC, ultimately impacting the differences seen in survival. If true, Reovirus signatures may prove to be important biomarkers for OPSCC prognosis, helping to explain the heterogeneity in clinical responses and resistance to therapy in some tumors. Additional biomarkers for prognosis could aid in decisions regarding adjuvant therapy and treatment de-escalation to avoid morbidity. Regardless, this hypothesis is speculative and outside of the scope of the current study.

The viromes for the negative-node and positive-node primary tumors (NPR and PPR) had similar distributions of viral signatures. It remains to be determined if certain viruses impact the primary tumor’s ability to metastasize regionally. We did not identify any correlations between viral signatures and survival or patterns of failure; however, the number of recurrences and deaths in our study was likely too small to detect any differences. Similarly, the association between any viral probes and pathologic features will need to be explored with a larger scale study. Further work may determine that certain viruses induce alterations in gene expression to impact oncogenesis or cancer progression. These alterations may manifest in pathologic features such as tumor invasion of lymphatics and nerve or spread outside of lymph nodes; conversely, changes may not be appreciated microscopically and instead exist at the molecular level.

Although associations between specific viral signatures and p16-positive OPSCC disease characteristics remain to be determined, this study’s description of the OPSCC virome serves as a foundation for additional micobiome analyses. Several separate studies in the cancer literature suggest that host microbiome interactions effect oncologic outcomes [16, 17, 40, 41]. Specifically, it has been shown that the breast microbiome interacts with host cells to impact signaling pathways, ultimately modulating breast cancer growth and metastatic progression [40]. Within the pancreatic cancer literature, research has demonstrated that differences in microbiome diversity and the presence of specific microbial signatures predict long-term pancreatic adenocarcinoma survival [41]. Moreover, the tumor microbiome can be therapeutically altered in pancreatic adenocarcinoma with fecal microbiota transplantation, leading to changes in tumor growth and immune system response. Like the above-mentioned studies, additional work on the OPSCC microbiome may determine that cancer microbial dysbiosis impacts risk stratification and predicts oncologic outcomes and treatment response.

Important limitations exist with the current study. For example, the analysis and generalizability may be impacted by the control group having approximately double the prevalence of smoking compared to the estimated U. S. prevalence of 15.1% [42]. This increased smoking prevalence is likely related to the fact that patients with OSA have higher rates of smoking compared to those without OSA [43]. Additionally, the control group had about half the prevalence of smoking compared to the cancer group which may have impacted the analysis.

In summary, specific viruses, including HPV16, are known to impact the tumor biology and clinical behavior of OPSCCs. The virome of HPV-positive OPSCC primary tumors and neck lymph nodes include the virus families Papillomaviridae, Herpesviridae, Baculoviridae, Reoviridae, Siphoviridae, Myoviridae, and Polydnaviridae. Additional studies are necessary to determine if the identified viral signatures correlate with tumor behavior.

MATERIALS AND METHODS

Experiments were performed in accordance with relevant guidelines and regulations with appropriate licensing and approvals by institutional committees at Perelman School of Medicine, University of Pennsylvania. The study was approved by the University of Pennsylvania institutional review board (protocol number 828709). The Cancer Center’s Clinical Trials Scientific Review and Monitoring Committee determined that this study met the requirements for review exemption.

PathoChip design

The PathoChip Array design has been previously described [11, 19]. The array is comprised of 60,000 probe sets representing reference sequences of all known viruses, and human pathogenic bacteria, fungi, and parasites in Genbank. The PathoChip array contains probes which target genomic regions that are conserved between viruses in the same family and unique probes for each specific virus or microorganism. The arrays are manufactured as SurePrint glass slide microarrays (Agilent Technologies Inc.) with 8 replicate arrays per slide. Each probe is a 60 nucleotide (nt) DNA oligomer that targets multiple genomic regions of viruses, prokaryotic, and eukaryotic microorganisms. Accession annotations are available in the Gene Expression Omnibus (http://www.ncbi.nlm.nih.gov/geo/) [17].

Study samples

Tonsil cancer and lymph node specimens were obtained from patients with p16 positive, T1 or T2 tonsil squamous cell carcinoma (staged with the American Joint Committee on Cancer (AJCC), 7th edition, TNM staging [44] and consistent with the current 8th edition [45]) treated with TORS and therapeutic or elective neck dissection. Only patients with no positive lymph nodes (N0) or a single positive level 2 lymph node less than 6 cm were included in the study (N1 by AJCC, 7th and 8th editions [44, 45]). All tumors were p16 positive determined by surgical biopsy or resection per the guidelines of the University of Pennsylvania Pathology Department and prior studies [46]. Specifically, samples are classified as p16 positive when p16 immunohistochemistry has greater than or equal to 70% of tumor nuclear and cytoplasmic staining.

Cases were divided into two groups based on the presence or absence of neck metastases after pathology review (Figure 1). For patients with no neck metastases, the primary tumor and matched negative neck lymph node specimens were denoted as NPR (“negative-node primary”) and NLN (“negative lymph node”), respectively. Alternatively, for patients with neck metastasis, the primary tumor and matched positive neck lymph node specimens were denoted as PPR (“positive-node primary”) and PLN (“positive lymph node”), respectively. The control group consisted of tonsil specimens from patients with obstructive sleep apnea (OSA) and no history of recurrent or chronic tonsillitis who underwent tonsillectomy as part of their treatment for OSA.

Sample preparation and microarray processing

PathoChip screening utilized DNA and RNA extracted from formalin-fixed paraffin-embedded (FFPE) tissues as described previously [10, 19]. De-identified FFPE samples were received as 10 μm sections on non-charged glass slides. Tumor and control tissue were read by the pathologist at the Department of Pathology at the Hospital of the University of Pennsylvania, Philadelphia, PA. The samples were cut in a sterile fashion using a sterilized microtome.

The PathoChip screen procedure has been described in prior work [9–11, 15]. Briefly, DNA and RNA are extracted from the FFPE samples. Quality of extracted nucleic acids was determined by agarose gel electrophoresis and A260/280 measurements. Whole genome and transcriptome amplification (WTA) using the TransPlex Complete Whole Transcriptome Amplification Kit (Sigma-Aldrich, St. Louis, MO) was performed on extracted DNA and RNA samples (50 ng of DNA and 100 ng RNA) in accordance with manufacturer protocol. Reference human DNA and RNA was extracted from the human B cell line, BJAB (obtained from ATCC). The cell lines were cultured for less than 6 months and 15 ng of DNA and RNA were used for WTA. The human DNA and RNA served as a reference for determining cross-hybridization of probes to human sample amplified genome. After purifying the WTA products (PCR purification kit, Qiagen, Germantown, MD, USA), quality was checked using agarose gel electrophoresis and 1 μg of the amplified products from the cancer and control tissues were labelled with Cy3; the human reference was labelled with Cy5 (SureTag labeling kit, Agilent Technologies, Santa Clara, CA). All labeled specimens were purified and hybridized to the PathoChip microarrays, as previously described [11, 19]. For each PathoChip array, a Cy3 labelled sample and a Cy5 labelled reference was combined and hybridized together in constant rotation for 40 hr at 65°C. The slides were then washed and scanned for visualization using an Agilent SureScan G4900DA array scanner.

Microarray data extraction and statistical analysis

The microarray data extraction and analyses techniques have been described previously [9–11, 19]. Briefly, Agilent Feature Extraction software was used for extraction of raw data from the microarray images, and the R program was used for data normalization and analyses [47]. The microarray screen data are available in Gene Expression Omnibus (Accession No. GSE111648). We calculated scale factors using the signals of green and red channels for human probes. Scale factors are the sum of green and sum of red signal ratios [∑(g)/∑(r)] of human probes. We used scale factors to obtain normalized signals for all other probes. For all probes except human probes, normalized signal is log2 transformed of green signals / scale factors modified red signals (log2 g – scale factor * log2 r). On the normalized signals, t-test is applied to select probes significantly present in cancer samples by comparing cancer samples versus controls and to select probes significantly present in the tonsil cancer samples versus the controls. The cut-off for significant detections in cancers versus the controls was log2 fold change > 1 and adjusted p value (with multiple testing corrections) < 0.05. Prevalence was calculated by counting the number of cancer cases with hybridization signal greater than the average signal of the negative control probes, and represented as a percentage.

Model-based analysis of tiling arrays (MAT) was used for detection of positive regions in the metagenome for each sample using the rMAT software [48]. This technique uses a sliding window to scroll through the entire metagenome of the array to detect positive hybridization signal (hot spots (highly expressed gene region) on the genome) [11, 19]. Analyses at the individual probe level (for both specific and conserved probes) and at the family level (taking into account all the probes per family) were performed.

The cancer samples were also subjected to unsupervised hierarchical clustering, based on the detection of microbial signatures in the samples (average hybridization signal per viral family or microbial genus), using the R program (Euclidean distance, complete linkage, non-adjusted values [47, 49], and the clusters were validated by Calinski and Harabasz index, which is implemented in R package as NbClust [50]. Calinski and Harabasz index is a cluster index that maximizes inter-cluster distances and minimizes intra-cluster distances. We calculated the possible cluster solution that would maximize the index values to achieve the best clustering of the data. The significant differences between the clusters observed by these methods were determined using two-sided t-test. The ANOVA test was carried out to find the common signatures significantly present in all the clusters.

Clinical data and analysis of oncologic outcomes

Patients were followed prospectively following their initial diagnosis and treatment. Clinical information was obtained by an honest broker at the Hospital of the University of Pennsylvania who had access to the patient identifiers, but not the microbiome data. The clinical data was delivered to the statistician who matched the clinical information to the corresponding samples. The clinical information obtained for all samples included patient age, gender, ethnicity, smoking history, and presence of diabetes or hypertension. Additional data was collected for cancer patients including tumor size and stage, node size and stage, presence of lymphovascular invasion (LVI), perineural invasion (PNI), extracapsular spread (ECS), positive margins, and adjuvant treatment (radiation and/or chemotherapy). The outcome measures included 2-year OS and 2-year DFS. Chi-square test and t-test were applied to test whether there were any significant differences disease characteristics in the tonsil cancer patients with and without positive neck disease.

We performed Kaplan–Meier analysis for viral probes to identify correlations with OS and DFS. Chi-square test was applied to test whether there were any significant differences of proportions of pathologic factors in different hierarchical clusters of cancer samples.