MicroRNA-Related Genetic Variants Associated with Survival of Head and Neck Squamous Cell Carcinoma

Head and neck squamous cell carcinoma (HNSCC) is the sixth most common cancer worldwide (1). Treatment can be effective in early-stage disease, but is much less so for advanced disease, and approximately two thirds of patients are diagnosed with advanced-stage disease (2). With the exception of recently approved immune checkpoint inhibitors as second-line therapy (3), treatment paradigms have improved little over the past three decades (4), and involve a combination of radiation, chemotherapy, and disfiguring surgery with severe morbidity (5). Alcohol consumption, tobacco use, and infection with high-risk subtypes of human papillomavirus (HPV) are well-established risk factors for HNSCC (6, 7). Several genetic risk factors have been identified (8–10); however, prognostic factors of HNSCC survival are less well studied. Recognizing this fact, in 2017, the International Head and Neck Consortium initiated a survival working group to address gaps in outcome research for HNSCC. Current models of HNSCC survival are based on histopathologic and staging criteria (11). Alcohol consumption, tobacco use, and HPV infection have also been associated with overall survival of HNSCC (12–14). In addition, genetic predictors of HNSCC survival have been identified, uncovering variants involved in DNA repair (15–17), cell-cycle control (18), and angiogenesis (19). Single nucleotide polymorphisms (SNPs) in miRNA target sites, miRNA seed sites, and miRNA processing proteins represent a class of genetic variation with gene-regulatory potential and have been associated with survival of several cancers (20), including HNSCC (21–24). Such genetic variation can alter miRNA-dependent gene regulation resulting in changes in cellular behavior (20). Although predictive utility of miRNA expression for HNSCC survival has been extensively studied (25, 26), genome-scale studies of the association between miRNA-related SNPs (miR-SNPs) and HNSCC survival are lacking. To date, associations between miR-SNPs and outcomes in HNSCC have used candidate SNP approaches (21, 22, 24, 27, 28). To address this gap in knowledge, we conducted a genome-scale investigation of the association between miR-SNPs and survival in HNSCC.

DNA was collected from study participants of two population-based case–control studies of HNSCC. The Massachusetts (n = 904 cases, n = 1,051 controls) and M.D. Anderson study populations (n = 1,338 cases, n = 1,356 controls) have been described previously (23, 29, 30). For the Massachusetts study, incident cases of HNSCC were identified at nine medical institutions in the Boston metropolitan area between 1999–2003 (phase I) and 2006–2011 (phase II). An independent study pathologist confirmed pathology report histology. Clinical information was obtained through medical chart review while demographic and exposure data were collected using self-administered questionnaires, reviewed by trained study interviewers. Publicly available databases were used to obtain overall survival times for cases in the Massachusetts study. In the M.D. Anderson study, incident cases of HNSCC were recruited at The University of Texas M.D. Anderson Cancer Center between 1996 and 2008. All subjects were treated for curative intent at this institution. Demographic and exposure data were also collected through completion of a self-administered questionnaire. Although both categorical and continuous measures of tobacco and alcohol use were collected, continuous measures were not available for most subjects from the M.D. Anderson study. Medical record review was used to obtain overall survival times, defined as the time between first appointment and death from any cause, or last contact date. Date of first appointment was used to define survival times, as patients are referred to the cancer center from across the United States after initial diagnosis. Importantly, subjects remained untreated until after their first appointment at the M.D. Anderson Cancer Center. International Classification of Disease, Ninth Revision (ICD-9) codings and pathologic analyses were used to assign cases to anatomic subdivisions within the head and neck oral cavity, pharyngeal, or laryngeal. These groups were used to conduct analyses stratified by tumor site. HNSCC cases were classified as ICD9 codes 141, 143–146, 148, 149, and 161. Site-specific ICD9 codings used are as follows: oral cavity—ICD9 codes 141.1–141.5, 141.8, 141.9, 143–145.2, 145.5–145.9, 149.8, and 149.9; pharyngeal—ICD-9 codes 141.0, 141.6, 145.3, 145.4, 146, 148, 149.0, 149.1; and laryngeal—ICD9 code 161.

