Introduction

Pathogenic mutations in the BRCA1 and BRCA2 genes substantially increase a woman’s lifetime risk of developing breast and ovarian cancers [14]. These risks vary significantly according to (a) age at disease diagnosis in carriers of identical mutations, (b) the cancer site in the individual who led to the family’s ascertainment, (c) the degree of family history of the disease [1, 4, 5], and (d) the type and location of BRCA1 and BRCA2 mutations [6]. These observations suggest that other factors, including lifestyle/hormonal factors [7] as well as other genetic factors, modify cancer risks in BRCA1 and BRCA2 mutation carriers. Direct evidence for such genetic modifiers of risk has been obtained through the association studies performed by the Consortium of Investigators of Modifiers of BRCA1/2 (CIMBA), which have shown that several common breast cancer susceptibility alleles identified through population-based genome-wide association studies (GWASs) are also associated with breast cancer risk among BRCA1 and BRCA2 mutation carriers [810].

Global analysis of GWAS data has shown that the vast majority of common variants associated with susceptibility to cancer lie within genomic non-coding regions and are predicted to account for cancer risk through regulation of gene expression [11, 12]. A recent expression quantitative trait loci (cis-eQTL) analysis for mRNA expression in 149 known cancer risk loci performed in five tumor types (breast, colon, kidney, lung, and prostate) has shown that approximately 30 % of such risk loci were significantly associated with eQTLs present in at least one gene within 500 kb [13]. These results suggest that additional cancer susceptibility loci may be identified through studying genetic variants that affect the regulation of gene expression. In the present study, we selected genes of interest for their known involvement in cancer etiology, identified 320 genetic variants in the vicinity of these genes with evidence of differential allelic expression (DAE), and then investigated the associations of these variants with breast and ovarian cancer risks among BRCA1 and BRCA2 mutation carriers. These included variants in genes involved in DNA repair (homologous recombination and DNA interstrand crosslink repair), interaction with and/or modulation of BRCA1 and BRCA2 cellular functions, cell cycle control, centrosome amplification and interaction with AURKA, apoptosis, ubiquitination, as well as known tumor suppressors, mitotic kinases, and other kinases, sex steroid action, and mammographic density.

Materials and methods

Subjects

All study participants were female carriers of a deleterious germline mutation in either BRCA1 or BRCA2 and aged 18 years or older [14]. Fifty-four collaborating CIMBA studies contributed a total of 23,463 samples (15,252 BRCA1 mutation carriers and 8211 BRCA2 mutation carriers) to this study, including 12,127 with breast cancer (7797 BRCA1 and 4330 BRCA2 carriers) and 3093 with ovarian cancer (2462 BRCA1 and 631 BRCA2 carriers). The number of samples included from each study is provided in Online Resource 1. The recruitment strategies, clinical, demographic, and phenotypic data collected from each participant have been previously reported [14].

Ethics statement

BRCA1 and BRCA2 mutation carriers were recruited through the CIMBA initiative, following approval of the corresponding protocol by the Institutional Review Board or Ethics Committee at each participating center (Online Resource 2); written informed consent was obtained from all study participants [8, 9].

SNP selection and differential allelic expression

