Genome-wide significant risk associations for mucinous ovarian carcinoma.

Genome-wide association studies have identified several risk associations for ovarian carcinomas but not for mucinous ovarian carcinomas (MOCs). Our analysis of 1,644 MOC cases and 21,693 controls with imputation identified 3 new risk associations: rs752590 at 2q13 (P = 3.3 × 10(-8)), rs711830 at 2q31.1 (P = 7.5 × 10(-12)) and rs688187 at 19q13.2 (P = 6.8 × 10(-13)). We identified significant expression quantitative trait locus (eQTL) associations for HOXD9 at 2q31.1 in ovarian (P = 4.95 × 10(-4), false discovery rate (FDR) = 0.003) and colorectal (P = 0.01, FDR = 0.09) tumors and for PAX8 at 2q13 in colorectal tumors (P = 0.03, FDR = 0.09). Chromosome conformation capture analysis identified interactions between the HOXD9 promoter and risk-associated SNPs at 2q31.1. Overexpressing HOXD9 in MOC cells augmented the neoplastic phenotype. These findings provide the first evidence for MOC susceptibility variants and insights into the underlying biology of the disease.

Ovarian carcinomas (OC) caused approximately 140,000 cancer deaths globally in 2008[1]. Germline mutations in genes conferring high (BRCA1 and BRCA2)[2] and more moderate risk (e.g., TP53, BRIP1, RAD51D, RAD51C, MLH1, MSH2, PMS1, PMS2 and MSH6)[3-56] of OCs are among the best-defined genetic risk factors but explain only 10-15 percent of all OCs[6-8]. More recently, genome-wide association studies (GWAS) have identified multiple regions of the genome harboring common variants (minor allele frequency [MAF] > 0.05) conferring low risk (odds ratios [OR] < 1.5) of invasive OC[9-17]. However, it is increasingly recognized that OCs encompass multiple distinct disease histotypes[18] that vary in epidemiologic[19-21] and genetic[22] risk factors, somatic alterations[23,24] and clinical response to platinum-taxane based therapy[18]. Most of the known common risk alleles for OC confer susceptibility to the most common histotype, high grade serous (HGSOC), with one region also associated with the clear cell histotype at genome-wide significance[12]. There are as yet no reports of confirmed genome-wide significant susceptibility loci for the other main histotypes, mucinous and endometrioid OCs.

Mucinous ovarian carcinomas (MOCs) are characterized by multicystic tumors with conspicious amounts of intracellular mucin (usually ≥ 50 percent of the cytoplasm) in more than 90 percent of the tumor cells[25]. Historically, MOCs have been estimated to account for about 12 percent of all OCs but recent refinements in the morphologic assessment indicate that primary invasive MOCs comprise approximately 3 percent of all OCs[26]. This lower prevalence is due to several reasons including consensus by pathologists to separate benign mucinous tumors from invasive MOCs[27] and from pathology guidelines[26,28-30] that aim to distinguish primary invasive MOCs from metastatic carcinomas involving the ovary, in which the majority derive from organs of the gastrointestinal system[26,31]. These criteria, along with the frequent inability to find a non-ovarian primary cancer, suggest that true MOCs develop de novo at the ovary and cannot be explained by metastatic lesions. This low incidence has made it challenging to study the etiology and pathogenesis of these tumors.

At the genetic level, MOCs are not associated with germline BRCA1/BRCA2 mutations. Unlike other OC histotypes, invasive MOCs usually harbor foci of benign or atypical (low malignant potential [LMP]) epithelium, with identical KRAS mutations frequently present[32-34], suggesting that this is an early somatic event in a multistep progression model. Normal mucin-secreting cells are not present in the ovary raising uncertainty regarding the cell at risk of transformation. It has been hypothesized that some MOCs originate from foci of benign endocervical-subtype Müllerian metaplasia of the surface epithelium or cortical inclusion cysts[35]. This subtype, however, may be less frequently associated with fully invasive MOCs, which comprise mostly the intestinal subtype[35]. To complicate further the etiology of MOCs, expression analysis of small numbers of MOCs (N = 3–9) associated these tumors more closely to colonic epithelium or colorectal carcinomas (CRC) than to ovarian surface epithelium[36,37], suggesting the pathogenesis of MOCs may be similar to colorectal carcinomas[38]. The current study reports the identification of genetic susceptibility alleles for MOCs, which may help to elucidate genes and biological pathways that are disregulated during MOC development.

Results

Genetic association analyses

We used genotypes from 16,038 ovarian cancer cases and 30,816 controls from various genotyping arrays providing genome-wide coverage (). Participating studies are listed in [10,12,39]. We imputed these genotypes into a reference panel from the 1000 Genomes Project to provide observed or imputed genotypes at 15,504,273 variants (Online Methods, ). Genotype re-imputation without pre-phasing was carried out for regions of interest to improve accuracy (see ). The primary association analyses reported in this paper were based on OCAC-COGS participants of European ancestry and those with invasive or LMP MOC, comprising 1,644 cases (1,003 invasive, 641 LMP) and 21,693 controls (). We identified SNPs in three different regions that were associated with MOC at genome-wide significance (, ). Two regions (2q13 and 19q13.2) have not been previously associated with risk for other OC histotypes; the third region (2q31.1) has been reported to be associated with HGSOC[10].