Whole blood or buccal cell DNA from discovery phase (Massachusetts study) was extracted using the QIAamp DNA Mini Kit (Qiagen) and genotyped using the standard Axiom miRNA Target Site Genotyping Array (Affymetrix). This array is specifically designed to interrogate miRNA-related genetic variation, and has previously been described (23). The array contains SNPs and insertions/deletions (indels) in approximately 238,000 miRNA-related loci including miRNAs, miRNA regulatory regions (such as promoters), miRNA processing proteins, and miRNA target sites. Five online databases containing data regarding miRNA-related genetic variation were used to construct the majority of content on the array; PolymiRTS (31), dPORE (32), Patrocles (33), miRNASNP (34), and microRNA.org (35). We performed quality control (QC) according the Axiom Best Practice Genotyping Analysis Workflow using data from all available cases and controls. Twenty-six samples with either a call rate <97% or a Dish QC (DCQ) metric below a threshold of 0.82 were excluded from the analysis. Minor allele frequencies (MAFs) were calculated using all available discovery phase study subjects. Analysis was restricted to variants with MAF ≥5%, leaving 40,286 markers. Genotyping of validation phase (M.D. Anderson study) subjects was performed using the MassARRAY iPLEX gold assay (Sequenom) for 114 variants selected from discovery analyses. Validation phase subjects with a call rate of <95% were removed from the analysis. Two variants that were monomorphic in the validation population and six variants with a call rate <95% were removed during QC. An additional six variants deviating from Hardy–Weinberg equilibrium (HWE) with P < 1 × 10⁻³ in healthy Caucasian control subjects were removed. Two heterozygote genotypes (AG, n = 4; CG, n = 856) were detected for rs2450137 in M.D. Anderson study subjects, in addition to CC (n = 113) and GG (n = 1,991) genotypes. Given the lack of subjects homozygous for rs2450137 with the A allele, subjects with AG genotypes were omitted from the final analysis for this SNP. All genotypes provided are from the + strand of hg19.

Multivariable cox proportional hazards regression was used to test the association between miR-SNP genotype and overall survival. R Statistical Software version 3.3.0 was used to calculate HRs, 95% confidence intervals (CIs), and P values, assuming a dominant mode of inheritance. Regression models for discovery phase analyses were adjusted for potential confounders, including age (≤50, >50–≤60, >60–≤70, >70), sex (male, female), race (Caucasian, other), HPV serology (positive/negative for listed serological subtypes), alcohol consumption (lifetime average number of drinks/week), tobacco use (lifetime pack-years), and tumor stage (low stage; I and II vs. high stage; III and IV, due to the clinical and prognostic similarity among early- and late-stage disease). Data from all available study subjects (cases and controls) were used to calculate population quartile distributions for alcohol consumption (≤1, >1–≤6, >6–≤31, >31 average number of drinks/week) and tobacco (≤2, >2–≤6, >6–≤14, >14 lifetime pack-years). Cases were ordered into discrete groups for the analysis, based on these quartile distributions. Complete covariate data were available in 847 subjects from the Massachusetts study (Table 1) and were used in the discovery phase analyses. Evidence of population stratification was assessed through calculation of genomic inflation factor λ using the GenABEL R-package. Markers were selected for validation genotyping based on P values from adjusted analyses and were pruned based on pairwise LD patterns (from all available study subjects) to select the two SNPs in lowest pairwise LD in each gene. In total, 123 associations (35 for overall HNSCC, 35 for oral cavity cancer, 29 for pharyngeal cancer, and 24 for laryngeal cancer) across 114 variants (due to overlap between some variants selected in the overall and site-specific analyses) were selected for replication testing. Of the markers selected for validation genotyping, only one variant, rs1971475, deviated from HWE in control subjects from the discovery population, and was removed from the final joint population analyses. After the QC steps described in the Genotyping section, 31, 20, 28, and 19 variants were left to be tested for replication of effect in the overall HNSCC, oral cavity cancer, pharyngeal cancer, and laryngeal cancer analyses, respectively. Analyses of validation phase subjects and the joint population analyses were adjusted for age, sex, race, smoking status (current, former, never), and tumor stage (as above). Complete data on these covariates were available for 1,236 cases from the M.D. Anderson population (Table 1), and were used in validation phase analyses. In the joint population analysis, genotyping data across studies were pooled and adjusted for age, sex, race, smoking status (current, former, never), and tumor stage. Thirty-seven cases from the Massachusetts study were missing smoking status (current, former, never) and were excluded from the joint analysis, leaving 2,046 subjects for the joint population analysis. To correct for inflated type I error resulting from multiple testing, we determined replication as those variants with P < 1 × 10⁻³ in the joint analysis, a lower joint analysis P value than discovery phase P value, and a consistent predicted direction of effect across both populations.

