Functional Analysis and Fine Mapping of the 9p22.2 Ovarian Cancer Susceptibility Locus.

Here, we utilized two independent approaches to identify a list of candidate causal SNPs ( Fig. 1 ). First, we conducted a comprehensive analysis to identify functional SNPs in linkage disequilibrium (LD, r 2 > 0.3) with rs3814113 with no prior assumption about their individual association to risk. Since all resided in non-coding regions, we hypothesized that SNP alleles determine the activity of regulatory elements in enhancers and promoter regions active in OSE and fallopian tube surface epithelial cells (FTSEC) ( 29 ). Second, we performed fine-mapping association analyses by densely genotyping over 15,000 ovarian cancer cases and 30,000 controls to identify a credible set of causal SNPs guide by association data. These parallel approaches identified the SNP most likely to be causal to ovarian cancer risk at the 9p22 locus. A total of 134 SNPs were chosen for functional analysis, based on their LD (r 2 ≥ 0.3) with rs3814113 in European 1000 Genomes Project data (v3 April 2012 release). They are distributed over an 82 kb region ranging from the first intron of BNC2 to ~44 kb centromeric to its transcription start site (TSS) ( Fig. 2a and Supplementary Table 3 ). Since all SNPs in the candidate functional set are in non-coding regions, several independent assays were used to identify transcriptional regulatory elements. First we analyzed data from FAIRE-Seq, and ChIP-Seq for H3K27Ac and H3K4Me1. FAIRE-Seq reveals regions of open chromatin while H3K27Ac or H3K4Me1 are markers for active chromatin and enhancers, respectively. The chromatin landscape profiles ( Fig. 2a ) were derived from iOSE and iFTSEC cells ( 16 ). Analysis of FAIRE- and ChIP-Seq data identified twelve regions with evidence of enhancer activity in at least one cell line ( Fig. 2a ). Twenty-two candidate causal SNPs ( Table 1 ) are located within five regions containing FAIRE or ChIP-Seq features suggesting that these SNPs might have a functional impact ( Fig. 2a ). The relatively lenient threshold for LD and criteria to consider a region as a putative enhancer was designed to favor sensitivity at the initial stage of analysis (with high specificity being achieved by the integration of the two approaches). To refine the analysis, we tested twelve genomic tiles (~2 kb each) ( Fig. 2a ), in both orientations, spanning the five candidate regions using a reporter assay to identify enhancer activity in iOSE4 cMYC ovarian cells ( 17 ). Although not present in a region with evidence of regulatory activity, we also tested one tile containing rs3814113 (Tile 12), the most significantly associated with high grade serous EOC in a previous study ( 2 ), and a control tile devoid of evidence for enhancer activity as judged by FAIRE and ChIP-Seq data ( Fig. 2a , Tile C). Tiles in regions 6 (T6), 7 (T7.2, T7.3, T7.6), and 8 (T8) contained nine candidate causal SNPs and showed significant activity (two tailed t- test p<0.05 compared to the control tile C; two replicates) in at least one orientation ( Fig. 2b ). Causal SNPs are hypothesized to display allele-specific effects. Therefore, we used site-directed mutagenesis in tiles T6, T7.2, T7.3, T7.6, and T8 to change each of the nine candidate causal SNPs from the reference to the effect allele and compared their activity. For tiles with multiple SNPs, the reference tile was the most common haplotype ( Supplementary Fig. 1 ) (All populations; 1000 Genomes Project). Individual SNPs were mutated to determine the contribution of each SNP, with other SNPs in the haplotype retaining the reference SNP allele. Seven SNPs in T6 (rs62541878), T7.2 (rs62541920, rs12379183), T7.3 (rs1092647), and T8 (rs77507622, rs10810657, rs12350739) demonstrated significantly different transcription activity (p<0.05) between the reference and effect allele in at least one replicate ( Fig. 2c-d ). These seven SNPs were retained for analysis. We conducted EMSAs using probes with both alleles for each of the seven SNPs in regions 6, 7 and 8 ( Fig. 2e ). Tiles that did not show activity in Fig. 2b were not tested. Tile 11 had significant transcription activity in only one reporter experiment but two SNPs within the region (rs113780397 and rs181552334) are correlated with the original SNP (r 2 of 0.818 and 0.5, respectively), and so four additional probes were tested. We also examined rs3814113, the most significant original GWAS SNP. EMSAs revealed allele specific nuclear extract binding for rs12379183, rs62541920 (Region 7), rs12350739, rs77507622 (Region 8) and rs181552334 (Region 11) ( Fig. 2e ) indicating these SNPs were strong causal candidates. Region 11 overlapped with an open chromatin region, according to FAIRE-Seq data obtained in ovarian cells, and one SNP showed allele-specific binding in EMSA experiments (rs181552334). However, this region lacked H3K4Me1 and H3K27Ac marks and luciferase assays showed weak evidence for enhancer activity in ovarian cells ( Fig. 2b ). Interestingly, the region is A/T rich (> 60%), a feature in regions that anchor the cell’s DNA to the nuclear scaffold/matrix ( 19 ). Moreover, Region 11 was predicted by MAR-Wiz to attach to the nuclear scaffold/matrix compared to the rest of the locus ( Supplementary Fig. 2 ). To determine whether Region 11 was attached to the nuclear scaffold in ovarian cells, we performed a nuclear scaffold extraction in iOSE11 cells ( 19 ), using HeLa cells as a control. Region 11 had significantly higher enrichment in the scaffold fraction of iOSE11 and HeLa cells than a previously defined negative control (ApoB Neg) ( 19 ) ( Fig. 3a-c ). A region previously defined as a S/MAR (ApoB S/MAR) ( 19 ) in HeLa cells did not have significantly higher enrichment in the scaffold fraction of iOSE11 cells than ApoB Neg ( Fig. 3a-c ) but had significantly higher enrichment in the scaffold fraction of HeLa cells than ApoB Neg. These results indicate that Region 11 acts as S/MAR in ovarian cells. Visual inspection of HiC (High dimensional chromosome conformation capture) data from seven cell lines suggested the presence of a 1.8 Mb (chr9:15,750,000–17,550,000) topologically associating domain (TAD) in which the S/MAR (Region 11, rs181552334) is situated close to one of its borders ( Fig. 3d-e ). This TAD includes TSS for BNC2, C9orf92 , and CNTLN. Two functional SNPs in Region 7 were located in an approximately 7kb region that includes the TSSs for two BNC2 transcripts ( Fig. 2a ) denoted by FAIRE-Seq and H3K4me1 ChIP-Seq data in ovarian cells, and ENCODE layered H3K4me3 (promoters) ChIP-Seq data ( Fig. 3f ). This region is the major BNC2 promoter, raising the hypothesis that BNC2 may act as the mediator of risk at the 9p22.2 locus. Region 8, containing two SNPs with allele specific activity in luciferase assays and EMSA, overlapped with FAIRE-Seq and ChIP-Seq data in ovarian cells with features indicative of an enhancer ( Fig. 2a ). To determine potential interacting promoters with the enhancer at region 8, we examined all genes ( c9orf92, BNC2, CNTLN and SH3GL2 ) within a stretch of 1 MB at either side of the region containing the candidate SNPs ( Fig. 3f ). First, guided by H3K4me3 marks in seven non-ovarian cell lines from ENCODE, we identified their promoters close to TSSs ( Fig. 3f ). Next, we inferred whether the gene was expressed in ovarian cell lines using H3K27ac as a marker of active promoters combined with analysis of transcript levels from RNA sequencing (RNA-Seq) data for ovarian and fallopian tube epithelial cells ( Fig. 3g ). This analysis indicated that BNC2 and CNTLN were expressed in ovarian cells, but c9orf92 or SH3GL2 were not ( Fig. 3g ). Next, we used Chromatin Conformation Capture (3C) to determine which promoters physically interacted with Region 8. In iOSE11 cells, Region 8, when compared to an adjacent site displayed two regions of frequent ( Fig. 3h ; blue arches) interactions. The interaction peak closer to the anchor is located upstream of the TSS but does not overlap with any known chromatin marks. The second interaction peak corresponds to region 7 considered to be the core promoter of BNC2 . No significant interaction was detected between Region 8 and the CNTLN TSS ( Fig. 3h ). As expected, no interaction was detected between the S/MAR in Region 11 and promoters in the region ( Fig. 3i ). The modules in Regions 7 and 8 appear to affect the major promoter of BNC2 and are a distal regulatory enhancer that physically interacts with the BNC2 promoter, respectively. Next, as part of our two-pronged approach, we conducted fine mapping of the 9p22 locus in 15,437 women diagnosed with invasive EOC and 30,845 controls ( Fig. 4a ). We evaluated the association between genotype and disease using logistic regression by estimating the associations with each additional copy of the minor allele (log-additive models) for 4,234 SNPs of which 2,418 were directly genotyped ( Supplementary Table 4 ). SNP rs3814113 remained the most statistically significant association (P = 2.10 × 10 −34 ) ( Fig. 4a ) with the minor allele [C] being protective. Next, we calculated the likelihood ratio of each SNP relative to the most significant SNP (rs3814113) of being the functional variant underlying the signal. For any given set of correlated associated SNPs, the strength of evidence was estimated by the log likelihood statistic from the logistic regression; thus difference in the log likelihood between the SNP with the strongest association and any other SNP provides a measure of the log odds in favor of the most significant SNP being the SNP that is truly driving the observed association. There were 40 SNPs with odds of 1:1000 or better and were considered to be credible candidates for mediating the observed association. They were all in strong LD (r 2 > 0.89). While 35 out of the 40 SNPs were part of the set of 134 SNPs assessed during functional analysis, five SNPs (rs34131140, rs112442786, rs113198237, rs199782476, and c9_pos16900214/rs62543587) were not ( Supplementary Fig. 3 ). SNPs rs199782476 and rs62543587 did not overlap with any biofeatures (FAIRE-seq, H3K4me1, H3K27Ac) suggesting they were not functionally relevant. The remaining three SNPs were part of Tile 11 (rs112442786), Tile 12 (rs34131140) or the control (TC) Tile (rs113198237). Tiles 12 and TC did not display significant activity in enhancer scanning ( Fig. 2b ) suggesting that rs112442786, which resided in Tile 11 and mapped to the S/MAR region, may be functionally relevant ( Supplementary Fig. 3 ). Among the 40 SNPs, six SNPs with odds ranging from better than 1:4 (rs112442786) to better than 1:200 (rs181552334) mapped to the S/MAR (Region 11) ( Fig. 4a ). We repeated the association analyses adjusting for rs3814113 to identify additional independent signals in the region ( Fig. 4b ). Nine SNPs were significant at P < 10 −5 ( Supplementary Table 4 ; Conditional tab) of which two (rs7848057, rs80039758) mapped to the S/MAR ( Fig. 4b ). In this group of 9 SNPs, rs10756825 and c9_pos16889285 were the next most significant associations (p = 2.3 × 10 −5 ; p = 6.23 × 10 -5 , respectively) and mapped close to the enhancer in Region 8 ( Fig. 4b ). Finally, several attempts to remove the S/MAR using CRISPR-based genome editing techniques were not successful, suggesting that deletion of this region may impact the viability of ovarian and fallopian tube cells. To identify eQTL associations for c9orf92, BNC2, CNTLN and SH3GL2 we searched the GTEx dataset for single gene eQTLs in all tissues (GTEx Analysis Release V7; dbGaP Accession phs000424.v7.p2; fallopian tube not included due to small sample size). Although all four genes displayed eQTL associations ( CNTLN = 11,039; C9ORF92 = 1; SH3GL2 = 361; BNC2 = 94) ( Supplementary Table 5 ) only BNC2 displayed eQTL associations with SNPs (rs10962662, rs10756823, and rs10124837; whole blood) present in our set of 40 credible candidate SNPs. Next, we searched for single SNP eQTL associations in all tissues for 40 credible candidate SNPs. The only three eQTL associations found were for SNPs rs10962662, rs10756823, and rs10124837 with BNC2. The data from the functional analysis and fine mapping data provide evidence that the candidate causal SNPs at the locus exert their effects in a 1.8Mb TAD with BNC2 as the most likely target gene at the locus. BNC2 has three pairs of C2H2 zinc fingers (ZF) raising the possibility that it recognizes specific DNA sequences and is involved in transcription regulation ( Fig. 5a ) ( 30 ). To identify DNA sequences recognized by BNC2, GST-tagged constructs of each ZF pair ( Supplementary Fig. 4 ) were expressed in bacteria and applied to a protein binding microarray (PBM) with overlapping, rationally randomized nucleotides, representing every possible motif up to 10 bp ( 21 , 22 ). When aligned the top ten scoring sequences for each ZF pair generated a sequence logo using position weight matrix scoring ( Fig. 5a ). The motifs for ZF1,2 and 5,6 were consistent with the predicted C2H2 “recognition code” ( 31 ). Binding for ZF3,4, which yielded lower-confidence data, did not match the recognition code predictions ( Fig. 5a ) ( 32 ). The 3’ end of the ZF1,2 and ZF5,6 binding motifs had the same nucleotides at the exact same position and weight, consistent with the similarity in amino acid residue positions between ZF2 and ZF6 ( Fig. 5a ). Notably, the BNC2 promoter region contains two BNC2 ZF5,6 PBM binding sequences ( Supplementary Fig. 4 ). We validated BNC2 binding sequences identified with the PBM by conducting ChIP in iOSE11 and iFTSEC283 cells for endogenous BNC2 ( Supplementary Fig. 4 ) at the PBM sites (−582 and −914 bp upstream of the TSS) at the BNC2 locus. A significantly larger amount of DNA was immunoprecipitated with the BNC2 antibody than with the IgG control at the −582 (iOSE11 p = 2.6 × 10 −3 , iFTSEC283 p = 8.3 × 10 −3 ) and −914 (iOSE11 p = 1.8 × 10 −4 , iFTSEC283 p = 2.0 × 10 −6 ) bp sites, but not at the −2184 bp site (negative control; Supplementary Fig. 4 ). These data provides evidence that the sites identified in the PBM experiment are recognized by endogenous BNC2. To identify genomic sites bound by BNC2 in ovarian cells, we used ChIP-Seq in iOSE11 and iFTSEC283 cells (see Extended Data ). MEME, a motif analysis tool, defined a motif centrally enriched in the ChIP-Seq peaks in both cell types ( Fig. 5b-c ). The motif identified by MEME appears to be a concatenation of the reverse complement motif for ZF1,2 and the motif for ZF5,6 with a 75% homology ( Fig. 5b ). The concatamer motif was significantly enriched in ChIP-Seq peak summits in iFTSEC283 and iOSE11 cells ( Fig. 5c ). ChIP-Seq data replicated BNC2 binding in the iOSE11 cells (chr9:16871799–16872039) at the −914 position tested in ChIP-qPCR ( Supplementary Fig. 4 ). To identify putative target genes regulated by regulatory elements containing BNC2 recognition sites, we generated a list of 995 genes/transcripts with TSS within 30 kb of the BNC2 ChIP-Seq peak centers found in both iOSE11 and iFTSEC283 cells ( Supplementary Table 2 ). Next, we used PANTHER ( 33 ) and found that several functional classes were statistically overrepresented in our set including system development (GO:0048731), anatomical structure development (GO:0048856), single-multicellular organism process (GO:0044707), multicellular organism development (GO:0007275), and tissue development (GO:0009888) ( Supplementary Table 2 ). From the above set, we selected a set of 87 genes that were: a) implicated in ovarian cancer; b) ovarian development; c) were part of KEGG pathways related to cancer; d) in which BNC2 ChIP-Seq peaks were found in their core promoter (within 1kb from the TSS) ( Supplementary Table 2 ) and tested the extent to which their expression (measured by Nanostring) was modulated by overexpression of BNC2 in HEK 293T.. Multiple unsuccessful attempts were performed to manipulate expression - silencing or ectopic overexpression - levels of BNC2 in ovarian cells, suggesting that BNC2 levels are tightly controlled. Several genes mapping to KEGG Focal Adhesion, ECM-receptor interaction or TGF-β Signaling Pathways and implicated in ovarian cancer or ovarian development showed significant changes in expression upon BNC2 overexpression ( Table 2 ; Supplementary Table 2 ). Although most genes showed a positive correlation with BNC2 overexpression, FEM1A and IGTB5 showed an inverse correlation suggesting that BNC2 modulation of expression is likely to be context dependent ( Table 2 ). Taken together, these experiments validate the BNC2 binding site in vivo and reveal putative downstream targets of BNC2 activity in ovarian cells.