At 2q13, the most strongly associated SNP, rs752590, was imputed (imputation r2 = 0.66, effect allele frequency, EAF = 0.21). It is located 347 bases upstream of PAX8 (paired box 8) and the effect allele was associated with increased risk for all MOC (OR = 1.34, 95% CI = 1.21–1.49, P = 3.3 × 10−8) (). The risk was similar for invasive and LMP cases (data not shown). At 19q13.2, the most strongly associated SNP, rs688187, was also imputed (imputation r2 = 0.55, EAF = 0.32). It lies approximately 489kb downstream of IFNL3 (interferon, lambda 3) and the effect allele was associated with decreased risk for all MOC (OR = 0.67, 95% CI = 0. 0.60–0.75, P = 6.8 × 10−13). Again there was little difference in risk between invasive and LMP cases (data not shown).

At 2q31.1, the most significantly associated SNP, rs711830 (EAF = 0.32), is located downstream of the 3' region of HOXD3 (homeobox D3). The effect allele was associated with increased risk for all MOC (OR = 1.30, 95% CI = 1.20–1.40, P = 7.5 × 10−12) () with similar risk for invasive and LMP MOC (data not shown). This SNP was also associated with invasive HGSOC (OR = 1.14, P = 1.9 × 10−13) (). It is highly correlated (r2 = 0.99) with rs2072590, the variant previously reported for HGSOC[10].

MOCs of extra-ovarian origin are more likely to be stage 3 tumors[26]. In our dataset, only 146/1,644 (8.9 %) MOC cases were stage 3 suggesting that the majority of diagnoses in this study are likely to be true primary ovarian MOCs. Risk estimates were also similar or larger (farther from the null) in women diagnosed with early-stage mucinous tumors compared to stage 3 tumors: OR = 1.39 (95% CI = 1.22–1.58, P = 5.4 × 10−7) for rs752590 at 2q13; OR = 1.28 (95% CI = 1.17–1.41, P = 6.7 × 10−8) for rs711830 at 2q31.1 and OR = 0.65 (95% CI = 0.56–0.75, P = 5.9 × 10−9) for rs688187 at 19q13.2, supporting our initial findings.

Assessment of imputation quality

We assessed the imputation accuracy between the SNPs achieving genome-wide significance by comparing the correlation between observed genotypes with estimated genotype doses from the imputation without pre-phasing for 2,739 OCAC cases of European ancestry (Online Methods). SNPs were selected for genotyping based on pre-phased imputation results. Where primer design failed and the original risk SNP could not be genotyped, a highly correlated (r2 = 1) alternate SNP was selected (). The correlation between observed and imputed genotypes was 0.59 for rs6542125 and rs6758928 (proxy SNPs for rs72831838, the top SNP at 2q13 for pre-phased imputed data) and 0.51 for rs35963157 (top SNP at 19q13.2 for pre-phased imputed data). These were close to the expected correlations based on the imputation r2 (0.64 and 0.55, respectively), and the similarity of the empirical to the imputation r2 indicates the results are unlikely to be false positives.

Risk associations based on observed genotypes for 2,739 cases could not be evaluated, as we did not genotype any control samples. However, SNPs at MOC risk regions at 2q13 and 19q13.2 were not associated with other histotypes of OC (). We, therefore, compared the results of a case-only analysis for 151 MOC cases and 2,588 cases with other invasive OC histotypes using imputed genotype doses and observed genotypes (). We confirmed associations for the two alternate genotyped SNPs at 2q13, rs6542125 (case-only OR = 1.42, 95% CI = 1.07–1.90, P = 0.01) and rs6758928 (case-only OR = 1.41, 95% CI = 1.06–1.88, P = 0.02). While the ORs from genotyped samples were somewhat attenuated compared to the results for the imputed data, the P-values were smaller suggesting that the imputed association was robust. However, for rs35963157, there was no case-only association using imputed data in this subset of genotyped cases (case-only OR = 0.98, 95% CI = 0.69–1.39, P = 0.93). For the same samples, the association using observed genotypes was in the direction expected based on the full dataset, although not significant (case-only OR = 0.88, 95% CI = 0.67–1.14, P = 0.33). The association was somewhat stronger when an additional 1,274 cases (59 MOC cases and 1,215 other OC histotypes) without imputed data were included (case-only OR = 0.80, 95% CI = 0.64–1.00, P = 0.05) suggesting the observed association for the imputed data is not due to artifacts of imputation.

Functional annotation of variants in risk regions

