Genetic variants in CYP2B6 and HSD17B12 associated with risk of squamous cell carcinoma of the head and neck.

As shown in the study workflow depicted in Figure 1 , this is a two-phase study design which included a discovery (MDACC) dataset and then a replication (OncoArray) dataset. Table 1 presents the distribution of characteristics in both the discovery and replication datasets. The mean age of cases and controls are 57.9 vs 50 years old in the discovery dataset and 59.7 vs 58.1 years old in the replication dataset, and there were more males in cases than in controls in both the discovery dataset(77.2% vs 55.1%) and replication dataset (74.2% vs 70.9%). There were 2,171 case and 2,169 controls with smoking and drinking data available in the discovery dataset. More smokers and drinkers were found in cases than in controls (68.4% vs. 47.0%, P < 0.0001; 73.5% vs. 54.6%, P < 0.0001, respectively). Of the cases in the discovery population, there were 631 (29.1%), 1,144 (52.7%), 316 (14.6%), and 78 (3.6%) patients with cancers of the oral cavity, oropharynx, laryngeal, and hypopharynx or overlapping cancer, respectively. In the replication dataset, there were 2,568 (49.5%), 2,328 (44.8%), 295 (5.7%) patients with oral cavity cancer, oropharyngeal cancer, and hypo-pharyngeal or overlapping cancer, respectively. Two and 14 patients in the discovery and replication datasets, respectively, were missing values for histological types. Single-locus analysis was performed for each of the SNPs with an imputation quality r 2 ≥ 0.3 and a minor allele frequency ≥ 0.01 in the MDACC discovery population. There were 212, 263, and 173 SNPs with P < 0.050 in the main effect analysis for SCCHN, and stratified analyses by oral and oropharyngeal, respectively ( Supplementary Table 2 - 4 ). After multiple test correction, there are 110, 145, and 110 SNPs passed BFDP threshold ≤ 0.8 in both the main effect analysis and the stratified analyses by organ site (i.e., oral, and oropharyngeal), respectively ( Supplementary Tables 2 - 4 ). No SNPs could pass the FDR correction (≤ 0.2). The Manhattan plots of these results are shown in Figure 2A - C . SNPs that passed BFDP correction were then selected for replication with the OncoArray dataset, with 2, 2, and 1 SNPs replicated with P 0.1; I-squared < 50%), a meta-analysis of discovery and replication results was performed with the fixed-effects model ( Table 2 ). In this model, 2 SNPs (rs145533669 and rs35246205) in CYP2B6 within the 19p13.2 loci exhibited a significant association with SCCHN risk. The variant T allele of the leading SNP CYP2B6 rs145533669 was associated with decreased SCCHN risk (OR = 0.60, 95% CI = 0.46–0.78, P = 1.57 × 10 −4 ). Two SNPs ( EPHX1 rs117522494 within the 8p21.2-p21.1 loci and CYP2B6 rs145533669 within the 19p13.2 loci) showed significant associations with oropharyngeal cancer risk, with the variant T alleles of SNPs rs117522494 and rs145533669 associated with increased (OR = 1.45, 95% CI: 1.16–1.82, P = 0.001) and decreased risk (OR = 0.62, 95% CI: 0.44–0.85, P = 0.004), respectively. In addition, one SNP (rs4359199) in HSD17B12 exhibited significant associations with the risk of oral cancer, with the variant allele C associated with increases cancer risk (OR = 1.11, 95% CI: 1.04–1.18, P = 0.0018). Considering only few SNPs identified in the above analyses, we calculated the observed power by using the online Genetic Association Study (GAS) Power Calculator ( https://csg.sph.umich.edu/abecasis/gas_power_calculator/ ) to exclude the possibility of false negative results. As shown in Supplementary Table 5 , the power in the oropharyngeal cancer group with 1,144 cases and 4,493 controls appears to be enough to detect the effects of SNPs with a minor allele frequency > 0.02 and a relative risk > 1.68. Otherwise, the power would be insufficient. For the analysis of the oral cancer group (with 631 cases and 4493 controls), the power was not adequate to exclude the type II error. Regional association plots of the four SNPs were presented in Supplementary Figure 1A - D . Linkage disequilibrium (LD) analysis was performed for these four identified SNPs ( Supplementary Figure 2 ). There is no LD observed for the two SNPs (rs145533669 and rs35246205) in CYP2B6. Two other SNPs (rs117522494 in EPHX2 and rs4359199 in HSD17B12 ) were located at different chromosome regions (8p21.2 and 11p11.2). As smoking and drinking are known risk factors for SCCHN, we also performed interaction analysis of the four identified SNPs with smoking and alcohol drinking. As shown in Supplementary Table 6 , we found a marginal significance was found for the interaction effect between SNP rs117522494 in gene EPHX2 and drinking status in oral cancer (interaction P = 0.080), and SNP rs4359199 in gene HSD17B12 region exhibited significant interaction effects with drinking status in overall SCCHN (interaction P = 0.028) and oropharyngeal cancer (interaction P = 0.033). We then performed stratified analysis by smoking and drinking status for this SNP. As shown in Table 3 , the variant C allele of SNP rs4359199 was associated with increased risk of SCCHN (OR = 1.10, 95%CI: 0.99–1.23, P = 0.085) and for oral cancer (OR = 1.22, 95%CI: 1.03–1.43, P = 0.019) in smokers. While a protective effect was observed for the T allele on oropharyngeal cancer risk in non-drinkers (OR = 0.82, 95%CI: 0.68–0.99, P = 0.043). No significant effects were observed in smokers or non-smokers for SCCHN or after stratification by oropharyngeal and oral cancers. To test the multiple-markers effect, gene-based analysis was performed with MAGMA, which is based on a multiple linear principal components regression model using an F-test to compute the gene p-value 29 . The results showed that the associations between ten metabolism-related genes (i.e., AKR1A1, EPHX1, SULT1B1, UGT3A1, NAT2, AKR1C1, HSD17B12, SULT1A2, SULT1A1, and CYP2A6 ) and SCCHN risk were statistically significant in the MDACC discovery dataset (P < 0.05, Supplementary Table 7 ). Two of the associations (for AKR1A1 and HSD17B12 ) were replicated in the OncoArray dataset ( P = 0.053 and 0.003, respectively). For those replicated SNPs, the association results from eQTL analysis in multiple tissues were shown in Figure 3 : A) rs35246205 and CYP2B6 , C) rs117522494 and EPHX2 , E) rs4359199 and HSD17B12 , respectively. Multiple pairs were found with significant correlations in different tissues. As the PAHs and TSNA were mainly metabolized in liver and immunity capacity can influence the risk of SCCHN, we were more interested in the eQTL results in these related tissues. Significant correlations were found between rs4359199 and the mRNA expression of HSD17B12 in blood cells ( P = 4.9E-21), monocytes ( P = 3.3E-9), and T cells ( P = 1.74E-08) ( Figure 3E ). sQTL associations in multiple tissues were presented in Figure 3 : B) rs35246205 and CYP2B6 ; D) of rs117522494 and EPHX2 ; F) rs4359199 and HSD17B12 , respectively. Notable significant correlations were found between rs35246205 and CYP2B6 in liver ( P = 8.71E-04 and 0.011 in Figure 3B ); rs4359199 and HSD17B12 in LCL, T-cells, and liver cells ( P = 1.0E-19, 1.1E-04, and 1.3E-04, respectively, in Figure 3F ). PheWAS results (with gene level P < 0.05) of the three identified genes and 600 traits in UK Biobank were presented in Figure 4A - C . As it shown, significant correlations were found between these genes and multiple traits. For example, there were significant association between CYP2B6 and cigarettes used per day ( P = 2.27E-11), squamous cell lung cancer risk ( P = 0.004) ( Figure 4A ). EPHX2 was also found to be correlated with traits including ever smoker ( P = 2.73E-05), breast cancer risk ( P = 0.012), Cadherin-15 contents in blood ( P = 4.67E-06) ( Figure 4B ). HSD17B12 was found to be associated with traits in the metabolic domain, including body Mass Index ( P = 2.95E-28), type 2 Diabetes ( P = 1.21E-10), number of cigarettes previously smoked daily ( P = 4.35E-04), and lung cancer risk ( P = 0.008) ( Figure 4C ). Further functional annotation using HaploReg V4.1 also revealed that CYP2B7P rs117522494 is located in the potential promoter region with some evidence from DNAse, histone, and protein bound ChIP-Seq experiments ( Supplementary Table 8 ). It should be noted that rs4359199 is located in the intron of gene HSD17B12 but was correlated with 50 eQTL traits in GTEx.

The discovery dataset included 2,171 SCCHN cases ascertained at the Head and Neck Surgery Clinic at the University of Texas MD Anderson Cancer Center (MDACC; Houston, TX) between December 1996 and July 2011 11 , 21 , 22 . All cases were individuals with newly diagnosed, histologically confirmed, previously untreated squamous cell carcinoma of the oral cavity, oropharynx, hypopharynx, or larynx. Blood samples were collected for genomic DNA extraction and genotyping was performed with Illumina HumanOmniExpress-12v1 BeadChip 11 . This study included 4,493 cancer-free controls, who were recruited from genetically unrelated visitors who accompanied cancer patients to MDACC outpatient clinics, or from the MDACC melanoma study with GWAS data deposited in database of Genotypes and Phenotypes (with dbGaP accession#: phs000187.v1.p1) 21 – 23 , or from the Study of Addiction, Genetics and Environment (SAGE) (SAGE; dbGaP accession #: phs000092.v1.p1) 24 . Of these controls, there were 1,149 cancer-free individuals recruited for the SCCHN study with genomic DNA genotyped by using Illumina HumanOmniExpress-12v1 BeadChip; 1,022 cancer-free individuals from the melanoma GWA study whose genomic DNA was genotyped by Illumina Omni1-Quad_v1–0_B BeadChip, and 2,322 cancer-free individuals of European descendent from the SAGE study 24 whose genotyping data were generated by Illumina Human1Mv1 BeadChip. All of these three GWAS datasets contained around 1 million genotyped SNPs. The SCCHN GWAS data is available in dbGaP (accession #: phs001173.v1.p1) 11 . The replication dataset was from the OncoArray study of oral and pharynx cancer, which is part of the International Head and Neck Cancer Epidemiology Consortium (INHANCE) and included 6,034 cases and 6,585 controls derived from 12 epidemiological studies in European populations 9 . Genomic DNA extracted from blood or buccal cells was genotyped at the Center for Inherited Disease Research (CIDR) with the Illumina OncoArray custom array. Related GWAS data were requested from dbGaP (accession#: phs001202.v1.p1) in which genotyping data were available for 6,034 cases and 4,062 controls. After quality control 11 , 5,205 cases and 3,232 controls of European ancestry were included in the present study. Based on previous publications, we selected 43 genes encoding enzymes that play a key role in the metabolism of either PAHs or TSNAs 7 , 14 , 25 . These genes are listed in Supplementary Table 1 . Imputation was performed on the Michigan imputation server ( https://imputationserver.sph.umich.edu ) with the Haplotype Reference Consortium (HRC) reference panel (Version r1.1 2016) consisting of 64,940 haplotypes for individuals of predominantly European ancestry. In the discovery stage, the genotyping/imputed data was extracted for 4,177 SNPs (555 genotyped and 3,622 imputed) with a minor allele frequency ≥ 0.01 and r-square ≥0.3 within 5 kb upstream and 5 kb downstream of the 43 metabolizing enzyme genes from the SCCHN GWAS of MDACC. Single-locus analysis with an unconditional logistic regression model was performed to estimate odds ratios (ORs) and 95% confidence intervals (CIs) per effect allele using the PLINK (v2.00) software with adjustment for the age, sex and top five significant principal components (PCs). Stratified analysis by tumor site (i.e., oral, oropharyngeal, and hypo-pharyngeal and laryngeal cancer) was performed using the false discovery rate (FDR) followed by Bayesian false-discovery probability (BFDP) for multiple-testing correction (with the setting of prior probability = 0.1 and the upper bounder of detectable OR = 3). Because the SNPs under investigation were mostly from imputation (3,622 out of 4,177 SNPs), those SNPs with BFDP ≤ 0.8 were then examined in the OncoArray study. The summary results of both discovery and replication datasets were then examined by a meta-analyses with the inverse variance-weighted average method. A random-effects model was applied as a Q-test P ≤ 0.10 or I 2 >50.0%; otherwise, a fixed-effects model was applied. eQTL and sQTL analyzes were applied to elucidate biological mechanisms of those identified genetic variants by using FIVEx ( https://fivex.sph.umich.edu/ ). Phenome-wide association analysis (PheWAS) analysis was performed to investigate the association between identified genes and 600 traits from UK Biobank release 2 by using the online tool: https://atlas.ctglab.nl/PheWAS . The related methods can be found elsewhere 26 . Briefly, the MAGMA gene analysis was performed to test the gene-trait correlation. Firstly, K SNPs are assigned to genes with 1-kb window both sides. SNP-wise mean models were then used to combine SNP Z statistics into a gene test-statistics T : T = ∑ j K Z j 2 , where Z J = Φ ( p j ), Φ is the cumulative normal distribution function; p j is the marginal P value of SNP j . Gene based p-value: P = P r ( T ≥ T obs ). In the evaluation of gene test-statistics, LD between SNPs in the gene is also be estimated to produce accurate results. We also used LocusZoom 27 and Haploview v4.2 28 to construct the regional association plots and linkage disequilibrium (LD) plots, respectively.

In this study, we performed a comprehensive candidate-pathway analysis of SNPs in 43 PAH- and TSNA-metabolizing genes and the risk of SCCHN by using a two-phase study design. We identified four SNPs in three gene regions (i.e., CYP2B6, EPHX2, and HSD17B12 ) associated with risk of SCCHN and its subtypes. We also revealed significant interaction effects between SNP rs4359199 in HSD17B12 and drinking status. In-silico eQTL and sQTL analyses also provided functional evidence for the identified SNPs (rs35246205 and rs4359199) in CYP2B6 and HSD17B12, respectively . Our results suggest that the variant alleles of the these identified SNPs might be associated with increased cancer risk of SCCHN and its subtypes by regulating corresponding gene expression or affecting mRNA alternative splicing. Cytochrome P450 (CYP) enzymes are crucial for the metabolic activation of PAHs and TSNA, which are important carcinogens in cancer development. Dysfunctions of the relevant genes have been found to the increased risks of multiple cancers 25 . In this study, we found SNPs located in the CYP2B6 region were associated with the risk of SCCHN and oropharyngeal cancer. Our eQTL and sQTL results showed that the identified genetic variant rs35246205 might be a tagging SNP of other functional variants influencing the alternative splicing of CYP2B6 in liver. Polymorphisms in this gene have also been associated with the risk of breast cancer and acute myeloid leukemia 30 , 31 . CYP2B6 and other family members, such as CYP2A6, are mainly responsible for converting of BaP and NNK to their intermediates that produce DNA adducts 32 . Although inconsistent results were also reported 33 , several studies reported a potential link between this gene and nicotine metabolism 34 – 36 , which were consistent with the PheWAS results and suggested that CYP2B6 might also contribute to squamous cell cancer risk by changing smoking behavior. Two SNPs in two other gene ( EPHX2 and HSD17B12) were also associated with the risk of oropharyngeal cancer and hypo-pharyngeal and laryngeal cancer, respectively, in the present study. EPHX2 is a member of the epoxide hydrolase family whose encoded enzyme can bind to specific epoxides (the intermediates from PAH) and convert them to the corresponding dihydrodiols 37 . One study has reported that the expression of EPHX2 was higher in an oropharyngeal HPV-positive cohort than in the negative cohort from The Cancer Genome Atlas (TCGA) 38 . Although previous reports suggested that EPHX2 might be involved in the regulation of smoking induced inflammation and autophagy, the underlying function of EPHX2 involved in oropharyngeal is still unknown 39 , 40 . Smoking and drinking are two important independent risk factor for the development of oral and pharyngeal cancers. In the present study, we found a significant interaction effect between HSD17B12 SNP rs4359199 and drinking status in oropharyngeal cancer. The variant allele C presented a protective effect on oropharyngeal cancer risk in non-drinkers, but increased the risk of oral cancer in drinkers. However, few reported studies have investigated the interaction effect of this genes and the underlying mechanism is still unclear. HSD17B12 play an important role in the regulation of sex steroid metabolism by catalyzing the conversion between 17-keto and 17-hydroxysteroids. The dysregulation of this gene has been implicated in multiple estrogen and androgen-related diseases including breast cancer, endometriosis, prostate and non-small cell lung cancer 41 , 42 . Polymorphisms in this gene have been linked to risk for multiple cancers including gastric cancer, endometrial cancer, breast cancer and prostate cancer 43 – 46 . PheWAS results showed that the identified SNP involves multiple metabolic-related phenotypes. This is consistent with the current findings that HSD17B12 may play its role in inflammation and cancer development by upregulating fatty acid synthesis in human cancers. 47 While the present study is the first to report on associations of genetic variants in these two gene and risk for SCCHN and its sub-types, further replications are required. There are several limitations to this study. Firstly, the underlying biological mechanism of the identified associations are still unclear. Although we have provided some in-silico functional evidence, most of them were from whole blood or other tissues but not the target tissue; additional evidence in target tissues or cell lines will be required. Second, the smoking and alcohol data are not available in the OncoArray study. The identified interaction effects in the MDACC study need to be further replicated in another independent study. Thirdly, HPV infection is one of the primary risk factors in oropharyngeal cancer. However, we did not have adequate HPV data available to detect the possible interactions between SNPs and HPV status. In addition, SNPs with minor effects or low frequency cannot be identified due to sample size limitations of the MDACC and OncoArray studies. In conclusion, we found that two SCCHN risk-associated SNPs (rs35246205 in CYP2B6 and rs4359199 in HSD17B12 ) were correlated with mRNA expression levels of the corresponding genes and that two SNPs, rs117522494 in CYP2B6 and rs4359199 in HSD17B12 , are located at the potential promoter regions and may regulate mRNA expression and alternative splicing. Further functional validation and population replication are warranted to substantiate the present findings.

Squamous cell carcinoma of the head and neck (SCCHN) starts from the mucosal epithelium in the oral cavity, pharynx and larynx. It is the sixth most common malignancy and seventh leading cause of cancer-related deaths worldwide 1 , 2 . In the United States, there are approximately 65,630 new SCCHN cases and 14,500 SCCHN-related deaths in 2020 3 . Epidemiology studies have revealed some important risk factors of SCCHN, such as exposure to tobacco and alcohol as well as infection with human papillomavirus (HPV) 4 – 6 . However, only a fraction of tobacco users, alcohol drinkers or individuals who contracted HPV develop SCCHN, suggesting an important role of genetic susceptibility in its etiology 7 , 8 . Two previous studies have reported that polymorphisms in alcohol-related genes including alcohol-dehydrogenase 1B ( ADH1B ) and ADH7 were associated with risk of head and neck cancer. Recent genome-wide association studies (GWASs) have also revealed multiple loci to be associated with risk of SCCHN and its subtypes (i.e., 2p23.3, 5p15.33, 5q14.3, 6q16.1, 6p21.33, 6p21.32, 9p15.3, 9q34.12, 10q26.13, 11p15.4, 11q12.2, 12q24.21 and 16p13.2) 9 – 11 . Polycyclic aromatic hydrocarbons (PAHs) and tobacco-specific nitrosamines (TSNA) are two of the most carcinogenic components of tobacco smoke. Accumulating evidence has demonstrated that variability in metabolic enzymes important in the bioactivation or inactivation of PAHs and TSNA is linked to risk of many cancer types 12 – 14 . Recent GWASs with large sample sizes also reported associations between genetic variants in metabolism-related genes and risk for multiple cancers [i.e., lung cancer 15 , 16 , bladder cancer 17 , colorectal cancer 18 ], but these hypothesis-free GWASs have not reported any SNPs in some well-known metabolic genes that were identified for head and neck cancer in previous molecular epidemiology studies using a candidate gene approach 7 , 19 , 20 . In the present study, we used a hypothesis-driven approach to comprehensively assess associations of genetic variants in 43 PAH and TSNA metabolism-related genes and the risk of SCCHN as well as its subtypes (e.g., oral and oropharyngeal cancers) in European populations using two published head and neck GWAS datasets, and explored interaction effects between identified SNPs and smoking/drinking status, as well as possible functional annotation of the identified SNPs.