We used two immortalized normal OSE cell lines, iOSE4 and iOSE11 ( 14 ), and three immortalized normal fallopian tube surface epithelial cells (iFTSEC33, iFTSEC246, and iFTSEC283), a normal epithelial ovarian cell line, iOSE4 CMYC , immortalized with hTERT and transformed with MYC ( 15 ), and HEK293FT cells. Cell line aliquots were tested for mycoplasma (PCR-based method) and authenticated using STR analysis before being used for experiments, which were conducted before 20 passages after thawing. FAIRE-Seq (Formaldehyde Assisted Isolation of Regulatory Elements followed by sequencing) and ChIP-Seq (Chromatin immunoprecipitation followed by sequencing) for Histone H3 Lysine 27 Acetylation (H3K27Ac) and Histone H3 Lysine 4 Monomethylation (H3K4me1) were performed in iOSE4, iOSE11, iFTSEC33, iFTSEC246, iFTSEC283 ( GSE68104 ) ( 16 ). We used an optimized method to identify genomic regions with enhancer activity ( 17 ). Genomic tiles of ~2 kb were generated by PCR using bacterial artificial chromosome (BAC) Clone RPCI-11–185E1 (Empire Genomics) as the template and cloned in forward and reverse orientations upstream in the firefly luciferase reporter vector designed to test for enhancer activity ( 17 ). Primers can be found in Supplementary Table 1 . Transfections included a plasmid expressing Renilla sp. luciferase as internal control and every tile was tested in two independent experiments. Tiles with significantly (two-tailed t-test; p <0.05) higher luciferase counts than the control tile (TC) were tested for allele specific effects. For allele-specific luciferase assays, tiles with the effect allele were considered significant if the luciferase counts were significantly higher or lower (p <0.05) in at least one independent experiment than the tile with the reference allele. Nuclear extracts were obtained from iOSE4 CMYC cells at 70–90% confluence and EMSAs were run as previously described ( 18 ). A lithium-based nuclear scaffold extraction was performed as previously described ( 19 ). Scaffold and genomic DNAs were quantified by qPCR using primers for Region 11, the ApoB S/MAR and the ApoB Neg regions ( Supplementary Table 1 ). Samples were run using Sybr Green Spectrum on Applied Biosystems 7900 HT Real-Time PCR System. Enrichment was calculated by dividing the quantity of the scaffold DNA by the quantity of the digested genomic DNA. A Z-score for the region 11 and ApoB S/MAR was calculated as described previously ( 19 ) [Z score = (average of S/MAR – average of ApoB Neg)/std dev of ApoB Neg)]. A Z-Score > 8 indicates a site positive for scaffold binding ( 19 ). Each experiment includes three technical replicates. 3C libraries were prepared as previously described ( 20 ). qPCR was performed by using Taq Polymerase PCR Kit (Qiagen) and Syto9 (Life Technologies). Samples were run using FAM Spectrum on an Applied Biosystems 7900 HT Fast Real-Time PCR System. EcoR1 digested BACs (RPCI-11–185E1 Empire Genomics, RPCI-11–179K24 Life Technologies, RPCI-11–106G11 Life Technologies) for the region were used for the standard curve. Interactions were calculated as a percentage of a restriction site directly adjacent to the bait restriction site. Sites with a significantly higher frequency of interaction than the site adjacent to the anchor were considered significant (p <0.05; two-tailed t -test). 3C was performed in two independent experiments and three technical replicates each. Fragments containing cDNAs of each of the zinc finger pairs were PCR amplified from a plasmid containing BNC2 cDNA (a gift from Dr. Philippe Djian) using primers containing Gateway recombination sites ( Supplementary Table 1 ). PCR products were cloned into pDONR221 using the BP recombination kit and transferred to pDEST15 as a fusion to Glutathione-S-transferase (GST) using LR recombination kit (Invitrogen). Purified GST-ZFs were eluted from beads with 50 mM reduced glutathione and 0.5 μg of each GST-ZF protein construct were applied individually to two differently designed arrays designated ME and HK as previously described ( 21 , 22 ). ZFs typically bind to degenerate motifs and have the potential to have more than one recognition sequence ( 21 ). Each DNA probe sequence is given an E-score which is similar to the Area under the ROC curve statistical metric and an E-score above 0.45 was considered significant. Chromatin immunoprecipitations were performed as previously described ( 23 ) using a validated BNC2 antibody (Sigma Atlas) (see Supplementary Data ). Real-time qPCR was performed using Sybr Green chemistry with primers at the −2184, −914, and −582 positions relative to the TSS ( Supplementary Table 1 ) in an Applied Biosystems 7900HT Fast Real-Time PCR System. ChIP for each cell line was performed in four biological replicates. Overrepresentation test (release 20170413) was conducted with PANTHER version 11.1 released 2016–10-24 using all genes in Homo sapiens database as a reference list and a Bonferroni correction for multiple testing. The uploaded list contained 965 genes of which 839 were mapped to GO-Slim. For BNC2 ChIP-Seq four individual ChIP samples were pooled for each cell line (iOSE11 and iFTSEC283) in two biological replicates. Immunoprecipitated DNA was used to generate a sequencing library using the NuGEN Ovation Ultralow Library System with indexed adapters (NuGEN, Inc., San Carlos, CA). The library was PCR amplified and size-selected using AxyPrep Fragment Select beads (Corning Life Sciences – Axygen Inc., Union City, CA). Each enriched DNA library was then sequenced on an Illumina HiScan SQ sequencer to generate 20–30 million 100-base paired-end reads. The raw sequence data was de-multiplexed using the Illumina CASAVA 1.8.2 software (Illumina, Inc., San Diego, CA) and binding sites were identified using the MACS2 software ( 24 ) using input DNA as a control. See Supplementary Data for further details. pNTAP-BNC2 (or the empty vector) was transfected with Fugene 6 into 293FT cells at 70% confluence. Cells were harvested after 24 h, RNA was isolated using Trizol RNA Isolation (Life Technologies), and cleaned using Qiagen RNeasy Mini Kit (Qiagen). The three biological replicates for HEK293FT cells with the empty vector or over-expressed BNC2 were applied to a Nanostring platform containing probes for 87 genes and 10 reference genes ( Supplementary Table 2 ) used to normalize the data in the NanoString nSolver Analysis Software v 1.1. These genes had a %CV < 50. Genes were considered to be differentially expressed if p <0.05 (two-tailed t-test). To refine the observed signal at rs3814113 ( 2 ), fine-mapping was conducted using a customized Illumina iSelect genotyping array (iCOGS). SNPs were selected based on data from 1000 Genomes Project (1000GP) ( 25 ) CEU (April 2010) and Hapmap III within a 1 Mb interval of rs3814113 (chr9: 16407967–17407967)( 26 ). We included tagging SNPs (r 2 > 0.1) with a minor allele called at least twice in the 1000GP and additional SNPs tagging remaining variation in the interval (r 2 > 0.9), requiring Illumina Design score > 0.8. The iCOGS array was used to genotype cases and controls from constituent studies of the Ovarian Cancer Association Consortium as previously described ( 6 ), supplementing with data from three independent ovarian cancer GWAS. In iCOGS, we excluded samples if they were not of European ancestry, had a genotyping call rate of <95%, showed low or high heterozygosity, were not female or had ambiguous sex, or were duplicates (cryptic or intended). SNPs were excluded if they were mono-morphic, had a call rate of < 95%, showed evidence of deviation from Hardy-Weinberg equilibrium in controls or had low concordance between duplicate pairs( 6 ). For two of the GWAS (from Mayo Clinic and the UK), we also excluded rare SNPs (MAF < 1% or allele count < 5, respectively). The final data set comprised 11,069 cases and 21,722 controls from iCOGS (‘OCAC-iCOGS’), 2,165 cases and 2,564 controls from a GWAS from North America ( 27 ), 1,762 cases and 6,118 controls from a United Kingdom-based GWAS ( 2 ), and 441 cases and 441 controls from the Mayo Clinic. All subjects included in this analysis provided written informed consent as well as data and blood samples in accordance to ethical guidelines under protocols approved by institutional review boards of their respective study sites. Overall, 43 studies from 11 countries provided data on 15,437 women diagnosed with invasive EOC, 9,627 of whom were diagnosed with serous EOC, and 30,845 controls from the general population. We imputed variants separately for the OCAC-iCOGS and each GWAS from 1000 Genomes Project data using the v3 April 2012 release as the reference panel using the IMPUTE2 software ( 28 ) without pre-phasing. The final data set comprised genotypes for 4,234 SNPs of which 2,418 had been directly genotyped. We evaluated the association between genotype and disease using logistic regression by estimating the associations with each additional copy of the minor allele (log-additive models). The analysis was adjusted for study and population substructure by including the eigenvectors of the first five European-specific principal components as covariates in the model. We used the same approach to evaluate SNP associations with serous ovarian cancer after excluding all cases with any other or unknown tumor subtype. For imputed SNPs, we used expected dosages in the logistic regression model to estimate SNP effect sizes and P values. We carried out analyses separately for OCAC-iCOGS and each GWAS and pooled data thereafter using a fixed-effects meta-analysis; thus, all results are based on the combined data. We also performed analyses adjusted for rs3814113 to evaluate evidence of independent signals.