At each of the three risk regions we identified all SNPs with a 1:100 or greater statistical odds of being the disease-causing variant (). We annotated these SNPs with respect to exons, introns and untranslated regions (UTR) and epigenetic marks from two ovarian epithelial cell lines and two OC cell lines. Given the biological similarities between some MOC and CRC[40], we also annotated these SNPs for epigenetic marks profiled in a CRC cell line (HCT-116) and normal colonic tissues. The vast majority of the variants lie within non-coding DNA, suggesting they influence the function of non-coding regulatory elements.

At 2q13, there were 55 candidate SNPs spanning a 78.6 kb region encompassing most of PAX8 and part of PSD4 (pleckstrin and Sec7 domain containing 4) (). Most risk-associated variants were in PAX8 introns, but the most statistically significant SNP (rs752590) and a moderately correlated variant (rs4849174, r2 = 0.73) in an RNA polymerase II binding site lie within the PAX8 proximal promoter. Three SNPs (rs874898, rs1478 and rs1479) lie within the 3' UTR of PAX8 and so could influence RNA stability. Two SNPs (rs6734610 and rs7585510) lie within the sequence of the PAX8-AS1 (PAX8 antisense RNA 1) long non-coding (lnc) transcript and so may impact its stability or function (). Eleven (20 percent) and 13 (24 percent) of SNPs lie in enhancer elements detected in ovarian and colonic cells, respectively. Of these, rs2305132, lies within a CTCF binding site detected in multiple cell types, suggesting this variant may be involved in the repression of PAX8 and/or PSD4 expression during MOC development.

At 19q13.2, there are 14 candidate SNPs located in and around IFNL3 and IFNL4 (). Rs11882871 lies within the 3' UTR in IFNL3 and rs4803222 lies within the 5' UTR in IFNL4 suggesting they may influence RNA stability[41]. Rs11322783 is a coding SNP and predicted to cause a frameshift change in IFNL4, while rs8103142 encodes a missense change (lysine to arginine) in IFNL3. Rs8103142 is predicted to be non-deleterious by PolyPhen v2 and SIFT. There were no overlaps between these 14 risk SNPs and regulatory DNA elements.

At 2q31.1, there were 19 candidate causal variants spanning ~27 kb encompassing HOXD3 and the lnc RNA HAGLR (HOXD antisense growth-associated long non-coding RNA) (). There was extensive overlap between SNPs and regulatory elements in this region. Eleven and eight SNPs, respectively, coincided with putative enhancers in ovarian and colonic cells. Rs1051929 encodes a synonymous change in HOXD3 and five SNPs lie within transcribed regions of the HAGLR and HAGLR-OS lnc RNA genes.

eQTL analysis

We evaluated associations between risk SNPs and transcript expression for all genes within a 100 kb window centered on the most risk-associated SNP in each region using publicly available data for 339 HGSOCs and 121 CRCs from The Cancer Genome Atlas (TCGA)[42,43]. No data were available for primary MOCs in TCGA. Where genotyping data were not available for a risk-associated SNP, correlated proxies (r2 > 0.7) were evaluated. At 2q13, we detected a significant eQTL association between rs6542127, a variant highly correlated (r2 > 0.9) with 6 risk-associated SNPs, and PAX8 expression in CRCs (P = 0.03, FDR = 0.09) (). Rs6542127 was associated with MOC risk (OR = 1.20, P = 8.81 × 10−6). The most significant eQTL association with PAX8 expression at this locus was for rs2863243 (P = 2.2 × 10−6) but this SNP was not associated with MOC risk. However, the third most significant eQTL SNP (rs3748916, P = 3.1 × 10−4) was associated with MOC risk (OR = 0.84, P = 9.37 × 10−6). There were no statistically significant eQTL associations with PAX8 expression in HGSOCs.

At 2q31.1, the most significant risk SNP (rs711830) in HOXD3 was significantly associated with HOXD3 expression in CRCs (P = 0.01, FDR = 0.09) but there was no eQTL association for HOXD3 in HGSOCs. HOXD9, approximately 49 kb upstream of rs711830, is a candidate susceptibility gene for HGSOCs (K. Lawrenson et al, unpublished data). Rs711830 genotype was significantly associated with HOXD9 expression in both HGSOCs (P = 4.95 × 10−4, FDR = 0.003) and CRCs (P = 0.01, FDR = 0.09) (). Another SNP in the region, rs10188827 (r2 = 0.59 with rs711830), showed a slightly stronger eQTL association in HGSOC (P = 2.05 × 10−4) but a slightly weaker association with MOC risk (OR = 1.29, P = 3.41 × 10−10 for all mucinous cases). There was also a stronger association between rs10188827 and HOXD9 expression in CRCs (P = 0.003), although the strongest eQTL association in CRCs was for another SNP, rs973456 (P = 5.30 × 10−5), which had not been imputed. There were no eQTL associations for genes in the 19q13.2 region in either the HGSOC or CRC tumor datasets.

