Introduction

Breast cancer is globally the leading cause of cancer-related mortality for women1. The average 5-year survival rate is 83–90% in the Western countries, substantially better than for many other cancers, but due to its high incidence, breast cancer still leads the statistics for cancer mortality in Europe and comes second in North America1,2. On an individual level, prognosis varies greatly, depending on both the inherent tumor biology and the stage of malignant progression at diagnosis. Women with early-stage, localized disease have a very good 5-year prognosis of ~97–99%, but for 10–15% of women, diagnosed with locally advanced disease, the expected 5- and 10-year survival rates range between 40–80% and 30–40%, respectively. Furthermore, the mortality associated with metastatic breast cancer is even higher, with median survival of <3 years2,3,4.

Currently, the prognosis of breast cancer patients is based on the tumor characteristics. Gene expression and copy number profiling5,6,7 or expression of marker proteins, like estrogen receptor (ER), progesterone receptor (PgR), human epidermal growth factor receptor 2, and marker of proliferation Ki-67, can be used to categorize breast cancers into biological subtypes with specific treatment recommendations and different estimates for patient survival8,9. Tumor grade is a histological measure of cellular growth pattern and the best individual prognostic factor10. The phase of tumor progression is approximated by clinical stage, which is based on the tumor size, as well as local and systemic spread of metastatic cells11. However, there is great variation in the survival rates between tumors with similar characteristics and stage. This variance has been suggested to have a heritable component, possibly consisting of genetic differences in metastatic potential and sensitivity to adjuvant therapy12,13,14,15,16,17,18. In addition, host factors, like tumor–microenvironment interaction, immune surveillance, and efficiency in drug metabolism may account for the genetic variability in breast cancer patient survival19,20,21.

Genetic determinants of patient prognosis and the treatment outcome prediction have been intensively sought using both candidate gene and genome-wide approaches, but only some of the discoveries have been successfully validated in subsequent studies. A meta-analysis of ten genome-wide association studies nominally validated 12 out of 45 earlier discoveries for survival from ER-positive or any breast cancer22. Most recently, the Breast Cancer Association Consortium (BCAC) reported two loci from chromosome 7, based on a well-powered meta-analysis of genome-wide association studies23.

In this study, we focused on the survival of women who carry germline pathogenic BRCA1 or BRCA2 variants. BRCA1 and BRCA2 are the two most important breast cancer susceptibility genes, with ~72 and 69% life-time risk of breast cancer and 44 and 17% risk of ovarian cancer, respectively24. BRCA1/2-deficient tumors have distinctive genomic aberration profiles, characterized by homologous recombination deficiency25, making them stand out as etiologically and phenotypically coherent groups of breast carcinomas. BRCA1 risk variants are associated with high grade and triple-negative breast cancer, which both predict poor outcome. BRCA2 risk variants predispose primarily to ER-positive breast cancer, but the risk of ER-negative breast cancer increases with age26,27. Meta-analyses on the survival of women with BRCA1/2 variants and breast cancer have suggested no significant difference in comparison to noncarriers with phenotypically similar tumors28,29,30. However, a recent study suggested ER-positive breast cancer as an adverse indicator for cases with BRCA2 variants, unlike for noncarriers27.

Results

Tumor characteristics

We investigated genetic survival associations in pathogenic variant carriers from the Consortium of Investigators of Modifiers of BRCA1/2 (CIMBA), genotyped on the OncoArray31,32. Twenty-one independent studies participating in CIMBA had survival data available for germline carriers of pathogenic BRCA1 variants (n = 3008) and 15 studies for carriers of BRCA2 variants (n = 2,009; Supplementary Table 1). Primarily, we analyzed patient survival in relation to all-cause death, because this was most complete across the participating studies. As a sensitivity analysis, the discovered survival variants were also always assessed for breast cancer-specific death.

Data on tumor characteristics was available for about two thirds of patients (Table 1). The distribution of the tumor characteristics of BRCA1 and BRCA2 carriers agreed with previous reports, so that small, high-grade, and early-onset tumors had a relatively high frequency26,33. Tumors from BRCA1 carriers were mostly ER-negative, while those from BRCA2 carriers were largely ER-positive. Therefore, the main survival analyses were performed in parallel in the following patient groups: all BRCA1 carriers, BRCA1 carriers with ER-negative tumors, all BRCA2 carriers, and BRCA2 carriers with ER-positive tumors.

Table 1 Tumor characteristics of the patients included in the survival analysis.

