A comprehensive gene-environment interaction analysis in Ovarian Cancer using genome-wide significant common variants.

The OCAC is an international multidisciplinary consortium formed in 2005 ( http://apps.ccge.medschl.cam.ac.uk/consortia/ocac/ ) with a goal of sharing data from worldwide ovarian cancer studies to establish reliable estimation of association between environmental and genetic factors related to risk of ovarian cancer 23 , 37 . Cases were defined as women with ovarian carcinoma (i.e., invasive epithelial ovarian cancers), fallopian tube cancer and primary peritoneal cancer. Controls were women without ovarian cancer and who had at least one ovary. For both cases and controls, individuals with prior cancers except non-melanoma skin cancers were excluded. In total, 28 single nucleotide polymorphisms (SNPs) previously identified through GWAS were included from 75 OCAC sites ( Table 1 ). The first 26 SNPs were found to be significantly associated with either ovarian cancer overall or one or more histotypes 36 . In addition, rs13255292 and rs10962643 were included because they were in the same region as two other significant SNPs but showed a strong independent association with ovarian cancer risk. The SNP at locus 15q26 (rs8037137), which was found to be genome-wide significant 13 , was not included because not enough non-carriers were present in our analytic dataset for examining interactions. The genetic data included both genotyped and imputed variants (imputation being carried out using phase 2 Hapmap reference panel). More details regarding genotyping and imputation of the genetic data have been previously described 12 , 17 , 18 , 20 . The methods for analyzing the SNP data in the OCAC have also been described previously 12 , 17 , 18 , 20 . Briefly, logistic regression models were fit to examine the association between ovarian cancer and each genetic variant under an additive model (using risk allele dosage). The models were adjusted for ethnicity, genotyping panel and the leading principal components for each ethnicity. The summary results are shown in Table 1 and are also available through the OCAC website (http:/apps.ccge.medschl.cam.ac.uk/consortia/ocac/). A total of seven established environmental risk factors for ovarian cancer were of primary interest ( Table 2 ), including four associated with decreased risk and three with increased risk for ovarian cancer or one specific histotype. These included: OCP use (measured as both ever/never and duration of OCP use (never users including <1 one year of use, 1-<5, 5+yr), tubal ligation (yes/no), breastfeeding (ever/never), parity (0, 1–2, 3+ full-term births (i.e., those lasting ≥6 months), type of menopausal hormone therapy use for more than 1 year after age 50 (never user, menopausal estrogen therapy only, any use of menopausal estrogen + progestin therapy), BMI (<25, 25-<30, 30+), and a history of endometriosis (yes/no). Four other environmental variables were included in our analysis, as covariates: baseline age (<50, 50-<55, 55-<60, 60-<65, 65–70, 70+ years), race (non-Hispanic white, Hispanic White, Black, Other), education (less than high school, high school graduate, some college, college graduate) and first-degree family history of ovarian cancer (yes/no). In addition to these four covariates, study site, OCP use, tubal ligation, parity, BMI and endometriosis were also included in all models for the environmental association analysis and gene by environment interaction analysis. A brief description of environmental data harmonization across OCAC study sites is provided in eMethod 1 in the Supplementary Material. To optimize power and enhance the chance for discovery, we carried out multiple imputation of the environmental data. The maximal amount of data was used for imputation (see eMethod 1 and eFigure 1 in the Supplementary Mate rial for details). A total of 19 studies comprising 13,722 cases and 22,975 controls with partially missing data were included for imputation. Of these 19 studies, 12 were from the US, 4 from Europe, 2 from Canada and 1 from Australia (see eTable 1 for a description of study sites). Further details for these 19 studies have been previously described (see Supplementary Material ). The environmental variables included in our analysis were multiply imputed by chained equations (MICE) to produce ten imputed datasets. See details of imputation model in eMethod 2.1 in the Supplementary Material . All analyses were performed on each of the ten imputed datasets, and coefficients/test statistics were properly combined to account for uncertainty due to imputation, following the recommended combination rule for multiply imputed datasets 38 (see details in eMethod 2.3 in the Supplementary Material ). Our marginal environmental association analysis was based on combined inference from the ten imputed versions of this harmonized E data. Logistic regression models were used for evaluating marginal associations between the environmental risk factors with ovarian cancer after adjusting for covariate. The estimated ORs, their 95% CIs, as well as two-sided Wald tests after accounting for imputation uncertainty are presented in Table 2 along with summary statistics of complete cases before imputation. Full results of the complete cases analysis using logistic regression models are presented in eTable 2 . After marginal analysis of the genetic and environmental risk factors, we considered gene by environment (G × E) interaction analysis both on the multiplicative (odds ratio/relative risk) and the additive (relative excess risk due to interaction/absolute risk) scale 39 . From the 19 studies with imputed environmental data, a subset of 17 case-control studies with 9,971 cases and 15,566 controls had available genetic data, thus G × E analyses were carried out on these 17 studies. Each imputed environmental dataset was merged with the genetic data for subsequent G × E analyses. Interaction analyses were then carried out separately on the ten imputed G × E datasets, and then all tests and coefficients reported were combined using appropriate multiple imputation combination rules 38 . For both multiplicative and additive interaction analysis, we started with global likelihood ratio tests (LRTs) for each G × E pair as several environmental factors had multiple categories resulting in tests for interactions with multiple degrees of freedom (df). These global joint tests, serving as a screening step for G × E interactions, were carried out for a total of 196 (7×28=196) G × E pairs. After the global tests, we then followed up on the suggestive interactions (with global test P-value < 0.2) and carried out a two-sided Wald test for interactions involving each separate category of an environmental risk factor. For the k- th SNP G k k = 1 , … , 28 , coded as a continuous allelic dosage, the j- th environmental risk factor E j j = 1 , … , 7 , and a set of confounders/covariates C q ( q = 1 , … , Q ) , the basic fitted model for the probability of ovarian cancer of the i- th subject, namely, π i , is of the following form: [M1] l o g i t ( π i | G k i , E j i , C 1 i , … , C Q i ) = β 0 + β G G k i + ∑ l = 1 L β E l I ( E j i = l ) + ∑ l = 1 L β G E l I ( E j i = l ) G k i + ∑ q = 1 Q ∑ m = 1 M q β C q m I ( C q i = m ) , where L = (levels of E j ) – 1, M q = (levels of C q ) – 1, and Q is the number of adjusted covariates. For testing the multiplicative interaction between G k and E j , we first used the global LRT with L degrees of freedom to test for the joint null hypothesis H 0 : β G E 1 = β G E 2 = … = β G E L = 0 . If the global test P-value < 0.2, we further assessed the multiplicative interaction at each level of E j by using a Wald test with one degree of freedom for the null hypothesis H 0 : β G E l = 0 for the l- th level. Due to limitations of existing software (CGEN) 40 for testing additive interactions with continuous dosage data, we used the maximal probable genotype for imputed SNPs. We further conducted the LRTs with binary collapsing of SNPs assuming a dominant genetic susceptibility model (given the constraints in software) 31 . For a given SNP G k and an environmental risk factor E j with L categories, a global LRT with L df was used for the following joint null hypothesis H 0 : { exp ( β E 1 ) + exp ( β G ) − 1 } exp ( β E 1 + β G ) = exp ( β G E 1 ) , … , { exp ( β E L ) + exp ( β G ) − 1 } exp ( β E L + β G ) = exp ( β G E L ) , where the regression coefficients ( β ) are log odds ratio parameters described in model [M1]. This null hypothesis is based on a rare disease assumption 41 , which is tenable for our study (lifetime risk of ovarian cancer in the US is approximately 1.3%) 42 . If the global LRT P-value < 0.2, we further assessed the additive interaction at each level of E j through the relative excess risk due to interaction (RERI) 41 . At the l -h level of E j , a Wald test with one degree of freedom (35) was used to test for the null hypothesis: H 0 : R E R I G E l = 0 ,  where  R E R I G E l = exp ( β E l + β G E l + β G ) − exp ( β E l ) − exp ( β G ) + 1. After the screening step, we further explored the structure of the most promising interactions (defined as global test P-value < 0.01). This was accomplished by exploring odds ratios corresponding to E in sub-groups defined by G (for the multiplicative interaction) or absolute risks for ovarian cancer in each configuration of the values of (G, E) (for the additive interaction). To better understand these two different scales of interaction, we also compared the observed joint ORs with the corresponding expected ORs under the multiplicative and the additive nulls. To estimate sub-group specific absolute risk (AR) for each stratum defined by a given SNP G k and environmental risk factor, we need the relative risk and the joint distribution of G k and E j . The former was estimated from the fitted model [M1], and the latter was empirically estimated from the observed joint frequency of E j and G k in the control population ( details in eMethod3 from the Supplementary Material ). Table 4 presents the bootstrap confidence intervals for the estimated ARs and the risk differences (RDs) (see details in eMethod4 in the Supplementary Material ). The results for G × E analysis are presented in Table 3 (multiplicative interaction), Table 4 (additive interaction) and eTable 5 (observed and expected joint OR under the two different nulls). All calculations were performed in the statistical software R 30 , 40 .

The marginal G analysis was carried out on 26,864 cases and 48,034 controls and the results are shown in Table 1 . These results are available through the OCAC website ( http://apps.ccge.medschl.cam.ac.uk/consortia/ocac/ ). A total of 36,697 women with 13,722 ovarian cancer cases from 19 sites were included in the marginal E analysis using the imputed datasets. All seven environmental risk factors were associated with ovarian cancer in the expected direction ( Table 2 ). OCP use for five or more years was associated with a 52% decrease in risk of ovarian cancer compared to never users (OR=0.48, 95%CI = 0.45 to 0.51). Tubal ligation (OR=0.73, 95%CI = 0.69 to 0.78) and breastfeeding (OR=0.76, 95%CI = 0.71 to 0.80) showed similar magnitudes of decreased risk. Also, having more than 3 children (versus none) was associated with a 50% (OR=0.5, 95%CI = 0.46 to 0.53) reduction in risk of ovarian cancer. Using menopausal estrogen therapy only for more than one year (OR=1.22, 95%CI = 1.12 to 1.34), being obese (OR=1.15, 95%CI = 1.08 to 1.22), and history of endometriosis (OR=1.60, 95%CI = 1.46 to 1.75) were all associated with increased risk of ovarian cancer. The inference remained robust before and after imputation ( eTable 2 .). The global LRT essentially serves as a screening approach to identify a list of potentially interesting interactions. All interactions with global LRT P-value < 0.2 (40 on multiplicative scale and 41 on additive scale) are listed in eTable 3 , while more detailed analysis of the top interactions, which showed the strongest significance (P-value < 0.01; 4 on multiplicative and 2 on additive scale), are shown in Table 3 and Table 4 , respectively. According to Global LRT results, the top interaction on the multiplicative scale was identified with the SNP rs13255292 and OCP use (ever and never use: P-value = 3.48 × 10 −4 ; duration of use [<1 yr, 1–5 yr, 5+ yr]: P-value = 7.26 × 10 −3 ) ( Table 3 ). None of the observed interactions were significant based on a Bonferroni threshold of 0.05/(28 × 7)= 2.55 × 10 −4 . For the most promising multiplicative interactions reported in Table 3 we carried out an in-depth analysis to better understand the structure of interactions by estimating the ORs (with accompanying Wald CIs and tests) corresponding to E in strata defined by G. For example, the OR for OCP use among women with the TT genotype for rs13255292 is estimated to be 0.53 (95%CI = 0.46 to 0.60), whereas for the CC genotype the estimated OR is 0.71 (95%CI = 0.66 to 0.77) suggesting a stronger protective effect of OCP use among TT genotypes ( Table 3 , Figure 1A ). When OCP use was further stratified by duration, we observed an interesting pattern in its interaction with rs13255292. The estimated OR corresponding to 1–5 year of OCP use vs < 1 year use in the TT genotype group was 0.58 (95%CI = 0.50 to 0.69) compared to an OR of 0.79 (95%CI = 0.72 to 0.87) among women with CC genotype, showing effect modification by the risk allele (C) of rs13255292 ( Table 3 , Figure 1B ). This is akin to the result with ever/never user. However, the OR corresponding to 5+ years of OCP use vs < 1 year of use for the TT genotype group was 0.43 (95%CI = 0.37 to 0.50) and for the CC genotype was 0.53 (95%CI = 0.49 to 0.58) ( Table 3 , Figure 1C ). With overlapping confidence intervals, there is no significant difference in the odds ratios for long-term OCP users across genotype sub-groups. Table 3 shows that the P-value of the Wald test for interaction of rs13255292 and 1–5 years of OCP use (vs < 1 yr) was lower (P-value = 4.74 × 10 −3 ), when compared to the P-value for interaction of the same variant with 5+ years of OCP use (vs < 1 yr) (P-value = 2.43 × 10 −2 ). For the most statistically significant additive interactions in Table 4 , we estimated the sub-group specific absolute risks (ARs) and risk differences (RDs) in each E by G stratum. For example, for the strongest additive interaction based on the global likelihood ratio tests in Table 4 , there was suggestive evidence that rs11658063 modified the effect of menopausal estrogen therapy use, compared to never use of menopausal hormone therapy (P-value = 3.01 × 10 −2 ). Among women with the GG genotype, never users of menopausal hormone therapy had an estimated AR of 1.33% (95%CI =1.26% to 1.40%) while women who used menopausal estrogen therapy had an estimated AR of 1.96% (95%CI = 1.59% to 2.33%), leading to an absolute risk increase of 0.63% (95%CI = 0.24% to 1.02%) ( Table 4 , eFigure 2 ). For women with the CC genotype, the estimated AR was 1.27% (95%CI = 1.23% to 1.32%) for never receiving menopausal hormone therapy and 1.36% (95%CI = 1.15% to 1.57%) for receiving menopausal estrogen only therapy. This implies virtually no increased risk from taking menopausal estrogen only therapy among women with the CC genotype (95%CI = −0.14% to 0.31%; Table 4 , eFigure 2 ). The results on the additive interactions were in general weaker in terms of the strength of P-values.

We have conducted a comprehensive multiplicative and additive interaction analysis of previously identified common genetic variants and environmental factors unequivocally associated with ovarian cancer risk. We observed six suggestive interactions (with P-value < 0.01), four on the multiplicative scale and two on the additive scale. The lack of statistical significance of interactions after multiple testing correction from a large collection of data and well-curated studies enable us to conclude that it is unlikely that there are substantive interactions with single variants and environmental factors regardless of the choice of scale. This is consistent with what has been observed for other cancers. One may argue that the Bonferroni threshold for multiple comparisons is likely to be conservative for this set of correlated environmental factors, but the general pattern of findings remains consistent with smaller magnitude of interaction effect sizes. However, there are several interesting findings from this analysis that may be worthwhile to follow-up in future G × E studies of ovarian cancer. In addition to guiding targeted prevention strategies, G × E analysis has the potential to provide mechanistic insight into the complex multifactorial structure of the underlying biological pathway. One issue complicating observed gene-environment interactions of even confirmed susceptibility loci is that the true casual alleles and the biological impact of the variants are unknown. Our top interaction is between OCP use and rs13255292. This variant lies in the 8q24 region which harbors several risk loci for ovarian cancer 18 and other cancers 43 , 44 . The SNP is in the PVT1 gene which interacts with the oncogene MYC 45 . MYC has long been reported to be at least in part under hormonal control 46 , 47 thus an interaction with OCP use is plausible. Conversely, our top additive interaction is between menopausal estrogen use and rs11658063 which falls in HNF1B . To our knowledge there is no relationship between HNF1B and hormones thus underscoring the difficulty of understanding these gene-environment interactions given our limited understanding of the function of the variants and even more broadly the biological role of the genes. The strongest interactions are observed with OCP use or menopausal estrogen use which are modifiable exposures. Our most promising finding is the potential interaction between SNP rs13255292 and OCP use. This finding, if replicated could potentially lead to improved understanding of exposure pathways. We present a comprehensive analytical framework to carry out post-GWAS G × E analysis on both multiplicative and additive scale. Our framework starting with data harmonization and imputation followed by Global likelihood ratio tests and single df Wald tests provides a principled analytic architecture for such analysis. Our analysis reiterates the well-known fact that testing the additive and multiplicative nulls are very similar when the marginal associations are weak but could depart when both marginal associations are large in magnitude and the sample size is finite. In eTable 5 , we present observed joint odds ratios for strata defined by G and E along with the expected odds ratios under the multiplicative null and the additive null. We use our top hit rs13255292 and OCP use (ever versus never) and length of OCP use (<1yr, 1-<5 yrs, 5+ yrs) as an illustration. One can note that the expected ORs are fairly close under both models. However, their estimated departure from the observed joint OR is more pronounced for the 1-<5 yrs sub-group when compared to 5+ yrs, explaining the suggestive evidence for rejecting the null. We discussed the multiplicative interaction results for rs13255292 and OCP use in the previous section. We now explore the structure of additive interaction for this G × E result ( Figure 2A-2C ). Marginally, without including any genetic information, from a pure environmental association analysis we observed a relationship between duration of OCP use and risk reduction for ovarian cancer. For 1–5 years of OCP use (vs <1 year) the estimated absolute risk difference was 0.47% (95%CI = 0.37% to 0.56%), while the estimated absolute risk difference for long-term use of OCPs (5+ year vs <1 year) was 0.84% (95%CI = 0.77% to 0.92%) ( Figure 2B-2C , eTable 4 ), in agreement with previous findings that longer duration of OCP use is associated with larger risk reduction in ovarian cancer 3 . However, when stratified by rs13255292 genotype, we observed an interesting pattern. Among individuals with TT genotype, the corresponding absolute risk difference estimate for 1–5 year of OCP use (vs <1 year) was 0.69% (95%CI = 0.49% to 0.88%), whereas among individuals with CC genotypes the corresponding risk reduction estimate was 0.36% (95%CI = 0.22% to 0.50%), implying potential effect modification by the C allele at locus rs13255292 (P-value = 1.12 × 10 −2 ) ( Figure 2B , eTable 4 ). In contrast, the absolute risk difference is estimated at 0.95% (95%CI = 0.78% to 1.12%) for women with TT genotype and at 0.79% (95%CI = 0.69% to 0.90%) in women with CC genotype. This indicates that longer OC use is associated with greater risk reduction overall and the risk reduction might be even greater for women with the TT genotype than those with the CC genotype. From Figure 2B-2C we observe the interplay between “nature vs nurture” with risk due to germline genetic mutations offset by long-term use of a modifiable protective factor. This analysis also highlights the benefit of measuring duration of exposure as opposed to a coarse indicator of ever/never use. Prior work in G × E for ovarian cancer has focused solely on multiplicative interactions. We previously reported no departures from a multiplicative model with the first six risk loci identified through GWAS with a reduced set of exposures 3 . Follow-up work identified an interaction with menopausal estrogen therapy use and rs10069690 in the TERT gene 48 , but that finding was not replicated in the present analysis which included a larger set of studies. Fridley and colleagues have reported on G × E taking a candidate gene approach with several promising findings 49 . There are several studies in other cancers examining G × E on the multiplicative scale with limited success in identifying interactions, but to our knowledge, only prostate cancer and bladder cancer have been studied on the additive scale. In prostate cancer, suggestive additive interactions between vitamin D, confirmed genetic variants and risk have been identified 50 . In bladder cancer, additive interaction has been explored between confirmed genetic loci and smoking with risk of disease 31 . In this work the authors were able to demonstrate that the absolute risk of bladder cancer for current smokers varied from 2.9% to 9.9% based on the polygenetic risk score quartile. These results are similar to our findings on the additive scale with absolute risk differing based on genetics and hormone therapy use; an interesting next step for our work is to consider the polygenetic risk score for all of these confirmed ovarian cancer susceptibility alleles. There are several limitations of the current analysis. Though we considered both multiplicative and additive interactions, the logistic model in (M1) is linear in covariates and exposures. We ignored potential non-linearity and exposure x exposure as well as exposure x covariate interactions. Similarly, we ignored any higher order interactions. A completely non-parametric machine learning approach, based on a recursive partition of the predictor space may avoid misspecification of the model, but would lack interpretability from an epidemiologic and public health perspective. We also acknowledge that this exploration of interaction is purely statistical, a more causal interpretation in a biological sense will require functional validation. One may also want to explore G × E interaction with loci that are not significant at genome-wide threshold but are significant at a less stringent threshold or even conduct genome-wide G × E scans. The associations between ovarian cancer risk and some of the variants included here were limited to specific histotypes of ovarian cancer, however we have only presented results for all epithelial ovarian cancers combined. Developing histotype-specific risk stratification approaches is not feasible because for any given histotype the absolute risk is unlikely to ever reach an actionable threshold on a population level. In addition, risk reducing strategies are the same across histotypes and thus there is little benefit to considering histotype specific results from a precision prevention perspective. Heterogeneous associations between environmental risk factors and ovarian cancer risk by histology has previously been well characterized 3 , 23 , 27 . There is value in understanding histotype associations for disease etiology and mechanisms and this will be the focus of future work. The analyses presented here offer insight into potential biological mechanisms, opportunities for ovarian cancer risk stratification, and approaches to studying gene-environment interactions. Ideally, replication for the six promising findings would be undertaken, but this is challenging with ovarian cancer given that most studies with the relevant data are included here. Functional studies for the regions harboring our most promising findings are underway and it is possible that the association described here may help inform those investigations 51 . Also, gene-environment interaction analyses can also be used to identify novel genetic associations 51 and thus a deeper evaluation of variants that are still borderline significant, but do not exactly achieve a genome-wide threshold is warranted for subsequent G × E analysis. Of particular interest will be to conduct risk stratification and risk prediction analysis using a summative polygenic risk score and to conduct an agnostic genome-wide search for G × E interaction. Despite the limitations the comprehensive framework of data harmonization, imputation, screening test followed by characterization of effect and risk estimates that has been used in this analysis can serve as a robust model for future gene-environment interaction analyses.

Ovarian carcinoma (cancer) is a disease with high mortality; most women are diagnosed with advanced stage disease where five-year survival is less than 50% 1 . Effective screening modalities have been elusive 2 , and therefore primary prevention strategies remain the most promising avenue to minimize the incidence and mortality of ovarian cancer. Several factors consistently associated with reduced or increased risk have been identified for ovarian cancer, including some that represent opportunities for chemoprevention or surgical intervention. Factors associated with reduced risk include oral contraceptive pill (OCP) 3 use aspirin use 4 , tubal ligation 5 , parity 3 , salpingectomy 6 – 9 and bilateral salpingo-oophorectomy (BSO). Common germline genetic variation 10 – 20 , first-degree family history of ovarian cancer 21 , 22 , menopausal hormone therapy use 23 – 25 , greater body mass index (BMI) 26 and endometriosis 27 are risk factors for the disease. OCPs and aspirin use represent feasible chemoprevention strategies whereas salpingectomy is now recommended by many gynecologic societies as an ovarian cancer prevention approach for women seeking tubal sterilization, having a hysterectomy, or having other pelvic surgery. Average lifetime risk of ovarian cancer diagnosis for women in the U.S. is 1.3% 28 , but this number varies greatly depending on the composite exposure history of risk factors 29 . Pearce et al. estimated the lifetime risk for women in the general population ranges from 0.35% (95%CI = 0.29% to 0.42%) to 8.8% (95%CI = 7.1% to 10.9%) depending on exposure history for six factors: OCP use, parity, tubal ligation, endometriosis, first degree family history of ovarian cancer and genetic risk score quintile 29 . However, these lifetime risk estimates were limited to six risk factors and did not consider their interaction with individual genetic variants identified through genome-wide association studies (GWAS) 28 . The multiplicative scale is commonly used for gene-environment interaction (G × E) analysis. Additive interaction analysis has been suggested for case-control studies in many recent papers for a more mechanistic interpretation 30 – 34 . Validity of a truly multiplicative model implies existence of additive interaction when the two factors under consideration have non-null main effects 35 . Thus, failure to detect G × E interaction on multiplicative scale may imply there exists interaction on additive scale, but the ability to detect it depends on the sample size and the main and interaction effect sizes 35 . We present here our efforts to evaluate both multiplicative and additive gene-environment interactions in ovarian cancer using data from the international Ovarian Cancer Association Consortium (OCAC) comprising 17 case-control studies. We have included 28 common genetic variants previously associated with risk of ovarian cancer in genome-wide association analyses for our G × E analyses 36 . Environmental factors included in our analysis are OCP use, parity, tubal ligation, breastfeeding, menopausal hormone therapy, usual adult BMI, and endometriosis. A small number of studies in OCAC had data available on aspirin use and thus we have not included this risk factor in our analysis here. Among our list of environmental factors, BMI, OCP use, tubal ligation, breastfeeding, and menopausal hormone therapy are of special interest because they are modifiable targets for prevention.