Functional characterization using in vitro models of MOC

The eQTL analyses in HGSOCs and CRCs suggested HOXD9 is a candidate susceptibility gene and target of MOC risk SNPs. We therefore evaluated the role of HOXD9 in MOC development. We used chromosome conformation capture (3C) to establish if any of the risk SNPs at 2q31 interacted physically with the HOXD9 promoter in MOC EFO-27 cell lines . We found interactions for DNA fragments containing rs2072590, rs2857532 and rs4972504. These interactions spanned 31–55kb of genomic DNA and were confirmed by sequencing (). These genotyped variants were highly correlated with rs711830 (r2 = 1, 0.98 and 0.89, respectively), the most significant risk SNP in the region. Of the three SNPs, rs2072590 showed the greatest overlap with epigenetic marks, coinciding with enhancer marks in OC cell lines, colon cancer cells and colonic crypts (). Taken together, the results of eQTL and 3C analyses indicated that DNA regions at 2q31.1 harboring MOC risk SNPs are involved in the regulation of HOXD9 expression. Future studies using genome-editing approaches to manipulate the different alleles of rs2072590 will be needed to evaluate the effects on both regulatory activity and HOXD9 expression and confirm the role of this SNP and HOXD9 in MOC development.

We also evaluated the effects of overexpressing HOXD9 in two MOC cell lines (EFO-27 and GTFR230) using lentiviral transduction of a full length HOXD9 GFP fusion construct. Overexpression of HOXD9 was confirmed by qPCR and immunofluorescence microscopy (). HOXD9 expression was only detected in the nucleus, whereas in control cells expressing GFP only, GFP signal was detected throughout the cell. HOXD9 overexpression induced a significant increase in anchorage-independent growth in both MOC cell lines (P = 0.02 in EFO-27 and P = 0.04 in GTFR230, ) indicating a role in neoplastic transformation. We observed no effect on cellular invasion and migration (data not shown).

Discussion

GWAS have identified common low-penetrance genetic susceptibility alleles for a multitude of common traits and diseases. As the size and scope of GWAS increase, so the power to identify risk alleles for rare disease subtypes has also increased. For OCs, the vast majority of confirmed risk associations from GWAS were for HGSOC, which accounts for almost two-thirds of all invasive OCs[24]. We report, for the first time, genome-wide significant risk associations for the rarer MOC histotype, identified as part of the largest genetic association study yet performed for OC.

Two of the three susceptibility regions identified for MOCs (2q13 and 19q13.2) are specific to this histotype, which may not be surprising given that MOCs are clinically and biologically distinct from the other OC histotypes. The third region associated with MOC risk (2q31.1) was previously reported as a susceptibility locus for HGSOC[10]. Similarly, the 17q12 risk region encompassing HNF1B was reported to be associated with HGSOC and the clear cell OC histotype[12] suggesting that the different OC histotypes have some degree of shared germline genetic etiology despite differences in somatic genetic alterations[22-24], epidemiologic risk factors[19-21] and response to standard chemotherapy[18]. This may reflect the influence of the site of tumor development (i.e., the ovary) and the possible functional role of risk alleles interacting with common processes involved in malignancy, such as the ovarian microenvironment.

We identified PAX8 at 2q13 and HOXD9 at 2q31.1 as candidate MOC susceptibility genes using eQTL analysis of primary HGSOCs and CRCs. CRCs share some molecular and histologic characteristics with MOC[40]; however, gene expression and the functional mechanisms of risk-associated SNPs may be tissue specific[44]. We were unable to perform eQTL analyses in normal tissues or primary MOCs due to the lack of publicly available datasets for this tumor histotype and uncertainty of the likely cell(s) of origin of MOCs. While the eQTL associations we identified were statistically significant, it is possible that eQTLs exist between MOC risk SNPs and other genes, either within these regions or regulated more distally. For example at 2q13 we also observed regional associations of similar statistical significance to PAX8 for PSD4 and a PAX8-antisense transcript LOC654433, although neither gene has been previously implicated in MOC development. We also found evidence of stronger eQTL associations for SNPs with weaker risk associations at both loci. In addition to disease heterogeneity, eQTL analyses are complicated by intra-tumor heterogeneity due to variation in copy number, methylation and gene expression. Thus, caution needs to be applied when interpreting eQTL data. Additional analyses in larger sample sizes and in tissues more relevant to MOC etiology will be needed to confirm the significance of HOXD9 and PAX8 as likely susceptibility genes for MOC at these loci.