rs57025206 predicts survival of BRCA1 carriers with ER-negative breast cancer

We analyzed the association of germline genetic variants with the overall survival after the first primary breast cancer diagnosis in carriers of pathogenic BRCA1 or BRCA2 variants, using Cox regression. The analyses were adjusted for age at breast cancer diagnosis and stratified by country group to account for underlying genetic and clinical differences between populations. A single variant, rs57025206 (minor allele frequency in European population, MAFEUR, 0.027; info score 0.97 for imputed variant), a nine base pair insertion in an intergenic region of 3p21.2, was identified as a highly significant survival marker for ER-negative breast cancer patients carrying BRCA1 pathogenic variants, with a hazard ratio (HR) of 4.37 (95% confidence interval (CI) 3.03–6.30, P = 3.1 × 10−9, Fig. 1). Furthermore, a multivariable analysis adjusted for tumor grade, size, PgR expression status, and lymph node involvement, suggested that the association is independent of these tumor characteristics, with rs57025206-associated HR: 6.19, 95% CI: 3.73–10.3. A similar trend was observed in the analysis of breast cancer-specific death, although the number of patients with available data was much smaller (Supplementary Table 2).

Fig. 1: Survival variant rs57025206.
figure 1

a Kaplan–Meier plot on the survival of ER-negative breast cancer patients carrying pathogenic BRCA1 variants, stratified by rs57025206. b Forest plot of HR associated with rs57025206 across country groups. (P.het: P-value against between-study heterogeneity). c Kaplan–Meier plot on the survival of ER-positive breast cancer patients carrying pathogenic BRCA1 variants, stratified by rs57025206.

The survival association was specific to women with ER-negative breast cancer, since in the analysis of all BRCA1 carriers the HR was attenuated (HR: 2.12, 95% CI: 1.55–2.91, Supplementary Fig. 1) and not seen in those with ER-positive breast cancer (n: 476, HR: 0.52, 95% CI: 0.17–1.60, Fig. 1c). In the analysis of BRCA2 carriers, rs57025206 was not associated with patient survival (HR: 0.94, 95% CI: 0.59–1.51).

The main analyses suggested nine further survival loci

We used a looser discovery threshold P-value, P < 5 × 10−7, than in the traditional genome-wide setting to characterize loci that potentially modify the patient survival and may form a basis for future research hypotheses. In the analysis of the overall survival of BRCA1 carriers, six variants exceeded the selected threshold (Table 2 and Supplementary Fig. 2). The risk magnitude was proportionally associated with the number of risk alleles for two variants, rs59010985 and rs537497819, whereas the other four variants followed a dominant inheritance pattern. The analysis of BRCA1 carriers with ER-negative breast cancer identified four further potential variants influencing survival (Table 2 includes also rs57025206, presented above). No variant exceeded the significance threshold in the analysis of BRCA2 carriers, alone. However, a meta-analysis of BRCA1 and BRCA2 carriers (see below) suggested germline variants associated with survival in both of these groups.

Table 2 Variants associated with the survival of breast cancer patients carrying pathogenic germline variants in BRCA1.

The reliability of the findings was estimated with a Bayesian false discovery probability (BFDP) and by assessing the between-strata (country group) heterogeneity. Ten out of the 11 BRCA1-associated survival variants had BFDP ≤ 0.33 for at least one of the three tested prior maximum effect sizes (Table 2). For these ten variants, there was little heterogeneity between countries (Supplementary Fig. 3, P against heterogeneity >0.1 for all variants). Moreover, the variant effect sizes were consistent in multivariable models adjusted for tumor pathologic characteristics and in the analyses of breast cancer-associated death (Supplementary Table 2), suggesting a consistent effect throughout the data.

Most of the identified survival variants complied with the proportional hazards assumption, suggesting a constant HR over time, but rs537497819 was especially associated with poor short-term prognosis, the effect leveling out after the first 5 years following the diagnosis (Table 2).

Meta-analysis suggested four variants with the consistent survival effect in BRCA1 and BRCA2 carriers

The correlation between variant effect sizes in the BRCA1 and BRCA2 survival analyses was quite low (Pearson’s R = 0.0065 [95% CI 0.0059–0.0072]), suggesting no overall trend for similar genetic survival effects in these patient groups. Nevertheless, a fixed-effects meta-analysis highlighted four candidate variants, which had consistent effect on patient survival in both groups (Table 3, Supplementary Table 3 and Supplementary Fig. 4).