To identify miRNA–mRNA interactions potentially disrupted by genetic variants associated with overall survival in HNSCC, we utilized data from several publicly available databases containing predicted miRNA target sites or predicted effects of genetic variation upon miRNA–mRNA interaction. We have previously described the analytic framework used for this approach (23). Briefly, to identify high confidence–predicted miRNA target sites, we intersected all available miRNA target site predictions from the microRNA.org database (August 2010 Release) with UCSC Genome Assembly hg19 coordinates for survival-associated SNPs. miRNA target site predictions in the microRNA.org database were made using the miRanda algorithm and scored for meaningful downregulation of mRNA targets by the mirSVR algorithm (more negative scores suggest greater potential for meaningful downregulation of its target). High-confidence predictions were considered those with a mirSVR score of <−0.1. Percentiles for mirSVR scores were calculated to help interpret the potential impact of a predicted miRNA target site. It should be noted that because more negative mirSVR scores are suggestive of greater potential for downregulation of a target mRNA, those with the lowest percentile ranks (i.e., 1st, 2nd, 3rd, etc.) represent those predictions with the most likely functional impact. To identify miRNA–mRNA interactions potentially disrupted by survival-associated genetic variants, we mined data from three online databases; PolymiRTSv3.0 (31), miRNASNP (34), and MirSNP (36). PolymiRTS (http://compbio.uthsc.edu/miRSNP/) uses context+ scores (generated by the TargetScan algorithm) to rank the predicted effects of genetic variation upon miRNA target site binding. TargetScan context+ scores predict binding of a miRNA over a 3′UTR. Greater likelihood of miRNA–mRNA disruption for a given variant is indicated by more negative context+ scores. The difference between context+ scores predicted for sequences containing each allele of a genetic variant provides a measure of disruption to predicted miRNA binding. We calculated percentile ranks from all available context+ score differences available from the PolymiRTS database. Lower percentile ranks indicate miRNA–mRNA interactions that are more likely to be disrupted by genetic variation. miRNASNP (http://www.bioguo.org/miRNASNP/index.php) utilizes TargetScan and miRanda algorithms to identify genetic variants likely to impact miRNA–mRNA interaction. MirSNP (http://bioinfo.bjmu.edu.cn/mirsnp/search/) determines predicted miRNA target sites for sequences containing each allele of a given genetic variant using the miRanda algorithm and assigns a predicted impact of genetic variation (enhance, decrease, create, or break) at this locus based on the results.

Clinical/demographic data and miRNA-seq data for subjects from The Cancer Genome Atlas (TCGA) HNSCC project were downloaded from at the Genomic Data Commons Data Portal (accessed 29 May, 2017; refs. 37, 38). Excluding subjects whose samples were from recurrent tumors or received neoadjuvant treatment, 523 primary HNSCCs (oral cavity cancer: n = 314; pharyngeal cancer: n = 86; laryngeal cancer: n = 119; other: n = 4) and 44 adjacent normal tissue samples were available for analysis. ICD-10 codes available in the TCGA dataset were converted to ICD-9 codings to assign TCGA subjects to oral cavity, pharyngeal, and laryngeal cancer groups, for concordance with the subsite classification used in our study, outlined in the “Study participants” section. Four subjects were removed from the TCGA dataset as they did not align with the subgroup classification criteria of our study. Conversions for ICD-10 to ICD-9 codes are summarized in Supplementary Table S1. To assess expression of miRNAs predicted to target regions containing miR-SNPs, average mapped reads/million from isoform-specific miRNA-seq data were calculated across all available subjects for either tumor or normal tissue samples, log₂ transformed, and converted to percentile ranks. For survival analyses, log₂ transformed FPKM (fragments per kilobase million) expression values were calculated using the overall miRNA expression (not isoform-specific) and RNA sequencing (RNA-seq) data. A total of 303 oral cavity cancer cases with available overall miRNA expression and complete covariate data for the variables of interest were used to test the association between survival time and MIR100HG-derived miRNAs. For mRNA expression analyses, RNA-seq data were available for 478 subjects with complete covariate data. BLID expression was undetectable in the majority of subjects with tumors of the oral cavity. Measurable BLID expression and complete covariate data were only available in 37 (12%) oral cavity cancer subjects; therefore, the association between BLID expression and overall survival was not tested. SH3BP4 expression and complete covariate data were available for 108 laryngeal tumors. Cox proportional hazards regression models adjusting for age at diagnosis, sex, race (white vs. nonwhite), and tumor stage (low grade; I + II vs. high-grade; III + IV). Gene expression was modeled as a continuous variable.

R Code used to perform the in silico and TCGA analyses has been deposited in the “HNSCC-miRSNP-Survival” repository on GitHub (https://github.com/Christensen-Lab-Dartmouth). Genotype and phenotype data will be made available under Database of Genotypes and Phenotypes (dbGAP) accession phs001668.v1.p1.

To identify miR-SNPs associated with HNSCC-specific survival, HNSCC cases from a population-based case–control study of HNSCC were genotyped using the Affymetrix Axiom miRNA Target Site Array, which interrogates approximately 238,000 genetic variants in miRNA target sites, miRNA genes, and genes in the miRNA biogenesis pathway. 40,286 genetic variants with a MAF ≥5% (in available cases and controls) were tested for association with overall survival in all HNSCC cases, or tumor site–specific (oral cavity, pharyngeal, laryngeal) analyses. Demographic and clinical characteristics of the 847 cases from the discovery population are in Table 1. We applied Cox proportional hazards regression adjusted for age, sex, race, HPV seropositivity, tobacco consumption, alcohol consumption, and tumor stage. The strongest associations were observed for rs4880198 (overall HNSCC, UAP1L1, HR, 1.54; 95% CI, 1.24–1.92; Supplementary Fig. S1A), rs2292842 (oral cavity cancer, PRDM8; HR, 3.40; 95% CI, 2.06–5.60; Supplementary Fig. S1B), rs72765149 (pharyngeal cancer, ZBED3-AS1; HR, 2.90; 95% CI, 1.79–4.71; Supplementary Fig. S1C), rs1568391 (laryngeal cancer, IRF8, HR, 0.25; 95% CI, 0.13–0.47; Supplementary Fig. S1D). Generally, the strongest associations were observed in the laryngeal cancer analysis, with 74 variants observed with an adjusted P < 1.0E−03 (Supplementary Fig. S1D). In the analyses of overall HNSCC, oral cavity cancer, and pharyngeal cancer, 39, 57, and 45 variants, respectively, had adjusted P < 1.0E−03 (Supplementary Fig. S1A–S1C). We observed minimal evidence of population stratification (Supplementary Fig. S2A–S2D). Association strength by genomic location for all variants is shown in Supplementary Fig. S1, where variants selected for genotyping in the validation phase are highlighted in green.

To validate our findings, we genotyped selected variants in an independent study population of incident HNSCC from The University of Texas M.D. Anderson Cancer Center (Houston, TX). Pairwise linkage patterns between the strongest associated variants in the discovery phase were used to prune and select variants for validation phase genotyping. After quality control, 31, 30, 28, and 19 variants were tested for replication of effect with overall HNSCC, oral cavity cancer, pharyngeal cancer, and laryngeal cancer, respectively. By pooling genotyping data across studies, we performed a joint analysis adjusting for age, sex, race, smoking status, and tumor stage, with a total of 2,046 cases (with complete covariate data) of HNSCC. In site-specific analyses 648, 1,081, and 317 cases were included in the oral cavity, pharyngeal, and laryngeal analyses, respectively. Summary statistics for all variants tested for replication are presented in Supplementary Table S2. Two of these variants, rs1816158 (oral cavity cancer survival) and rs56161233 (laryngeal cancer survival), had P < 1.0E−03 in the joint analysis, a lower joint analysis P value than discovery phase P value, and a consistent direction of effect (Table 2). Carriers of the minor allele for rs1816158 and rs56161233 had worse probability of survival than those homozygous for the major allele (Fig. 1A and B). rs1816158 is within intron 1 of the long noncoding RNA (lncRNA) MIR100HG was significantly associated with oral cancer survival (joint analysis oral cavity cancer; HR, 1.56; 95% CI, 1.21–2.00; Table 2). In addition, rs1816158 overlaps with promoter and enhancer histone marks, as well as DNase hypersensitivity sites, in several tissues (Supplementary Table S3). No nearby genotyped variants were in strong linkage with rs1816158 (Supplementary Fig. S3A). rs56161233 is located within the 3′UTR of SH3BP4 and was validated for its association with laryngeal cancer survival (joint analysis HR, 2.87; 95% CI, 1.73–4.75; Table 2). rs4589703 was the only genotyped variant in strong linkage (r² > 0.8) with rs56161233 and had a modest association with laryngeal cancer survival in the discovery phase (Supplementary Fig. S3B), which was dependent upon rs56161233 in conditional analyses (Supplementary Fig. S4). Neither rs56161233 or rs56161233 showed evidence of an association in the analysis of overall HNSCC (Supplementary Table S4), suggesting their effects are site specific. Although no variants from the overall HNSCC and pharyngeal cancer analyses met our replication criteria, several variants approached replication (Table 2), including RPL28 variants rs56312243 and rs45494801, both of which have been identified as eQTL loci for RPL28 in multiple tissues (Supplementary Table S5). Using the Schoenfeld residual test to assess the proportional hazards assumption, we observed no evidence for violation of this assumption in these analyses. Generally, similar results were obtained in a sensitivity analysis restricted to Caucasian subjects (Supplementary Table S6). In addition, we conducted a sensitivity analysis limiting follow-up to 5-year survival and observed consistent effect sizes (Supplementary Table S7).

As genetic variation at HNSCC-associated miR-SNPs may result in gene expression changes that modify cellular behavior and ultimately survival times, we tested the association between expression of genes containing miR-SNPs that met our replication criteria (rs1816158, MIR100HG; rs56161233, SH3BP4), and overall survival, in HNSCC subjects from the TCGA HNSCC project (Table 3). To determine the potential effects of MIR100HG-derived miRNAs on overall survival in oral cavity cancer, we utilized miRNA-seq data from available oral cavity cancer cases (n = 303). Expression of miR-100 was significantly associated with oral cancer survival (HR, 1.25; 95% CI, 1.06–1.49), in multivariable Cox proportional hazards regression models adjusting for age at diagnosis, sex, race, and tumor stage (Table 3). These findings suggest increased miR-100 expression may alter cellular behavior in such a way that modifies survival times. We did not observe a significant association between SH3BP4 expression and overall survival in laryngeal cancer cases (n = 108; HR, 1.13; 95% CI, 0.65–1.97; Table 3).

Genetic variation within miRNA target sites can create, destroy, or alter the binding affinity for complimentary miRNAs. To identify miRNA–mRNA interactions potentially disrupted by genetic variation, we intersected predicted miRNA target sites made using the miRanda algorithm with genomic coordinates of SNPs outlined in Table 2. Primary tumor miRNA expression data from the TCGA HNSCC project was also used to assess the potential functionality of implicated miRNAs. We identified multiple high-confidence miRNA–mRNA interactions overlapping each of rs16988668 (ZSCAN22), rs3831960 (ZSWIM5), rs56312243 (RPL28), rs77506493 (ZNF766), and rs56161233 (SH3BP4; Table 4). Several of the identified miRNAs are highly expressed in normal or tumor tissue from HNSCC subjects (Table 4). The highest confidence miRNA target site observed in this analysis was the interaction between miR-2110 and SH3BP4 overlapping laryngeal cancer–associated miR-SNP rs56161233, which was predicted with a mirSVR score in among the top 1% of all mirSVR predictions genome wide. In addition, miR-548m, miR-593-3p, and miR-3150a-3p are also predicted to target SH3BP4 within the region of the 3′UTR containing rs56161233 (all among the top quartile of mirSVR predictions). Multiple miRNA target sites overlapped each of miR-SNPs, rs3831960 and rs16988668, that approached replication in the association analysis of overall HNSCC survival (Table 4). miR-96-3p was predicted to target ZSWIM5 with notably high confidence (mirSVR score 6th percentile), while ZSCAN22 was predicted to be targeted by miR-339-5p, which is abundantly expressed in both normal and tumor tissue (90th and 92nd percentiles, respectively) from HNSCC subjects. Predicted miRNA target sites were also identified within RPL28 and ZNF766 that overlap rs56312243 and rs77506493, respectively.

We next sought to identify miRNA target sites whose binding is disrupted by genetic variation at these loci. We mined three publicly available databases (PolymiRTSv3.0, miRNASNP, and MirSNP) that contain predicted effects of genetic variation at miRNA-related loci. As above, for each miRNA–mRNA relationship identified, miRNA expression was calculated using normal and tumor tissue samples from TCGA HNSCC subjects. Details regarding the predicted effects of genetic variation upon the miRNA–mRNA interactions from each database are provided in Supplementary Tables S8–S10. TargetScan context+ score differences between each allele of a specified variant are provided in the PolymiRTS database and provide a likelihood measure that the variant will affect miRNA binding. miRNA–mRNA interactions with context+ score differences in at least the top 20% of all available scores were identified for rs4127682 (RGS3), rs16988668 (ZSCAN22), rs56312243 (RPL28), rs77506493 (ZNF766), and rs56161233 (SH3BP4; Supplementary Table S8). Several miRNAs whose binding is predicted to be disrupted by genetic variation at these loci are expressed at high levels in normal and tumor tissue from TCGA HNSCC subjects. Contrastingly, predictions stored in the miRNASNP database are made using the TargetScan and miRanda prediction algorithms simultaneously to identify miRNA target sites created or destroyed by genetic variation. From the miRNASNP database, we identified target sites created or destroyed by genetic variation at rs4127682 (RGS3), rs16988668 (ZSCAN22), rs75820821 (CCDC97), rs56312243 (RPL28), and rs56161233 (SH3BP4; Supplementary Table S9). Notably, six miRNA target sites are predicted to be created by presence of the T allele for rs56312243 (RPL28), while six target sites are predicted to be destroyed by the T allele. miR-486-3p, one of the target sites predicted to be created by presence of the T allele of rs56312243, was expressed in the 83rd and 73rd percentiles of normal and tumor tissue samples, respectively. Finally, for predictions in the MirSNP database, the miRanda algorithm was used to identify miRNA-binding sites created, destroyed, enhanced, or decreased by genetic variation. Predicted effects for several of the target sites predicted to be created or destroyed by rs4127682 (RGS3), rs16988668 (ZSCAN22), rs75820821 (CCDC97), rs56312243 (RPL28), rs45494801 (RPL28), and rs56161233 (SH3BP4) were similar between the miRNASNP and MirSNP databases (Supplementary Tables S9 and S10). Presence of the G allele for rs4127682 was predicted to enhance and decrease binding of miR-1909-3p and miR-4763-3p to RGS3, respectively, across five RGS3 transcript variants (Supplementary Table S10). Many of these miRNAs were expressed at high levels in normal and tumor tissue from TCGA HNSCC subjects. miR-2355-5p, predicted to be disrupted by genetic variation of rs45494801, was among the most abundantly expressed miRNAs, in the 89th percentile in normal tissue, and 93rd percentile in tumor tissue. These results highlight miRNA–mRNA interactions with potential to be disrupted by the miR-SNPs tested for association with overall survival of HNSCC.

Studies of the association between germline genetic variation and clinical outcomes in cancer can reveal insights into the complex regulatory networks related with disease progression, as well as identify biomarkers that may possess clinical utility. Previous studies of miR-SNPs and overall survival in HNSCC (21, 22, 24) have used candidate-based approaches. We present a genome-scale evaluation of the contribution of common miR-SNPs to survival in overall HNSCC, and site-specific disease, validating associations of MIR100HG variant rs1816158 with overall survival in oral cavity cancer, and SH3BP4 3′UTR variant rs56161233 with overall survival of laryngeal cancer. Moreover, providing mechanistic insight into the association of rs1816158 with oral cavity cancer, we demonstrated that expression of MIR100HG-derived miR-100 is associated with head and neck cancer survival in TCGA tumors. We also identify multiple high-confidence miRNA target sites predicted to bind the 3′UTR of SH3BP4 in the region containing rs56161233. In addition, we describe miRNA target sites located in the 3′UTR of SH3BP4 that are likely to be disrupted by genetic variation at rs56161233. Finally, we identify several additional miR-SNPs with putative survival associations for head and neck cancer where there is also high confidence for miRNA target sites and disruption of miRNA–mRNA interaction introduced by the miR-SNP: oral cavity cancer, rs3831960 (ZSWIM5, 3′UTR), rs16988668 (ZSCAN22, 3′UTR); laryngeal cancer, rs3777710 (FAM8A1, 3′UTR), rs12280753 (BUD13, 3′UTR). Together, these findings demonstrate miR-SNPs are associated with overall survival in HNSCC and provide evidence that genetic variation at these loci result in functional changes to gene regulation.

MIR100HG is a lncRNA from which 3 miRNAs (miR-100, let-7a-2, miR-125b-1) and one coding gene (BLID) are derived. Multiple studies suggest that MIR100HG and its derived genes contribute to various cancers. MIR100HG has been identified as an oncogene in acute megakaryoblastic leukemia (39) and implicated in regulation of lymph node metastases in early-stage cervical cancer (40). mir-100 induces epithelial–mesenchymal transition in mammary epithelial cells, while suppressing tumorigenesis, migration, and invasion (41). Further evidence of potential tumor-suppressive effects exists in the context of adult stem cells (42), although the effects of MIR100HG and its derived genes in cancer are likely highly complex and tissue dependent, warranting further study. Most recently, miR-100 and miR-125b were found to be overexpressed in HNSCC cell lines, and to mediate resistance to cetuximab (43). Given these data, we would expect the C allele of rs1816158, associated with increased risk of death in our study (Table 2), to increase miR-100 expression; however, this will require validation in future studies. Let-7a-2, one of the three isoforms that produce mature let-7a, is a member of the highly conserved let-7 miRNA family. Let-7a expression is downregulated in several cancers and negatively regulates oncogenes RAS (44), HMGA2, (45), and MYC (46). Expression of let-7a and its family members are tumor suppressive across a range of cancer types (47). Although let-7a-2 expression was not associated with overall survival of HNSCC in this study, given the known functions of let-7a, let-7a-2 cannot be discounted as a potential mediator of the observed association between rs1816158 overall survival in HNSCC. Overall, our findings add to the growing body of data supporting a complex role for MIR100HG in tumor biology.

SH3BP4 was among the first proteins identified to mediate cargo-specific control of clathrin-mediated endocytosis (48). Further studies have indicated that SH3BP4 negatively regulates amino acid-Rag GTPase-mTORC1 signaling (49). As hyperactivation of mTORC1 has been observed in several cancer types (50), it has been postulated SH3BP4 exerts tumor-suppressive effects over mTORC1 signaling (49). SH3BP4 has also been identified as a required component of a molecular complex, including FGFR2b and PI3K that controls receptor recycling and cellular outcomes in response to FGF10 stimulation (51). The diverse functions of SH3BP4 provide a rationale for how miRNA-related genetic variation in SH3BP4 may modify survival in HNSCC. Although we did not observe an association between laryngeal cancer survival and SH3BP4 expression, this does not preclude the possibility that rs56161233 affects survival through modification of miRNA-mediated gene regulation. miRNA binding can result in decreased mRNA translation rather than promoting enzymatic cleavage of the target (20). Future studies involving proteomic data collection and analysis will be required to delineate such instances. Furthermore, in recent years, development of the competing endogenous RNA (ceRNA) hypothesis has posited that competition for miRNA binding between transcripts is a substantial determinant of gene regulation (52). Disruption of miRNA binding through genetic variation within a miRNA target site may therefore alter regulation of a gene that shares this target site, while potentially leaving expression of the gene containing the variant mostly unchanged. Such phenomena have not been largely explored and deserve study.

Recent technological advances have facilitated the identification of predicted and experimentally validated miRNA target sites throughout the genome. Here, we leverage a novel genotyping platform to assay miRNA-related genetic variation at the genome-scale level. Our study represents a substantial improvement upon the existing candidate-based approaches. Furthermore, our study leverages detailed follow-up information on subject-specific survival to identify miR-SNPs associated with overall survival of HNSCC in two population-based studies. However, our study is limited by somewhat modest sample size for the detection of germline genetic variants associated with outcomes of a complex trait such as overall survival of HNSCC. Collection of detailed follow-up data in additional population-based studies of HNSCC is needed to increase statistical power in future work, particularly for tumor site-specific analyses. In addition, detailed exposure information regarding HPV infection and alcohol consumption was not available in the validation population, preventing adjustment for these factors in joint population analysis. HPV infection is an important predictor of overall survival in oropharyngeal cancers and alcohol consumption is a prognostic factor for overall survival in HNSCC (53). Lack of available data likely contributed limited replication of miR-SNPs with putative survival associations observed in the discovery population. Furthermore, use of HPV infection as a prognostic factor in HNSCC relates to tumor cell HPV infection; therefore, potential discordance between the HPV seropositivity detection methods used here and tumor HPV status could introduce bias (54). Although the presence of HPV-specific antibodies has been associated with HPV infection in tumor cells and improved prognosis in multiple studies (55, 56), future studies should take this into consideration in both study design and interpretation of results. Cross-study differences in variable structure may have also contributed to lack of replication. In addition, continuous tobacco exposure was available for discovery phase subjects, although only categorical exposure data (current, former, never smokers) were available in validation phase subjects. Our results highlight the importance of adjustment for clinical and demographic factors that strongly predict survival in genetic association studies. Although treatment data were also not considered in these analyses, future studies capable of collecting such data should aim to evaluate the potential effects of different treatment regimens on the hazards conferred by miR-SNPs. Finally, with continued development of miRNA prediction algorithms in recent years, and generation of new datasets of experimentally validated miRNA target sites, more up-to-date lists of miRNA target sites may be available for future studies aiming to evaluate genome-wide miRNA-related genetic variation. Detailed collection of specific exposures in future studies will improve statistical power for the detection of genetic variants associated with survival in diseases with major environmental contributors, such as HNSCC. Overall, our genome-scale analysis of miRNA-related genetic variation with survival in HNSCC highlights the power of using functional annotation to guide genetic association studies, as well as the need for more powerful evaluation of the association between miR-SNPs and survival in HNSCC.

J. Gui is a consultant/advisory board member for White River Junction VA Medical Center and Celdara Medical LLC. No potential conflicts of interest were disclosed by the other authors.

Conception and design: O.M. Wilkins, G. Li, K.T. Kelsey, B.C. Christensen

Development of methodology: O.M. Wilkins, L.A. Salas, E.M. Sturgis, G. Li, B.C. Christensen

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): G. Li, K.T. Kelsey, B.C. Christensen

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): O.M. Wilkins, A.J. Titus, L.A. Salas, J. Gui, M. Eliot, E.M. Sturgis, G. Li, K.T. Kelsey, B.C. Christensen

Writing, review, and/or revision of the manuscript: O.M. Wilkins, A.J. Titus, L.A. Salas, R.A. Butler, E.M. Sturgis, K.T. Kelsey, B.C. Christensen

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): O.M. Wilkins, R.A. Butler, K.T. Kelsey

Study supervision: E.M. Sturgis, K.T. Kelsey, B.C. Christensen

This work was supported by NIH grants R01DE022772 and R01CA216265 to B.C. Christensen and R01CA078609 to K.T. Kelsey. We acknowledge funding contributions from The University of Texas M.D. Anderson Christopher and Susan Damico Chair in Viral Associated Malignancies, National Institute of Environmental Health Sciences grants R01ES11740 and R01CA131274, and NIH grant P30CA016672 (to The University of Texas M.D. Anderson Cancer Center). A.J. Titus was supported in part by T32LM012204. We would like to thank Dr. Qingyi Wei for agreeing to collaborate on R01DE022772 before he moved on from The University of Texas M.D. Anderson Cancer Center. We are grateful to Kevin C. Johnson for discussions that improved this manuscript.

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.