Functional studies suggest that HOXD9 is the target MOC susceptibility gene at 2q31.1, through its interaction with three different regions harboring MOC risk SNPs and its ability to enhance neoplastic phenotypes when overexpressed in MOC cells. HOXD9 is also a candidate susceptibility gene for HGSOC (K. Lawrenson et al, unpublished data). The results from 3C analysis showed that one of the three HOXD9 interacting regions in MOC cells (containing rs4972404) also interacts with HOXD9 in HGSOC cells. This suggests that similar functional mechanisms regulating HOXD9 expression are acting in both MOC and HGSOC, but that the other two interacting regions were tissue-specific to MOCs, indicating regions that control regulation of HOXD9 in MOCs but not in HGSOCs. HOXD9 is a member of the HOX family of transcription factors that are only expressed during embryonic development to control patterning and differentiation. HOXD9 has not been well characterized in the context of cancer development, although the gene was aberrantly expressed in cervical cancer[45] and has been implicated as a marker of cancer stem cells in glioma[46].

The 2q13 MOC risk region has not previously been associated with risk of other diseases or traits. PAX8 is a plausible candidate susceptibility gene target at this locus. It encodes a transcription factor important in the development of the Müllerian duct[47] and may be a cell-type lineage marker that distinguishes carcinomas of gynecologic origin (e.g., ovary, uterus, peritoneum and fallopian) from other sites such as the gastrointestinal tract[23,47,48]. PAX8 is overexpressed in the majority of HGSOCs compared to normal ovarian epithelial cells[23], due partly from gene amplification[49], but is expressed in 10 percent[23] to 25 percent (L.E. Kelemen et al, unpublished data) of LMP and invasive MOCs and is not expressed in CRC cell lines[49]. Although the precise role of PAX8 in cancer development is unclear, PAX8 expression may be important for acquiring characteristics that maintain a malignant state, including repression of differentiation programs for specific tissue lineages[49-51].

The 19q13.2 risk region has been associated with impaired clearance of hepatitis C virus[52] and variation in response to hepatitis C therapy in Asians involving IFNL3 (also known as IL-28B)[53]. The most likely functioning risk SNP in this region is upstream of IFNL3 in the coding region of IFNL4: rs11322783 (also annotated as rs368234815 in dbSNP 141) is a predicted truncating variant, suggesting it has a loss-of-function role. The insertion allele turns IFNL4 into a polymorphic pseudogene and abolishes its activity[54]. The variant rs8103142 causes a nonsynonymous coding change (Arg70Lys) in IFNL3 but this change was predicted to be non-deleterious in vitro[55]. There are no reports implicating IFNL4 or IFNL3 specifically in the development of OCs or CRCs, although multiple reports have indicated a role for interleukins more broadly in OC. The 19q13 region has also been associated with structural rearrangements in OCs[56].

The novel risk associations we found at 2q13 and 19q13.2 were identified using imputed genotypes and were based on an estimated imputation r2 that was moderate for both SNPs. The imputation r2 is an estimation of the expected correlation between imputed genotypes and actual genotypes. Confirmation genotyping of imputed SNPs in a subset of the samples showed that the estimated imputation correlation was similar to the correlation between imputed and observed genotypes. Furthermore, case-only associations for these SNPs based on observed genotypes were associated with smaller P-values, providing support for the imputed genotype associations. While it is possible that imputation may be sensitive to small genotyping errors and differential with respect to case-control status, we would expect this bias to apply equally to all cases and not specifically to MOC cases. Since we did not observe significant associations with the other histotypes, the collective findings suggest that the associations with MOC were not due to biases of the imputation process.

The relatively large number of invasive and LMP MOCs in this study represents a significant strength. We combined genotyping data from patients diagnosed with invasive and LMP MOCs because these tumors are thought to evolve along a morphologic continuum[32-34]. In the three risk regions, the statistical significance of the risk associations was stronger in the combined analysis compared to the LMP dataset alone. However, molecular epidemiologic studies are limited by access to details on clinical presentation and by the difficulty to perform centralized histologic review. Although it is reasonable to assume that most LMP tumors arose primarily in the ovary, review of the histology and clinical records of the cancers might have led to exclusion of some cases that were metastases to the ovaries from non-ovarian primary cancers. To address this, we evaluated risk associations for early-stage MOCs separately and observed similar or larger effect estimates. However, even when all the relevant clinical information is available and immunohistochemistry is performed, it is sometimes impossible to be certain whether an MOC arose in the ovary or elsewhere[40]. Notably, we found no overlap between the risk associations we reported for MOC and those discovered in GWAS of gastrointestinal cancers[57] suggesting that the invasive MOCs in this study were mostly primary OCs rather than metastases.

In summary, we reported the first genome-wide significant alleles to be identified for MOC. The power to detect risk associations for MOCs has so far been limited by the small numbers of MOC cases collected through OCAC. The experiences of GWAS for more common cancers (e.g., breast[58] and prostate[59]) indicate that with larger number of MOC cases, we would expect to identify additional susceptibility alleles for MOC. The functional evaluation that we performed for the MOC risk-associated regions also suggests that future studies are likely to provide new insights into the understanding of the biology of MOC. Finally, because MOC and HGSOC are distinct OC histotypes and can be considered separate diseases[23,24], the identification of HOXD9 as a potential gene target showing oncogenic characteristics in both MOC and HGSOC can be considered independent evidence for the general role of this gene in oncogenesis.