Table 3 Variants associated with the survival of breast cancer patients in the meta-analysis of BRCA1 and BRCA2 carriers.

Four variants may have age-dependent survival effects for BRCA2 carriers

The lack of positive findings in the BRCA2-specific analyses prompted us to consider the possibility that there may be confounding factors, so that the genetic survival associations would depend on tumor and patient characteristics. We had tumor pathology data available only for a subgroup of the study subjects, and lacked statistical power to investigate the interaction between genetic variants and tumor characteristics. However, the age at diagnosis was available for all cases, and to investigate potential age-dependent genetic survival effects, we included in the Cox regression model a covariate coded as the product of the variant genotype and diagnosis age (continuous). The model containing the interaction term was tested against a nested model without interaction. Likelihood-ratio test P-values (Table 4) were corrected for multiple testing using the Benjamini–Hochberg method and variants with false discovery rate (FDR) ≤ 0.33 were considered as potential discoveries. Four variants, all from the analysis of BRCA2 carriers, passed the threshold. For illustrative purposes, the HRs in Table 4 are presented separately for two age groups, even though the regression model suggested that the interaction between the variant survival effect and the diagnosis age is continuous (Supplementary Fig. 5). The variant effect sizes were consistent in multivariable models, including tumor pathologic characteristics, as well as in the analysis of breast cancer-specific death (Supplementary Table 4).

Table 4 Variants, whose association on patient survival was dependent on the age at the breast cancer diagnosis.

ER-negative breast cancer is more common for postmenopausal than premenopausal BRCA2 carriers26,33. To test whether the age-dependent survival association of the four variants (Table 4), was a hidden association with the tumor subtype, we analyzed the survival separately in age- and subtype-stratified subgroups, including only patients with data on ER expression available (Table 1). This analysis suggested that the ER status did not explain the age-dependent survival association. For two variants (rs2109815 and rs35431863), the age-dependent survival association was similar in ER-positive and ER-negative patient groups. For the other two variants (rs372812554 and rs11255420), the survival association did not vary by age of diagnosis in the ER-negative patient group, and the group resembled the younger ER-positive patient group (Supplementary Table 5).

The positional and eQTL analyses reveal potential target genes

None of the variants we found associated with survival was located on a protein-coding region. Therefore, we took two parallel approaches, positional and gene expression based, to identify potential target genes. A target gene was considered to have a strong positional evidence, if the survival variant or any regulatory variant in linkage disequilibrium was located on an experimentally validated enhancer of the target gene in a relevant tissue (mammary, ovarian, blood, and adipose)34, based on the 1000 genomes and Encode data35,36. This is indicated in the Table 5, column “Positional target tissue”. If no other evidence was available, the target gene was assumed the closest gene (no “Positional target tissue” mentioned in Table 5). For the gene expression-based approach, we used eQTL data from public databases covering various human tissues and cell types (Supplementary Data 1 and 2).

Table 5 Target genes.

rs57025206 and its proxy variants in high linkage disequilibrium (R2 > 0.8 in the European population, green track named “lead variant and proxies with R2 > 0.8” in Fig. 2) are co-located with a cluster of IQCF-genes (IQCF1–6), expressed exclusively in the testis. These variants have eQTL associations with altogether five genes in different tissues, for example, with MIR135A in adipose tissue. However, the region does not contain annotated regulatory elements. Furthermore, it is not in physical contact with any transcription start site (TSS) in human mammary epithelial cells (HMEC), suggesting that this region does not contain any HMEC-specific enhancer (see, for example, position “g” in the heatmap of Fig. 2). Furthermore, the MAFs of these variants range from 2.2 to 3.2% in the European population, but between 12 and 27% globally, questioning whether the causal variant would be rs57025206 or any of the proxies with R2 > 0.8.

Fig. 2: Potential target genes on a survival locus 3p21.
figure 2