Here, we started from early findings from GWAS for EOC risk and delineated a mechanistic hypothesis for susceptibility at the 9p22.2 locus. Using a two-pronged approach combining functional analysis and fine mapping we identified three genomic features (and enhancer to BNC2 , the BNC2 promoter, and a Substrate/Matrix Attachment Region) harboring twelve potentially functional SNPs ( Fig. 1 ). Based on the likelihood statistics, the most likely causal SNPs were in a Substrate/Matrix Attachment Region (S/MAR) located in a 1.8 Mb topologically associating domain (TAD). Also, this TAD includes associated SNPs revealed in conditional analysis (adjusting for rs3814113) that locate to an enhancer region that interacts with the BNC2 promoter in ovarian cells. Taken together our data implicate multiple candidate causal SNPs at the locus that converge to regulate BNC2 in ovarian cancer susceptibility. Our functional analysis revealed two SNPs in the BNC2 promoter, two SNPs in an enhancer that physically interacts with the BNC2 promoter, and a functional SNP in a S/MAR with allele-specific effects. Of these five SNPs, the strongest genetic evidence for causality is for rs181552334 in the S/MAR. An additional SNP, rs112442786 (r 2 to rs181552334 = 0.9556), located in the same region, which emerged in our fine-mapping approach was not directly tested and may also contribute to risk. S/MARs are thought to help maintain the local 3D chromatin structure by contributing to looping and modulate gene expression ( 34 ). Polymorphisms in S/MARs can regulate, in an allelic-specific manner, attachment to the nuclear scaffold/matrix ( 35 ). Interestingly, our EMSA experiments suggest allelic specific binding to nuclear proteins. Several attempts to remove the S/MAR using CRISPR-based genome editing techniques were not successful suggesting that ovarian cells may not be viable without this S/MAR. CRISPR-based deletion of a region including the S/MAR in 293T HEK cells led to a two-fold reduction of BNC2 expression ( 36 ), with no changes in CNTLN expression, implicating the S/MAR in BNC2 regulation. Notably, all three cis-eQTL associations detected for the 40 credible SNPs from fine mapping were with BNC2 ; and all cis-eQTL associations for the four genes at the locus only BNC2 showed eQTL with the set of credible SNPs. High dimensional chromatin interaction (HiC) data from seven cell lines indicates the presence of a 1.8 Mb topologically associating domain (TAD) in which the S/MAR is situated close to one of its borders. Our data suggest a role for several regulatory interactions, defined by a TAD containing multiple non-coding elements which target BNC2 ( Fig. 1 ). Although future research will further delineate the relationship between BNC2 and ovarian biology, recent reports support our findings. Hnisz et al. ( 37 ) identified a super enhancer in ovarian cells near BNC2 , consistent with BNC2 representing a cell identity gene or master regulator in ovarian cells. The bonaparte zebrafish ( Bnc2) mutants display skin pigmentation defects (no body stripes), stunted growth, and dysmorphic ovaries coupled to infertility ( 38 ). In mice, Bnc2 is expressed in ovarian theca cells, and female mice nullizygous for Bnc2 display an excessive number of stromal cells combined with a reduced number of oocytes ( 39 ). Interestingly, rs12379183, in Region 7, is associated with sonographically detectable abnormalities in the ovaries ( 40 ). Moreover, a network-based integration of GWAS and gene expression in ovarian cancer focusing on transcription factors identified BNC2 using a combination of coexpression and enrichment analysis as a gene contributing to a HOX -centric network associated with serous ovarian cancer risk ( 41 ). Finally, a recent analysis of genetic interactions between germline polymorphisms and tumor formation in specific tissues revealed a significant association between rs3814113 and ovarian cancer ( 42 ). Genome-wide and candidate gene association studies suggest that this locus may also be pleiotropic in humans with effects on ovary, skin, and skeletal biology. SNPs in the 9p22.2 locus have been associated with skin pigmentation in Europeans ( 43 ) and Asians ( 44 ), with freckling ( 45 ), and height ( 46 ).Functional analysis revealed rs12350739 as the likely causal variant contributing to saturation of skin color ( 47 ). SNP rs12350739 was identified in the present study mapping to a Region 8 (a candidate BNC2 enhancer). An introgressed region of Neanderthal DNA (Chr9: 16,720,121–16,786,930) proposed to confer adaptive advantage to colder climates through changes in skin pigmentation is also present at this locus ( 48 , 49 ). Finally, this locus has also been shown to modify ovarian cancer risk in carriers of BRCA1 and BRCA2 pathogenic variants ( 50 ). We acknowledge limitations of this work including the that regulatory networks may be significantly altered during development ( 51 ), the incomplete knowledge of the regulatory landscape in ovarian cells ( e.g. lack of data on CTCF repressor marks and of information on other non-coding RNA elements) and the possibility of missing rare alleles that contribute to the phenotype that could be revealed using the larger Haplotype Reference Consortium data for imputation. Despite these limitations, our data identify plausible and likely biological mechanisms operating to modulate ovarian cancer risk. In summary, we confirmed the region as a highly associated susceptibility locus and propose that the mechanism of ovarian cancer susceptibility at the 9p22.2 locus is likely mediated by changes in a transcriptional regulatory network involving several regulatory elements (enhancers and S/MAR) acting on BNC2.