Online Methods

Genetic association studies

We used genotypes from samples of European ancestry available from several Ovarian Cancer Association Consortium (OCAC) genotyping projects. Data were available for five population-based GWAS of ovarian carcinoma. These included 1,785 cases (206 MOC) and 6,118 controls from a UK-based GWAS (“UK GWAS”)[9], 2,166 cases (116 MOC) and 2,564 controls from a GWAS from North America (“US GWAS”)[61] including two smaller GWAS from New England-Brigham and Women’s Hospital (“NEC/BWH”) and the National Cancer Institute-Polish study (“NCI-POL”), and 467 cases (36 MOC) and 441 controls from another GWAS from North America (“Mayo GWAS”)[17]. An additional 11,620 cases (1,286 MOC) and 21,693 controls from 41 OCAC studies were genotyped using the COGS array (“OCAC-COGS”)[12]. The UK and US GWAS comprised several independent case-control studies, and samples from some of these studies were also subsequently genotyped using the COGS array. After removing duplicate samples, remaining samples represented 43 studies from 11 countries including 16,038 women diagnosed with invasive ovarian carcinoma, 1,644 of whom were diagnosed with MOC, and 30,816 controls from the general population. Informed consent was obtained in each of the individual studies and local human research investigations committees approved each study. Details of the genotyping design for each dataset are shown in . Further details of the component studies are found in .

Genotyping and quality control

The final sample sizes for each of the five GWAS and OCAC-COGS are shown in . Details for genotyping and quality control for each dataset are found in the .

Imputation

To impute unobserved genotypes of common variation across the entire genome, we performed imputation separately for OCAC-COGS samples and each of the OCAC GWAS samples that passed QC ( and ). We imputed variants from the 1000 Genomes Project reference panel (Integrated Phase 1, version 3, March 2012 release). To improve computation efficiency we initially used a two-step procedure, which involved pre-phasing in the first step and imputation of the phased data in the second. We carried out pre-phasing using the SHAPEIT software[62]. Unobserved genotypes were inferred probabilistically with IMPUTE2 version 2[63]. To perform the imputation we divided the data into nonoverlapping segments of approximately 5Mb each. We excluded SNPs from the association analysis if their imputation accuracy was r2 < 0.25 or their MAF was < 0.005. In total, of 15,504,273 SNPs were imputed successfully ().

Re-imputation without pre-phasing

Following detection of a genome-wide significant association, we re-imputed regions of 1 Mb surrounding the SNP using IMPUTE2 software[63] without pre-phasing and the most recent 1000 Genomes Project reference panel (June 2014 data release) to improve accuracy. To increase the imputation accuracy even further, we changed some of the default parameters in the imputation procedure. These included an increase of the Markov Chain Monte Carlo iterations to 90 (out of which the first 15 were used as burn-in), an increase of the buffer region to 500 Mb and an increase of the number of haplotypes used as templates when phasing observed genotypes to 100. These changes were applied consistently for all datasets. This two-stage process resulted in different SNPs achieving the most significant P value between pre-phased and re-imputed data. The main findings focus on the associations obtained from re-imputation.

Imputation quality

Lymphocyte DNA from 4,013 OCAC cases of European ancestry, of which 2,739 were represented in COGS, was genotyped using the iPLEX Gold assay (Sequenom, Inc) in order to compare the accuracy of imputed and actual genotype doses. SNPs showing genome-wide significance were selected for genotyping from the pre-phased imputed results. Where primer design failed for the imputed pre-phased GWAS SNP, an alternate SNP was selected for genotyping within a region of high linkage disequilibrium (LD). Using the r2 coefficient, we compared the imputation accuracy between genotype doses derived from imputation using both the pre-phased and re-imputed results and those derived from actual genotyping.

Statistical analysis

Association testing was restricted to OCAC-COGS participants of European ancestry and those with invasive or LMP MOC, resulting in 1,644 cases (1,003 invasive-only and 641 LMP) and 21,693 controls. Genotypes (both typed and imputed) for each SNP were used to estimate allele frequencies and pair-wise LD between two variants was estimated with r2 values based on the 1000 Genomes Project reference panel using an online program (see URLs).

We estimated per-allele log odds ratios (OR) and 95% confidence intervals (CI) between each SNP and MOC risk using unconditional logistic regression, where number of variant alleles carried was treated as an ordinal variable (log-additive, co-dominant model). The likelihood ratio statistic was used to examine association; this statistic has been shown to have greater power for rare variants than alternatives such as the Wald test or the score test [64]. The logistic regression model was adjusted for study stratum and population substructure by including the first five eigenvalues from principal components analyses (see ref. [12]). Analyses were performed separately for combined invasive and LMP MOC and invasive-only MOC. Statistical P-values were two-sided and, unless stated otherwise, were implemented with STATA version 13.0 (StataCorp LP) and SAS version 9 (SAS Institute). P-values < 5 × 10−8 were considered to be genome-wide significant.