The heatmap visualizes the strength of association between two genomic loci in chromosome conformation capture data from human mammary epithelial cells. The graph has been created and the topologically associated domains (TAD) estimated with 3D genome browser, promoter.bx.psu.edu/hi-c/. The UCSC browser tracks visualize genomic annotation for the region. The three custom tracks indicate the positions of variants in linkage disequilibrium with rs57025206. The vertical dashed lines indicate the positions of the lead variant rs57025206 and three linked regulatory variants. UCSC genes track visualizes the positions of exons and introns for all genes on the region. The GeneHancer track displays the interactions of distant enhancers and gene transcription start sites (TSSs). The histone mark track shows the positions of H3K27Ac from the Encode project. In the heatmap, seven positions (ae) have been marked in the heatmap to indicate the probability of physical interaction between variants linked with rs57025206 and their target genes based on positional and eQTL analysis (Supplementary Data 1 and 2) as follows: (a) RAD54L2 TSS and rs72945708 on enhancer GH03J051383 targeting RAD54L2. (b) GRM2 TSS and GRM2 eQTL variant rs72945708 on enhancer GH03J051383. (c) DCAF1 TSS and rs9824779 on enhancer GH03J051668 targeting DCAF1. (d) DCAF1 TSS and rs56942057 on enhancer GH03J051684 targeting DCAF1. (e) BAP1 TSS and BAP1 eQTL variant rs9824779 on enhancer GH03J051668. (f) BAP1 TSS and BAP1 eQTL variant rs60497133. (g) MIR135A TSS and MIR135A eQTL variants rs1476290 and rs1605067.

The rs57025206 haploblock (with D′ > 0.8) covers a larger genomic region with functional variants located on enhancers regulating several protein-coding genes (turquoise and green tracks in Fig. 2, Table 5 and Supplementary Data 1). For example, rs56942057 (MAFEUR 0.1%, MAFGLOBAL 4.6%, DEUR 1.0, and DGLOBAL 0.34) and rs9824779 (MAFEUR 0.1%, MAFGLOBAL 3.9%, DEUR 1.0, and DGLOBAL 0.32) are located in enhancers of DCAF1, GH03J051684, and GH03J051668, respectively. Both variants affect conserved sequences of transcription factor binding motifs. Furthermore, the variant positions and DCAF1 TSS are in physical contact in HMEC (positions “c” and “d” in Fig. 2). These two (rs56942057 and rs9824779) and rs60497133 were eQTL variants for BAP1 in thyroid (Supplementary Data 2), but in HMEC there was no contact between the variant positions and BAP1 TSS (Fig. 2, positions “e” and “f”), suggesting that the regulatory associations of these variants and BAP1 may not exists in mammary epithelial cells. A third potentially causal variant, rs72945708 (MAFEUR 2.0%, MAFGLOBAL 18%, DEUR 0.95, and DGLOBAL 0.71), is located on a RAD54L2 enhancer GH03J051383 and associated with the expression of GRM2 in adipose tissue. However, this variant does not change conserved sequence of any known transcription factor binding motif, and the topological association between the variant locus and RAD54L2 or GRM2 TSS cannot be seen in HMEC (positions “a” and “b” in Fig. 2, respectively). Thus, the plausible culprit gene of the rs57025206 locus is DCAF1, an E3 ubiquitin ligase substrate receptor targeting, e.g., TP53, ER-alpha, and EZH2.

We performed a similar target gene analysis for the 17 further survival candidate variants (Tables 24) using positional and eQTL data (Table 5 and Supplementary Data 1 and 2). For only one variant, rs59010985, there was a single clear target gene, GAS7, which was supported by both analyses. The majority of the other variants could be divided into three categories. Some variants had multiple potential target genes. For example, variants in linkage disequilibrium with rs4879914, detected in the BRCA1BRCA2 meta-analysis, were eQTLs for several genes in white blood cells. Positional evidence supported the same array of target genes, in both mammary epithelial cells and white blood cells. For some other loci, like ASPH and CREB5, the target gene prediction was based merely on the regulatory variants located on active enhancers in relevant tissues, with no eQTL evidence. Furthermore, for six loci, the target gene was assumed the closest gene. One of the 17 survival variants was located in an intergenic region with no apparent target gene.

To detect common pathways or recurrent cellular or molecular functions, we performed a systematic literature review (Fig. 3 and Supplementary Table 6). Relevant literature was available for 26 out of the 40 target genes listed in Table 5, when searched from PubMed with the keyword combinations described in the “Methods” section. Twelve of the genes had previously been connected with breast cancer patient survival, based on either mRNA or protein expression, germline genetic variation, or copy number change in mammary tumors (Fig. 3 and Supplementary Table 6). Fourteen of the target genes were associated with proliferative and migratory capacity of mammary or other epithelial cells. This association was mediated via two primary mechanisms: regulation of MAPK/ERK pathway activity or modulation of cytoskeletal proteins. Target genes in two out of the four loci associated with the survival of BRCA1 carriers with ER-negative breast cancer (DCAF1, ZRANB1, and CTBP2) and two genes from the BRCA1BRCA2 meta-analysis (ZNF644 and TFCP2L1) affect the regulation of chromatin state either directly or via polycomb repressor complex 2 (PCR2). Eight target genes have been suggested to affect the response to adjuvant chemotherapy, whereas three target genes have been associated with the response to adjuvant endocrine therapy. The latter included two of the four loci (ASPH and GATA3) associated with age-dependent survival effect in BRCA2 carriers. Seven genes from six loci had been suggested to be involved in the regulation of immune response. However, four of these genes were also associated with regulation of mammary cell differentiation.