Epithelial ovarian cancer (EOC) is a poorly understood disease often diagnosed at late stages and with low 5-year survival rates. Although it used to be widely acknowledged that the ovarian surface epithelium (OSE) was the likely tissue of origin of EOC, recent evidence supports the notion that the epithelial lining of the fallopian tube and benign endometriosis contribute to the origin of invasive EOCs. Invasive EOCs may also originate from ectopic Müllerian tissue due to endosalpingiosis. The diverse cellular origins of EOC subtypes, in part, underlie the heterogeneity that characterizes ovarian cancer. Less than half of all familial ovarian cancer cases and less than 15% of high grade serous EOC are due to highly penetrant pathogenic alleles of genes such as BRCA1 and BRCA2 . However, exhaustive family-based linkage studies have not identified additional highly penetrant EOC susceptibility genes ( 1 ). The excess familial risk of EOC may be explained, at least in part, by common variants with low to moderate penetrance. Genome-wide association studies (GWAS) have identified ~40 common variant loci associated with risk of EOC ( 2 – 13 ). Delineation of the mechanisms and likely causal variants at GWAS-identified loci may reveal novel chemoprevention and therapeutic strategies. To evaluate the mechanisms by which single nucleotide polymorphisms (SNPs) may contribute to EOC, we conducted a functional dissection at the 9p22.2 locus, the first ovarian cancer risk locus identified through GWAS of European ancestry women ( 2 ). The SNP most significantly associated with high grade serous EOC risk was rs3814113, which is located 44 kb centromeric and 220 kb telomeric to the BNC2 and CNTLN transcription start sites (TSS), respectively ( 2 ). The minor allele [C; MAF = 0.323] was associated with reduced risk of high grade serous EOC (combined data OR = 0.82; 95%CI = 0.79–0.86; P = 2.5 × 10 −17 ).