Molecular and functional analyses

Analysis of epigenetic biofeatures

We used in-house FAIREseq and H3K27ac plus H3K3me1 ChIPseq data for “normal” ovarian epithelial cell lines and serous ovarian cancer cell lines (G. Coetzee et al, unpublished data). HCT116 and primary tissue data (colon and ovarian) were downloaded as BED files from Hnisz et al[60]. CTCF data were obtained from ENCODE using the UCSC Genome Browser (see URLs).

Expression quantitative trait locus (eQTL) analyses

For each MOC SNP, we first identified correlated variants (r2 > 0.7) in the 1000 Genomes European ancestry population. Germline genotypes of 339 high grade serous ovarian cystadenocarcinoma and 121 colorectal cancer samples were downloaded from The Cancer Genome Atlas (TCGA) data portal and samples were selected for inclusion in the eQTL analyses using EIGENSTRAT[65]. For each sample, tumor gene expression profiles, somatic copy number and CpG methylation data were downloaded from TCGA portal (see URLs). Expression profiles were adjusted for somatic copy number and CpG methylation variation as previously described[66,67]. For each SNP we evaluated the correlation between the genotype and adjusted expression levels of candidate genes. A false-discovery-rate (FDR) below 0.1 was considered to be a significant cis-association.

Cell culture

The EFO-27 cell line is commercially available and was originally derived from a solid metastasis of a mucinous papillary ovarian carcinoma[68]. The GTFR230 mucinous ovarian cancer cell line was created in house and derived from a stage IC low grade primary mucinous ovarian cancer collected as part of the Gynecological Tissue and Fluid Repository (GTFR) at the University of Southern California (USC). We routinely profile cell lines using STR profiling (using the Powerplex 16 panel) at the University of Arizona Genetics Core facility and screen for mycoplasma using a mycoplasma-specific PCR. A piece of the tissue removed during surgery was transferred to the cell culture laboratory, minced into small pieces (1–2 mm diameter) and placed into NOSE-CM[69] which consists of Medium 199 mixed in a 1:1 ratio with MCDB105 (Sigma), 15% fetal bovine serum (FBS, Hyclone), 10 ng/ml epidermal growth factor (EGF), 34 μg/ml bovine pituitary extract (Life Technologies), 500ng/ml hydrocortisone and 5 μg/ml insulin (Sigma). Cells were confirmed to be epithelial in origin by staining for cytokeratin.

Chromosome conformation capture (3C)

3C was performed as previously described[70]. A sample of 10 million EFO-27 cells was fixed with 1% formaldehyde, then lysed in 10 mM Tris-HCl pH 8, 10 mM NaCl, 0.2% Nonidet P-40. Nuclei were pelleted and then resuspended in 1X restriction enzyme buffer with 0.1% SDS and 1.6% Triton-X. Fifteen hundred units of Csp6i (Fermentas) enzyme were added. Samples were incubated at 37 °C overnight to digest after which 1.5% SDS was added before de-crosslinking samples by incubating at 65 °C for 30 minutes. Digested genomic DNA was then ligated by incubating samples in ligation buffer plus 4000U T4 DNA ligase (NEB) for 24 hr at 16°C. Control samples received no ligase. Following ligation, samples were de-crosslinked by overnight incubation at 65 °C with proteinase K. Chromosome conformation capture libraries were extracted using phenol/chloroform, the DNA precipitated using ethanol, and desalted using Microcon Ultra Cell YM –100 columns. Two primers were designed at the HOXD9 promoter and one primer for each restriction fragment of interest (Supplementary Table 9). PCR was performed using HotStar Taq polymerase (Qiagen), using the following conditions: 5 min at 94 °C, 42 cycles of (20s at 94 °C, 20s at 57–59 °C, and 30s at 72 °C), then 10 min at 72 °C. PCR products were run on a 1.2% agarose gel. PCR products were either cleaned up using the QIAgen QIAquick PCR Purification kit or were gel purified using the QIAgen QIAquick Gel Extraction kit. Cleaned products were sequenced from both ends by Genewiz.

In vitro modeling of HOXD9 overexpression