Fig. 3: Functional clustering of the target genes based on published literature.
figure 3

For details and references, see Supplementary Table 6. PRC2 polycomb repressor complex 2, BC breast cancer, CNA copy number aberration.

Comparison of the results to survival associations from a general breast cancer population

None of the survival associations discovered in the BRCA1/2 carriers were detected in the BCAC data of unselected and familial non-BRCA1/2 breast cancer patients, when our results were compared to those of Escala-Garcia et al. (Supplementary Table 7)23.

Discussion

We studied germline genetic variants associated with the survival of breast cancer patients carrying pathogenic variants in the high-risk BRCA1 and BRCA2 genes. We identified one variant, rs57025206, which was highly significantly associated with poor survival of BRCA1 carriers with ER-negative breast cancer, with HR = 4.37, 95% CI 3.03–6.30, P = 3.1 × 10−9. Furthermore, we discovered 17 additional candidate loci, which could enhance the understanding of the biological processes modifying the survival of breast cancer patients and form a basis for future research. It was noteworthy that the single significant discovery was made in a subgroup of women with BRCA1 pathogenic variants and ER-negative breast cancer, and this was attenuated in the wider analysis of all BRCA1 carriers, because inclusion of the ER-positive breast cancer cases, where the association did not exist, diluted the effect in the combined analysis. Thus, the strength of the discovery analysis came from a phenotypically homogeneous group of cases with shared disease etiology, rather than a large number of study subjects.

Minor alleles of two rare variants, rs9824779 and rs56942057, are in complete linkage (D′ = 1.0, R2 = 0.036) with the rare, poor-survival, allele of rs57025206. These two variants affect conserved transcription factor binding sequences in enhancers targeting the DCAF1 TSS, and the enhancer–TSS interaction is present in HMEC (Fig. 2). Based on our analysis, these variants would be the best functional candidates causing the association signal, although their frequency was too low to allow survival analysis in the CIMBA data.

DCAF1 is an ubiquitin ligase substrate receptor, which recognizes and binds substrates for the CLR4 and HECT-type ubiquitin ligases, thus regulating the substrate half-life. Furthermore, the casein kinase-like and Lis1 homology domains of the DCAF1 protein have been shown to modulate histones by phosphorylation and deacetylation, respectively37,38. According to in vitro models, DCAF1 silencing reduced proliferation and colony formation of MCF7 breast cancer and DU145 prostate cancer cell lines39,40, whereas high DCAF1 expression blocked the expression of tumor suppressor genes via H2A phosphorylation, and induced a proliferative gene expression signature together with FOXM1 (refs 38,40). In mammary epithelial progenitor cells, DCAF1 is involved in the regulation of the Hippo pathway and ER-alpha, contributing to the maturation of luminal and basal lineages41.

The genomic locus of rs57025206, 3p21, is frequently deleted in mammary carcinomas, especially in ER-negative and high-grade tumors42,43,44. Subsequent research has not been able to name an unequivocal cancer driver gene, even though several candidates have been suggested, including BAP1 and MIR135A45,46, which were detected as eQTL variants in our target gene analysis (Supplementary Data 2). However, the topological data suggested that these eQTLs might represent tissue-specific regulatory associations not present in HMEC, and that if the rs57025206-related survival effect is mediated via BAP1 or MIR135A, it probably reflects host–tumor interactions.

To explore the molecular mechanisms possibly contributing to the survival of breast cancer patients with pathological BRCA1 and BRCA2 variants, we included in a target gene analysis and functional literature review 17 sub-genome-wide significant variants with reasonably low false discovery probability. Four variants came from the BRCA1BRCA2 meta-analysis, four variants had age-dependent survival association for BRCA2 carriers, and the remaining nine variants were associated with the survival of BRCA1 carriers (Tables 24). With the selected thresholds for reporting, five to six loci are expected to be false discoveries, while the rest are likely to be true survival loci. Significant between-strata heterogeneity was not detected for any tested variant.