Antibody Chromatin immunoprecipitations were performed as previously described 1 using a validated Prestige® BNC2 antibody (Sigma Atlas; cat.no. HPA018525). Prestige Antibodies were developed supported by the Human Protein Atlas ( proteinatlas.org ). According to the manufacturer Prestige Antibodies are tested in a series of validation steps. The BNC2 antibody was able to immunoprecipitate ectopically expressed CBP-tagged BNC2 showing that it can also specifically recognize native BNC2 ( Extended Data Figure 1 ). Lysates of these cells were immunoprecipitated with either Rabbit IgG or the Prestige antibody for BNC2 (Sigma). Immunoprecipitates (IP; #1, 1μg of BNC2 antibody; #2, 2μg of BNC2 antibody) undergo Western Blot for CBP. A band for BNC2 between the 150 kDA and 250 kDA mark appears in the input and BNC2 IP for over expressed BNC2 but not in the input and BNC2 IP for over expressed GFP nor in the IgG IP. ChIP-Seq ChIP-Seq was performed on the endogenous BNC2. In brief, iOSE11 or iFTSEC283 cells at 70% confluence were cross-linked with 1% Formaldehyde in PBS. Crosslinking was quenched by adding Glycine to a concentration of 0.125 M. After washing, cells were collected in Szaks’ RIPA buffer [150 mM NaCl, 1% NP-40, 0.5% deoxycholate, 0.1% SDS, 50 mM Tris HCl pH8, 5 mM EDTA, Protease Inhibitors, 50 mM NaF, 0.2 mM sodium orthovanadate, 0.5 mM PMSF] and the lysate was brought to approximately 1 mg/mL. The lysate was then sonicated in Biogenode Sonicating Water Bath for 12 cycles of 30 sec on and 30 sec off for 8 min. One mg of protein was then mixed with 40 μL of 50% slurry protein A/G agarose beads (Santa Cruz) previously washed in Szaks’ RIPA buffer and pre-cleared for 1–2 h at 4°C. We prepared one lysate per cell line, referred to as OSE_input and FTE_input. Next, pre-cleared lysate was mixed with 5 μg of BNC2 antibody (Sigma Atlas) and 40 μL of 50% slurry protein A/G agarose beads previously washed in Szaks’ RIPA buffer and saturated with 1 mg/mL BSA. The mix was incubated overnight at 4°C while rotating. Beads were then washed twice with Szaks’ RIPA Buffer, four times with Szaks’ IP wash buffer [100 mM Tris HCl pH 8.5, 500 mM LiCl, 1% NP–40, 1% deoxycholate], twice again with Szak’ RIPA Buffer and twice with cold TE. Immunocomplexes were eluted by incubating samples at 65°C for 10 min in 1.5X Talianidis Elution Buffer [70 mM Tris HCl pH 8, 1 mM EDTA, 1.5% SDS]. Crosslinks were reversed by bringing samples to 200 mM NaCl solution and incubating at 65°C for 5 h. DNA was purified by phenol-chloroform extraction and re-suspended in 50 μL 10 mM Tris pH 8.0. For BNC2 ChIP-Seq four individual ChIP samples (from each input lysate) were pooled for each cell line (iOSE11 and iFTSEC283) in two biological replicates, referred to as OSE1, OSE2, FTE1, and FTE2. Immunoprecipitated DNA was used to generate a sequencing library using the NuGEN Ovation Ultralow Library System with indexed adapters (NuGEN, Inc., San Carlos, CA). The library was PCR amplified and size-selected using AxyPrep Fragment Select beads (Corning Life Sciences - Axygen Inc., Union City, CA). The size and quality of the library was evaluated using the Agilent BioAnalyzer, and the library was quantitated with the Kapa Library Quantification Kit (Kapa Biosystems, Woburn MA). Each enriched DNA library was then sequenced on an Illumina HiScan SQ sequencer to generate 100-base paired-end reads. The raw sequence data was de-multiplexed using the Illumina CASAVA 1.8.2 software (Illumina, Inc., San Diego, CA) and binding sites were identified using the MACS2 software 2 using input DNA as a control and callpeak function without building the shifting model, minimum FDR as 0.01. The .bam and .wig files were visualized and inspected using the UCSC genome browser 3 . The number of reads for each sample and their quality metrics are shown in Extended Data Figure 2 . All samples had >70% of reads with Q30 or better and 2% or less of duplicates. For peak calling −log 10 (q value) > 2 (corresponding to an 1% FDR) was used as a cut-off. The number of paired end reads ranged from ~50M to ~69M per sample above the ENCODE minimum requirement of 20M for point-source (ChIP-Seq) experiments 4 . iFTSEC283 cells had a total of 5,687 (FTE1) and 5,730 (FTE2) peaks with 3,396 overlapping peaks and iOSE11 cells had a total of 5,492 (OSE1) and 9,818 (OSE2) with 3,205 overlapping peaks. Peaks used for identification of potential target genes had an intensity greater than 0.05 (reads/length), number of reads greater than 50, and a fold change compared to the input greater than 10 for a total of 2,012 peaks for iFTSEC283 cells and 544 peaks for iOSE11 cells. Median enrichment ranged from 5.2 to 6.9 considered within the norm for ENCODE experiments 4 . Typical peaks are illustrated in Extended Data Figure 3 . A. Sample description and quality metrics. Rawnum and cleanreadsnum, number of raw and clean map reads, respectively. Paired and paired_perc, number and percent of paired mapped reads. ≥Q30(%), percent of map reads Q30 and above. B. Distribution of map reads according to their quality metrics (Q bins). Red dashed line indicates threshold of ≥Q30. C. Peak overlaps between replicates of the same cell lines and overlapped between the two samples. A. FAM49B peak. B. TGBR3 peak. C. Jun peak. Peaks located in these regions and found both in FTE and OSE samples. The total length of the called peak is shown as the blue highlight.