HOXD9-GFP and GFP lentiviral vectors were purchased from Genecopeia. Lentiviral supernatants were produced by co-transfecting HEK293T cells with vectors of interest, plus lentiviral packaging vectors. Supernatants were snap frozen and stored at –80°C until use. EFO-27 and GTFR230 cells were transduced with HOXD9-GFP and GFP viral supernatants in the presence of 4 g/ml polybrene (Sigma). Cells expressing GFP alone were controls. RNA was extracted from selected cells using the Zymo Quick-RNA kit, reverse transcribed using the MMLV enzyme (Promega) and gene expression analysis performed using TaqMan probes (Life Technologies). Gene expression for HOXD9 was normalized to expression of GAPDH and beta-actin using the delta-delta Ct method. For immunofluorescence, cells were fixed on glass coverslips using 4% paraformaldehdye and nuclei stained using Hoechst. Cells were imaged using a Zeiss Axio Imager Z1 fluorescent microscope. Positive stained cells were selected with puromycin at concentrations of 1 mg/ml (EFO-27) and 400 ng/ml (GTFR230). Anchorage-independent growth assays were performed in 6-well plates by re-suspending 20,000 cells per well in media containing 0.33% Noble agar and 1 mg/ml bacto-peptone (both Sigma) and plating atop 3 ml of media containing 0.6% Noble agar and 1 mg/ml bacto-peptone. Assays were incubated at 37 °C for 28 days before fixing and staining with 1 mg/ml p-iodo tetrazolium violet (Sigma) dissolved in 100% methanol (VWR). For invasion and migration assays, Cultrex 96-well invasion and migration kits were used according to the manufacturer’s instructions.

Table 1

Summary of genotyping datasets used for imputation*, European samples

Study Set[†]	Controls, N	All Ovarian Cases, N	Mucinous Cases All (Invasive-only)	Genotyping Platform	Genotyping Center	Number SNPs passing QC
U.K. GWAS	6,118	1,785	206 (184)	Illumina 610K (cases) Illumina 1.2M and 550K (controls)	Illumina Corporation Sanger Centre	507,094 507,094
U.S. GWAS	1,867	1,813	100 (99)	Illumina 610K	Mayo Clinic Medical Genome Facility	556,480
NEC/BWH GWAS	142	132	0	Illumina 317K	National Cancer Institute	305,690
NCI-POL GWAS	555	221	16 (13)	Illumina 550K	National Cancer Institute	527,435
Mayo GWAS	441	467	36 (11)	HumanOmni 2.5 BeadChip	Mayo Clinic Medical Genome Facility	1,587,042
OCAC- COGS	21,693‡	11,620‡	1,286 (696)‡	Illumina custom iSelect ~211K	McGill University and Gέnome Quέbec Innovation Centre and Mayo Clinic Medical Genome Facility	199,570
Total	30,816	16,038	1,644 (1,003)

All datasets were used for imputation; however, association analysis was based on all mucinous ovarian cancer cases (N=1,644) and controls genotyped in OCAC-COGS (N=21,693)

See Supplementary Table 1 for details of individual studies

Number of unique samples after exclusion of duplicates also genotyped in the other five GWAS

Table 2

Association testing in OCAC samples participating in COGS, 1,644 MOC cases and 21,693 controls

						Ref/Alt Allele	Invasive+LMP MOC N=1,644			Invasive-only MOC N=1,003
SNP	Locus	Position	Nearest Gene	EAF*	Imputation r2 †		OR	95% CI	P	OR	95% CI	P
rs752590	2q13	113972945	PAX8	0.21	0.66	A/G	1.34	1.21–1.49	3.3×10−8	1.31	1.17–1.46	2.9×10−6
rs72831838 ‡	2q13	114016401	PAX8	0.14	0.65	C/T	1.38	1.23–1.56	1.2×10−7	1.39	1.22–1.58	7.5×10−7
rs6542125 (alternate) §	2q13	114018826	PAX8	0.15	0.64	G/A	1.38	1.23–1.56	9.0×10−8	1.37	1.21–1.56	1.2×10−6
rs6758928 (alternate) §	2q13	114019350	PAX8	0.15	0.64	G/A	1.38	1.23–1.55	9.4×10−8	1.37	1.21–1.56	1.3×10−6
rs711830 ¶	2q31.1	177037311	HOXD3	0.32	1	C/T	1.30	1.20–1.40	7.5×10−12	1.26	1.15–1.39	7.7×10−7
rs688187	19q13.2	39732752	IFNL3	0.32	0.55	G/A	0.67	0.60–0.75	6.8×10−13	0.71	0.63–0.80	1.3×10−8
rs35963157 ǀǀ	19q13.2	39745695	IFNL4	0.34	0.55	G/GC	0.68	0.61–0.76	3.0×10−12	0.72	0.64–0.81	2.9×10−8

Effect allele frequency

Imputation r2 was estimated without pre-phasing using the 1000 Genomes Project reference panel

Most significant SNP on 2q13 in pre-phased imputation results; ranked #19 without pre-phasing

Alternate SNPs on 2q13 were correlated with rs72831838 (r2=1) and selected for genotyping in lymphocyte DNA due to failed primer design for rs72831838; rs6542125 is ranked #13 and rs6758928 is ranked #14 without pre-phasing

Genotyped SNP

Most significant SNP on 19q13.2 in pre-phased imputation results; ranked #8 without pre-phasing