A literature review of the plausible target genes indicated that similar biological functions underscored the genetic survival associations of BRCA1 and BRCA2 carrier breast cancer patients (Supplementary Table 6 and Fig. 3). The most prevalent function was the regulation of epithelial cell proliferative and migratory capacities, suggesting that germline variants modifying the liability of phenotypic changes in epithelial cells could be key determinants for the outcome of breast cancer. Another interesting observation was the accumulation of genes regulating mammary epithelial cell differentiation and stemness on the loci associated with the survival of BRCA1 carriers with ER-negative breast cancer. In addition to DCAF1, described above, these candidate genes include ZRANB1 and CTBP2 located on 10q26.13. ZRANB1 is a deubiquitinase targeting, e.g., EZH2, a component of the PRC2 and direct regulator of BRCA1 activation47,48, whereas CTBP2 is a transcriptional corepressor priming target gene promoters, including BRCA1, for PRC2-dependent silencing49,50.

The functional review highlighted treatment response as another important mechanism contributing to the survival of BRCA1/2 pathogenic variant carriers. The GATA3 locus from the analysis of BRCA2 carriers, the TPM2/GBA2 locus from the BRCA1–BRCA2 meta-analysis, and five out of ten loci from the analyses of BRCA1 carriers have previously been reported to affect the response to adjuvant chemotherapy (Supplementary Table 6 and Fig. 3). Interestingly, KIF26B, GAS7, MIR135A, and CTBP2, all potential target genes of variants associated with the survival of BRCA1 carriers, have been suggested to modify the response to platinum-based chemotherapy51,52,53,54,55,56. Furthermore, CTBP2 has been shown to affect the sensitivity to PARP inhibitors by targeting the BRCA1 promoter for epigenetic silencing50. Platinum compounds and PARP inhibitors are effective agents for adjuvant therapy of BRCA1 carriers57,58. Therefore, the information on genetic variation in these potential response-modifier loci could add to design of targeted therapy trials.

In the analysis of BRCA2 carriers, we did not find any variants exceeding the discovery threshold. This may be partly explained by the smaller number of BRCA2 carriers, and consequently events, in comparison to BRCA1 carrier analysis. However, it was intriguing that with similar thresholds for P-value and false discovery probability, we discovered four variants with age-dependent survival effects. The majority of the BRCA2-related breast cancers are hormone receptor positive (Table 1), and rely on the endogenous supply of estrogen and possibly progesterone for proliferation59. The level of both of these steroid hormones decreases between the ages of 35 and 55 years60, consistent with the diagnosis age of the majority of the BRCA2 carrier cancers in these data (Table 1). The target gene functional review indicated that two loci with age-dependent survival effect in BRCA2 carriers, GATA3 and ASPH, are associated with differences in response to adjuvant endocrine treatment61,62, further supporting the hypothesis that in hormone receptor-positive breast cancer, the age-related differences in hormone secretion and sensitivity may affect the survival of the patients.

We can draw two conclusions from our analyses, which support the hypothesis that the tumor etiology and subtype influence the genetic survival associations. First, the strongest BRCA1 survival variants were not associated with the survival of BRCA2 variant carriers. Second, we did not find any corroboration for our discoveries from the BCAC data (Supplementary Table 7). BRCA1 and BRCA2 carrier tumors have characteristic mutation profiles, accompanied by homologous recombination deficiency, which distinguishes these from unselected breast carcinomas on a molecular level25. Notwithstanding, the mutations in these two genes are associated with different breast cancer subtypes26,27. Taken together, our results encourage accuracy in definition of the outcome of interest and stringency in collecting comparable study subjects, in order to improve the future survival studies.

Like many earlier genetic survival association studies, this study had a limited discovery power, despite the fact that the analyzed high-risk variant carrier cohorts were enriched with poor-prognosis cancers, and therefore contained more events than any equally sized cohort of unselected patients would contain. With two phenotypically coherent groups of breast cancer patients, which represent the largest collections of BRCA1 and BRCA2 pathologic variant carriers in the world, we were able to identify a single genome-wide significant locus, but were probably underpowered to detect loci with more modest effect sizes. We estimate that better coverage of treatment and pathology would have greatly improved our analyses, especially as the functional literature review suggested that the survival differences could be mediated by differential response to adjuvant treatment. Unfortunately, we did not have enough well-documented treatment data available to test the hypothesis. The retrospective nature of our data brought its own limitations. The data included a notable proportion of prevalent cases, enrolled into the participating studies after a prolonged time period after the initial breast cancer diagnosis. The choice of the overall survival as the basis of the analyses impairs the clinical interpretation of the results to some degree. However, most of the follow-up information came from registries, and therefore the all-cause death was the best available indicator for poor prognosis. We tried to overcome these shortcomings with appropriate analysis methods for the prevalent data, with posterior likelihood estimation and sensitivity analyses with breast cancer-specific death as an end point, as well as by testing for internal consistency.