SNP selection was performed by first identifying a list of 175 genes of interest involved in cancer-related pathways and/or mechanisms. The list of genes was established by analyzing published results and by using available public databases such as the Kyoto encyclopedia of genes and genomes (http://www.genome.jp/kegg/). Next, DAE SNPs located within these gene regions were identified using previously reported data on allelic expression cis-associations, derived using (1) the lllumina Human1M-duo BeadChip for lymphoblastoid cell lines from Caucasians (CEU population) (n = 53) [15], the Illumina Human 1M Omni-quad for primary skin fibroblasts derived from Caucasian donors (n = 62) [13, 16], and the Illumina Infinium II assay with Human 1.2M Duo custom BeadChip v1 for human primary monocytes (n = 188) [17]. Briefly, 1000 Genomes project data were used as a reference set (release 1000G Phase I v3) for the imputation of genotypes from HapMap individuals. Genotypes were inferred using algorithms implemented in IMPUTE2 [18]. The unrelated fibroblast panel consisted of 31 parent–offspring trios, in which the genotypes of offspring were used to permit accurate phasing. Mapping of each allelic expression trait was carried out by first normalizing allelic expression ratios at each SNP using a polynomial method [19] and then calculating average phased allelic expression scores across annotated transcripts, followed by correlation of these scores to local (transcript ± 500 kb) SNP genotypes in fibroblasts as described earlier [16]. A total of 355 genetic variants were selected on the basis of evidence of association with DAE in the selected 175 genes (see Online Resource 3 for a complete list of SNPs and genes). Following the selection process, SNPs were submitted for design and inclusion on a custom-made Illumina Infinium array (iCOGS) as previously described [8, 9]. Following probe design and post-genotyping quality control, 316 and 317 SNPs were available for association analysis in BRCA1 and BRCA2 mutation carriers, respectively. Genotyping and quality control procedures have been described in detail elsewhere [8, 9].

Statistical analysis

Associations between genotypes and breast and ovarian cancer risks were evaluated within a survival analysis framework, using a one degree-of-freedom score test statistic based on modeling the retrospective likelihood of the observed genotypes conditional on the disease phenotypes [20, 21]. To estimate the magnitude of the associations [hazard ratios (HRs)], we maximized the retrospective likelihood, which was parameterized in terms of the per-allele HR. All analyses were stratified by country of residence and using calendar year and cohort-specific incidence rates of breast and ovarian cancers for mutation carriers. Given 320 tests, the cutoff value for significance after a Bonferroni adjustment for multiple testing was p < 1.5 × 10−4.

The associations between the genotypes and tumor subtypes were evaluated using an extension of the retrospective likelihood approach that models the association with two or more subtypes simultaneously [22].

Imputation was performed separately for BRCA1 and BRCA2 mutation carriers to estimate genotypes for other common variants across a ±50-kb region centered around the 12 most strongly associated SNPs (following the NCBI Build 37 assembly), using the March 2012 release of the 1000 Genomes Project as the reference panel and the IMPUTE v.2.2 software [18]. In all analyses, only SNPs with an imputation accuracy coefficient r 2 >0.30 were considered [8, 9].

Functional annotation

Publicly available genomic data were used to annotate the SNPs most strongly associated with breast cancer risk at locus 11q22.3. The following regulatory features were obtained for breast cell types from ENCODE and NIH Roadmap Epigenomics data through the UCSC Genome Browser: DNase I hypersensitivity sites, chromatin hidden Markov modeling (ChromHMM) states, and histone modifications of epigenetic markers, more specifically commonly used marks associated with enhancers (H3K4Me1 and H3K27Ac) and promoters (H3K4Me3 and H3K9Ac). To identify putative target genes, we examined potential functional chromatin interactions between distal and proximal regulatory transcription factor-binding sites and the promoters at the risk loci, using the chromatin interaction analysis by paired end tag (ChiA-PET) and genome conformation capture (Hi-C, 3C, and 5C) datasets downloaded from GEO and from 4D-genome [23]. Maps of active mammary super-enhancer regions in human mammary epithelial cells (HMECs) were obtained from Hnisz et al. [24]. Enhancer–promoter specific interactions were predicted from the integrated method for predicting enhancer targets (IM-PETs) [25]. RNA-Seq data from ENCODE was used to evaluate the expression of exons across the 11q22.3 locus in MCF7 and HMEC cell lines. For MCF7 and HMEC, alignment files from 19 and 4 expression datasets, respectively, were downloaded from ENCODE using a rest API wrapper (ENCODExplorer R package) [26] in the bam format and processed using metagene R packages [27] to normalize in Reads per Millions aligned and to convert into coverages.

eQTL analyses

The influence of germline genetic variations on gene expression was assessed using a linear regression model, as implemented in the R library eMAP (http://www.bios.unc.edu/~weisun/software.htm). An additive effect was inferred by modeling subjects’ copy number of the rare allele, i.e., 0, 1, or 2 for a given genotype. Only relationships in cis (defined as those for which the SNP is located at <1 Mb upstream or downstream from the center of the transcript) were investigated. The eQTL analyses were performed on both normal and tumor breast tissues (see Online Resource 4 for the list and description of datasets, as well as the sources of genotype and expression data). For all sample sets, the genotyping data were processed as follows: SNPs with call rates <0.95 or minor allele frequencies, MAFs (<0.05) were excluded, as were SNPs out of Hardy–Weinberg equilibrium with P < 10−13. All samples with a call rate <80 % were excluded. Identity by state was computed using the R GenABEL package [28], and samples from closely related individuals whose identity by state was lower than 0.95 were removed. The SNP and sample filtration criteria were applied iteratively until all samples and SNPs met the set thresholds.

Results

From the 175 genes selected for their involvement in cancer-related pathways and/or mechanisms, we identified a set of 355 genetic variants showing evidence of association with DAE (see Online Resource 3 for the complete list of genes and SNPs). Of those, 39 and 38 SNPs were excluded because of low Illumina design scores, low call rates, and/or evidence of deviation from Hardy–Weinberg equilibrium (P value <10−7), for BRCA1 and BRCA2 analyses, respectively. A total of 316 and 317 SNPs (representing 227 independent SNPs with a pairwise r 2 <0.1) were successfully genotyped in 15,252 BRCA1 and 8211 BRCA2 mutation carriers, respectively. Association results for breast and ovarian cancer risks for all SNPs are presented in Online Resource 5.

Breast cancer association analysis

Evidence of association with breast cancer risk (at p < 10−2) was observed for nine SNPs in BRCA1 mutation carriers and three SNPs in BRCA2 mutation carriers (Table 1). The strongest association with breast cancer risk among BRCA1 carriers was observed for rs6589007, located at 11q22.3 in intron 15 of the NPAT gene (p = 4.6 × 10−3) at approximately 54 kb upstream of the ATM gene. Similar associations were observed for two other highly correlated variants (r 2 >0.8) on chromosome 11, namely rs183459 (p = 5.7 × 10−3) also located within NPAT and rs228592 (p = 5.5 × 10−3) located in intron 11 of ATM. No association was observed between SNPs at this locus and breast cancer risk for BRCA2 carriers (Online Resource 5).

Table 1 Associations with breast cancer risk in BRCA1 and BRCA2 mutation carriers for SNPs observed at p < 10−2

The strongest evidence of association with breast cancer risk in BRCA2 mutation carriers was observed for rs6982040, located at 8q11.21 in intron 74 of the PRKDC gene (p = 2.7 × 10−3). However, this variant had a very low frequency in affected and unaffected individuals (MAF values of 0.002 and 0.006, respectively). No association was observed for this locus in BRCA1 carriers (Online Resource 5).

Of the nine SNPs associated with breast cancer risk in BRCA1 mutation carriers, three were primarily associated with estrogen receptor (ER)-negative breast cancer: rs11806633 at 1q42.13 in the CDC42BPA gene (p = 9.0 × 10−3), rs6721310 at 2p23.2 in the BRE gene (p = 3.0 × 10−3), and rs2305354 at 2q11.2 in the REV1 gene (p = 1.0 × 10−3), although the differences between ER-positive and ER-negative disease associations were not statistically significant (Table 2). Of the three BRCA2-associated loci, only rs9468322 at 6p22.1 was associated with ER-positive disease (p = 5.0 × 10−4), although the differences in HRs between ER-positive and ER-negative tumors were not statistically significant (Table 2).

Table 2 Associations with breast cancer risk by tumor subtype in BRCA1 and BRCA2 mutation carriers

Although evidence of association with breast cancer risk was observed for the above-described loci in BRCA1 and BRCA2 mutation carriers, none of these associations reached significance after a Bonferroni adjustment for multiple testing. Imputation using the 1000 Genomes data (encompassing ± 50 kb centered on each of the 12 associated variants, Online Resource 6) identified several SNPs with significant associations in BRCA1 mutation carriers at the 11q22.3 locus (with SNP rs228595 as the most significant, p = 7.38 × 10−6), and which were partly correlated with the genotyped SNPs (r 2 <0.4, Fig. 1). After imputation, we also found associations (albeit not statistically significant after multiple testing adjustments), between one imputed SNP at locus 12p13 (rs2255390, p = 5.0 × 10−4) and breast cancer risk for BRCA1 carriers, and two SNPs and breast cancer risk for BRCA2 carriers, namely 6p22 (chr6:28226644:I, p = 9.0 × 10−4) and 8q11 (rs189286892, p = 2.0 × 10−4).

Fig. 1
figure 1

Manhattan plot depicting the strength of association between breast cancer risk in BRCA1 mutation carriers and all imputed and genotyped SNPs located across the 11q22.3 locus bound by hg19 coordinates chr11:107990104_108173189. Directly genotyped SNPs are represented as triangles and imputed SNPs (r 2 > 0.3, MAF > 0.02) are represented as circles. The linkage disequilibrium (r 2) for the most strongly associated genotyped SNP with each SNP was computed based on subjects of European ancestry that were included in the 1000 Genome Mar 2012 EUR release. Pairwise r 2 values are plotted using a red scale, where white and red means r 2 = 0 and 1, respectively. SNPs are plotted according to their chromosomal position: physical locations are based on the GRCh37/hg19 map. SNP rs228606 was genotyped in the iCOGS array but was not included in our original hypothesis of association with DAE. Gene annotation is based on the NCBI RefSeq gene descriptors from the UCSC genome browser

Ovarian cancer association analyses

Evidence of association with ovarian cancer risk (p < 10−2) was observed for six SNPs in BRCA1 mutation carriers and three SNPs in BRCA2 mutation carriers (Table 3). The strongest association with ovarian cancer risk in BRCA1 carriers was observed for rs12025623 located at 1p36.12 (p = 7 × 10−3) in an intron of the ALPL gene. Another correlated variant (r 2 >0.7) on chromosome 1 was also genotyped, namely rs1767429 (p = 9 × 10−3), which was also located within ALPL. The strongest evidence of association with ovarian cancer risk in BRCA2 mutation carriers was observed for rs2233025 (p = 5 × 10−3), located at 1p32.22 within the MAD2L2 gene. None of these associations remained statistically significant after multiple testing adjustments. Imputed genotypes of SNPs in a region encompassing ± 50 kb centered on each of the nine associated variants did not identify stronger associations.

Table 3 Associations with ovarian cancer risk in BRCA1 and BRCA2 mutation carriers for SNPs observed at p < 10−2

eQTL analysis in breast tissue

To identify the genes influenced via the observed associations with breast cancer at locus 11q22.3, eQTL analysis was performed using gene expression data from tumor and normal breast tissues (for detailed descriptions of datasets, refer to Online Resource 4), and all genotyped as well as imputed SNPs within a 1-Mb region on either side of the most significant genotyped SNP. eQTL associations were observed in both normal and tumor breast tissues in this region, although none of those were correlated with our most significant risk SNPs (Online Resource 7). The strongest eQTL associations were observed in the breast cancer tissue dataset BC241 for the SLC35F2 gene (rs181187590, p = 1.4 × 10−5, r 2 = 0.08, i.e., 8 % of the variation in SLC35F2 expression was attributable to this SNP). Other eQTLs observed in this dataset included ELMOD1 (rs181187590, p = 1.3 × 10−4, r 2 = 0.06), EXPH5 (rs181187590, p = 3 × 10−4, r 2 = 0.054), and ATM (rs4987915, p = 3.7 × 10−4, r 2 = 0.05). In The Cancer Genome Atlas (TCGA) BC765 breast cancer dataset, the strongest associations with gene expression were observed for the non-coding RNA lLOC643923 (rs183293362, p = 2.3 × 10−4, r 2 = 0.02), ATM (rs4987924, p = 8.3 × 10−4, r 2 = 0.015), and KDELC2 (rs4753834, p = 8.6 × 10−4, r 2 = 0.015) loci. The eQTL analysis performed for the TCGA normal breast tissue dataset (NB93) showed an association between SNP chr11:108075271:D and ACAT1 gene expression level (p = 6.5 × 10−3, r 2 = 0.08). No association was observed in the normal breast tissue dataset NB116.

Functional annotation

In order to assess the potential functional role of the most significant risk SNPs in the 11q22.3 region, ENCODE chromatin biological features were evaluated in available breast cells, namely HMECs, breast myoepithelial cells, and MCF7 breast cancer cells. We observed some overlap between features of interest and candidate SNPs within the 11q22.3 region (Fig. 2). The most interesting variant was rs228606, which overlapped a monomethylated H3K4 mark in HMECs. Analysis of data from the Roadmap Epigenomics project also showed overlap with a monomethylated H3K4 mark and with an acetylated H3K9 mark in primary breast myoepithelial cells. From ChiA-PET data, chromosomal interactions were found in the NPAT and ATM genes in MCF7 cells, located mainly in the vicinity of the promoter regions of these genes, which encompassed a strongly associated imputed SNP at this locus, namely chr11:108098459_TAA_T. Lastly, although super-enhancers and predicted enhancer–promoter interactions mapped to the 11q22.3 locus in HMECs, none overlapped with our top candidate SNPs (Fig. 2).

Fig. 2
figure 2

Functional annotation of the 11q22.3 locus. Upper panel functional annotations using data from the ENCODE and NIH Roadmap Epigenomics projects. From top to bottom, epigenetic signals evaluated included DNase clusters in MCF7 cells and HMECs, chromatin state segmentation by hidden Markov model (ChromHMM) in HMECs, breast myoepithelial cells, and variant human mammary epithelial cells (vHMECs), where red represents an active promoter region, orange a strong enhancer, and yellow a poised enhancer (the detailed color scheme of chromatin states is described in the UCSC browser), and histone modifications in MCF7 and HMEC cell lines. All tracks were generated by the UCSC genome browser (hg 19 release). Lower panel long-range chromatin interactions: from top to bottom, ChiA-PET interactions for RNA polymerase II in MCF-7 cells identified through ENCODE and 4D-genome. The ChiA-PET raw data available from the GEO database under the following accession (GSE33664, GSE39495) were processed with the GenomicRanges package. Maps of mammary cell super-enhancer locations as defined in Hnisz et al. [24] are shown in HMECs. Predicted enhancer–promoter determined interactions in HMECs, as defined by the integrated method for predicting enhancer targets (IM-PET), are shown. The annotation was obtained through the Bioconductor annotation package TxDb.Hsapiens.UCSC.hg19.knownGene. The tracks have been generated using ggplot2 and ggbio library in R

Discussion

DAE is a common phenomenon in human genes, which represents a new approach to identifying cis-acting mechanisms of gene regulation. It offers a new avenue for the study of GWAS variants significantly associated with various diseases/traits. Indeed, the majority of GWAS hits localize outside known protein-coding regions [11, 12], suggesting a regulatory role for these variants. In the present study, we have assessed the association between 320 SNPs associated with DAE and breast/ovarian cancer risk among BRCA1 and BRCA2 mutation carriers. Using this approach, we found evidence of association for a region at 11q22.3, with breast cancer risk in BRCA1 mutation carriers. Analysis of imputed SNPs across a 185-kb region (±50 kb from the center of each of the three genotyped SNPs at this locus) revealed a set of five strongly correlated SNPs that were significantly associated with breast cancer risk. This region contains several genes including ACAT1, NPAT, and ATM. ACAT1 (acetyl-CoA acetyltransferase 1) encodes a mitochondrial enzyme that catalyzes the reversible formation of acetoacetyl-CoA from two molecules of acetyl-CoA. Defects in this gene are associated with ketothiolase deficiency, an inborn error of isoleucine catabolism [29]. NPAT (nuclear protein, co-activator of histone transcription) is required for progression through the G1 and S phases of the cell cycle, for S phase entry [30], and for the activation of the transcription of histones H2A, H2B, H3, and H4 [31]. NPAT germline mutations have been associated with Hodgkin lymphoma [32]. Finally, ATM (ataxia telangiectasia mutated) encodes an important cell cycle checkpoint kinase that is required for cellular response to DNA damage and for genome stability. Mutations in this gene are associated with ataxia telangiectasia, an autosomal recessive disorder [33]. ATM is also an intermediate-risk breast cancer susceptibility gene, with rare heterozygous variants being associated with increased risk of developing the disease [34]. Although several studies have assessed the role of the most common ATM variants in breast cancer susceptibility, the results obtained are inconsistent [35]. A recent study had identified an association between an ATM haplotype and breast cancer risk in BRCA1 mutation carriers with a false discovery rate-adjusted p value of 0.029 for overall association of the haplotype [36]. Four of the five SNPs making up the haplotype were almost perfectly correlated (r 2 >0.9) with the three originally genotyped SNPs of the present study. These SNPs were, however, only moderately correlated (r 2 >0.4) with the most significant risk SNPs (p = 10−6), identified later through imputation.

Although eQTL analysis has identified cis-eQTL associations between several variants and ACAT1, ATM as well as other neighboring genes in both breast carcinoma and normal breast tissues, none of these associations involved the most significantly associated risk SNPs. Furthermore, the correlation between eQTLs and the most significant risk SNPs was weak. The lack of consistency between the eQTL results and the allelic imbalance data originally used for SNP selection in the design of the present study can probably be explained by the differences between the cell types used in these analyses. The list of allelic imbalance-associated SNPs was selected from studies performed in lymphoblastoid cell lines [15], primary skin fibroblasts [13, 16], and primary monocytes [17], while eQTLs were analyzed in breast carcinoma and normal breast tissue. This tissue heterogeneity in the data sources used represents one of the limitations of this study, although no such data were available in mammary cells when this study was originally designed.

The identification of a region at 11q22.3 that is associated specifically with breast cancer risk in BRCA1 mutation carriers may explain why association studies performed using breast cancer cases from the general population have so far yielded conflicting results with regard to common variants at this locus. The majority of tumors arising in BRCA1 carriers show either low or absent ER expression, while the majority of BRCA2-associated tumors are ER positive, as in most sporadic cancers arising in the general population. Large-scale studies using only ER-negative or triple-negative (i.e., ER-, progesterone receptor-, and HER2-negative) cases could therefore be helpful to confirm the association of this locus with breast cancer risk.