In conclusion, we report a survival locus for BRCA1 carriers with ER-negative breast cancer. Furthermore, the results from the exploratory analyses suggest that the survival in women with breast cancer is influenced by a complex action between clinical characteristics, tumor biology, and the germline genetic landscape.

Methods

Study subjects

The study subjects included women of European ancestry diagnosed with invasive breast cancer before the age of 70 years, enrolled in studies participating in the CIMBA (Supplementary Table 1). A CIMBA study was included in the analyses if it provided sufficient amount of follow-up data, defined as at least 15 study subjects at risk during the time when five events occurred. The study-wise inclusion criterion was applied separately for the main analyses and the ER-specific subgroup analyses. This selection yielded survival data on 3008 women carrying pathogenic germline variants in BRCA1 from 21 studies and 2009 women carrying variants in BRCA2 from 15 studies. Tumor characteristics of the study subjects are provided in Table 1. The BRCA1 carriers were collected from 2664 families: 2391 families with one study subject, 220 families with two study subjects, and 53 families with more than two study subjects. The BRCA2 carriers were collected from 1713 families: 1486 families with one study subject, 182 families with two study subjects, and 49 families with more than two study subjects. All participating studies were approved by their appropriate ethics review boards, and all subjects provided informed consent.

Genotype data

The study subjects were genotyped with a custom-made array as a part of the OncoArray project31. Details on the variant selection and data quality control have been published elsewhere. In short, the genotyping array included a GWAS backbone (Illumina HumanCore) and potentially cancer-associated variants nominated by the six participating consortia. Genotyping of the CIMBA samples was conducted in six independent laboratories, which used the same HapMap reference samples and a common genotype-clustering file to ensure interlaboratory comparability31. The data were imputed using the 1000 genomes as a reference panel, as described previously32. In the analyses, we included variants with at least 60 carriers, corresponding to MAFs 0.01 and 0.015 for BRCA1 and BRCA2 carriers, respectively. This yielded data on ~9.7 million SNPs for BRCA1 carriers and 9.1 million SNPs for BRCA2 carriers.

Survival analysis

The patients were followed from the diagnosis of the first primary breast cancer until death of any cause and censored after 15 years or when lost from follow-up. Left truncation was applied to account for delayed study entry. In the BRCA1 carrier analysis, the number of person years was 16,056, the number of deaths 461, and the 15-year survival rate 0.66. The maximum number of study subjects at risk was 1336 and this was reached 2.8 years after the baseline. The BRCA2 carrier dataset covered 10,712 person years with 311 deaths leading to a 15-year survival rate of 0.65. Here, the highest number of study subjects at risk, 930, was reached 3.9 years after the baseline.

The genome-wide analysis of association between genetic variants and all-cause mortality was performed with Cox regression as implemented in the survival library of the R environment for statistical computing version 3.5.1 and 3.6.0 (refs 63,64,65,66). Linear per-allele risk model and dominant inheritance model were estimated in parallel. Analyses were adjusted for diagnosis age, allowing for variant–age interaction, and stratified by country group to account for population-specific genetic and clinical features (Supplementary Table 1). Likelihood-ratio test was used as a measure of significant association with alpha risk P < 5 × 10−7 (two-sided). The interaction model, i.e., Cox regression model including an interaction term, coded as a product of the number of the effect alleles (0, 1, 2) and diagnosis age (continuous), was tested against a nested model containing the variant and age without interaction. Robust variance estimation was used to account for relatedness of study subjects from the same families. Parallel genome-wide analyses were performed also within the biologically homogeneous patient groups: BRCA1 carriers with ER-negative tumors (1385 patients from 17 studies, see above the study-wise inclusion criterion for study-stratified analysis) and BRCA2 carriers with ER-positive tumors (1050 patients from 14 studies). The results from analysis of all BRCA1 were compared to results from analysis of all BRCA2 carriers with Pearson correlation, and combined using a fixed-effects meta-analysis as implemented in R library metafor67. For the meta-analysis, the standard errors were recalculated using the likelihood-ratio test statistic to avoid inflation caused by rare variants, as suggested previously23.

Using R-library powerSurvEpi68, we estimated that we had sufficient statistical power, with the selected alpha-risk (5 × 10−7) and beta-risk 0.2, to detect significant risk associated with common variants (MAF > 0.10), if the HR was >1.8, whereas for rare variants (MAF 0.03–0.10), the HR should be >2.5 for a discovery. Since we selected an alpha-level lower than the commonly accepted genome-wide significance threshold 5 × 10−8, we calculated a BFDP with R library gap for nominal variant effects and FDR for the interaction models to estimate the validity of our findings69,70,71,72. In the BFDP analysis, the prior discovery probability was set to 0.0001 and HR alternatively to 1.3, 1.8, or 2.5, as in Escala-Garcia et al.23, and according to the estimated thresholds for discovery from the power analyses. SNPs with BFDP or FDR less than one-in-three were considered as interesting discoveries.

We performed additional survival analyses for the newly discovered SNPs only. These included multivariable survival model adjusted for tumor pathologic characteristics and analysis of breast cancer-specific death. Complete pathologic data were available for 1104 (36.7%) BRCA1 and 743 (37.0%) BRCA2 carriers, and the data on breast cancer-associated death for 2066 (68.7%) BRCA1 and 1591 (79.2%) BRCA2 carriers. However, data on both pathology and cause of death were available only for 683 (22.7%) BRCA1 and 544 (27.1%) BRCA2 carriers. The pathologic covariates in the multivariable model were coded as follows: tumor ER expression: categorical—no expression/positive expression (not included in the analyses of ER-stratified patient groups); tumor PgR expression: categorical—no expression/positive expression; tumor grade: linear—1, 2, 3; tumor size: linear—1 = less than 2 cm in diameter, 2 = diameter between 2 cm and 5 cm, 3 = larger than 5 cm in diameter; and lymph node status: categorical—affected/not affected. The validity of the proportional hazards assumption was evaluated for all discovered variants.

Candidate gene identification

The newly discovered survival SNPs were characterized in silico utilizing data from the 1000 genomes35 and Encode36 projects as integrated in databases LDlink73,74, RegulomeDB75, and GeneCards76,77. We retrieved all short genomic variants linked with the discovered survival variants in the 1000 genomes European data with D′ > 0.8. The proxy variants were subjected to RegulomeDB analysis and all variants with scores 1a–2f were considered as regulatory variants potentially contributing to the survival signal (Supplementary Data 1). The variant positions were matched to locations of protein-coding genes and validated enhancer regions to identify the plausibly target genes. Furthermore, we retrieved eQTL genes from GTEx78,79 releases V6, V6p, V7, and Westra et al. blood eQTL dataset80 for all functional proxy SNPs (D′ > 0.8) and for all strongly linked proxy SNPs (R2 > 0.8) irrespective of their functional annotation (Supplementary Data 2). The topological chromatin status at the rs57025206 locus was investigated in Rao et al.81 data from HMEC using the 3D genome browser82.

For a functional summary, we performed a systematic literature search using the gene symbol and any of the keywords “breast cancer”, “mammary”, “estrogen”, “immune”, “leukocyte”, and “lymphocyte”. If none of these queries returned relevant publications, the search was continued with gene symbol and “cancer” or with gene symbol alone. This was the first stage of literature search. In the second stage, the queries consisted of the gene symbol and any of the recurrent functional terms or interacting proteins detected in the first stage (Fig. 3).

The survival associations of the candidate genes’ mRNA expression were tested in the KM plotter database for breast cancer (http://kmplot.com/analysis/)83, restricting the analysis in the relevant group of mammary tumors: ER-negative tumors for target genes from the analyses of BRCA1 carriers, ER-positive tumors for genes from BRCA2 analysis, and all tumors for genes from the meta-analysis. The best cutoff for the categorical survival analysis was selected automatically and results with FDR ≤ 5% were reported. Furthermore, the linear association between candidate genes’ mRNA expression and patient survival was tested in the METABRIC data (EGAD00010000434, 1302 breast cancer patients) using Cox regression. Expression data were available for 29 of the 39 candidate genes. For each of the candidate genes, the samples were split into three categories based on 33 and 67% percentiles of the expression values, and the categories analyzed for 10-year overall survival. The results are included in the functional summary in Supplementary Table 6.

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.