Introduction

Breast cancer is the most common cancer in women worldwide and accounted for 2.1 million new cases and 627,000 deaths in 20181. Studies have shown a significant contribution of genetic factors to breast cancer risk2,3, yet the landscape of this contribution has not been fully elucidated. Mutations in high- and moderate-penetrance genes confer relatively high risks of breast cancer but are rare in the population and account for <5–10% of cases4. Genome-wide association studies (GWAS) have been successful in identifying common low-penetrance genetic variation and approximately 200 risk loci have now been identified5,6,7. The risk loci so far identified have provided clues to elucidating breast cancer tumorigenesis through previously unknown mechanisms. Additionally, when combined into risk scores, these polymorphisms can be used for breast cancer risk prediction8.

Despite the usefulness of GWAS, the majority of the GWAS studies have been performed among European ancestry populations9,10,11,12,13, it is unclear whether the same genetic risk factors are also important in other populations, which may limit the applicability of the findings to other groups14. The earliest GWAS conducted in African ancestry populations identified genetic variants at 5p15.33 (TERT-CLPTM1L) associated with estrogen receptor (ER) negative breast cancer15. A larger analysis of African ancestry individuals which included several consortia identified a SNP at 3q26.21 also associated with ER-negative breast cancer16. Some common susceptibility loci are shared across populations, and the shared disease-associated variants are more likely to be causal6,9,14.

Here we present, using a cross-ancestry GWAS approach in 248,000 women, genetic risk variants at 1p13.3, 5q31.1, 15q24, and 15q26.3 for overall breast cancer, and at 1q41 and 7q11.23 for ER-negative disease. The consistency of the directions of the risk for these loci in African and European samples increases the likelihood of their being causal variants.

Results

We discovered six loci containing seven SNPs significantly associated with breast cancer at P < 5 × 10−8 on cross-ancestry meta-analysis, with odds ratios (OR) ranging from 0.95 to 1.05 (Tables 1, 2; Supplementary Figs. 1, 2). Five SNPs were associated with overall breast cancer risk (rs17024628 at 1p13.3, rs2522057 at 5q31.1, rs1869959 at 15q24.1, rs60381548 at 15q24.2, rs181337095 at 15q26.3) and two were associated with ER-negative breast cancer (rs67931591 at 1q14 and rs1637365 at 7q11.2). The two SNPs at the 15q24 region were about 582 kb apart and independently associated with breast cancer risk. Four SNPs were within genes (rs67931591 in KCNK2, rs2522057 in C5orf56, rs1869959 in SCAMP2, and rs60381548 in SIN3A) and the others were in intergenic regions. The direction of the associations was consistent for the pooled African and European estimates. The estimates for overall and ER-negative breast cancer were generally consistent across the five contributing studies of African ancestry participants (Supplementary Table 2) and the BCAC European datasets (Supplementary Table 3).

Table 1 Novel breast cancer risk loci identified by cross-ancestry meta-analysis of African and European populations.
Full size table
Table 2 Association analysis of novel SNPs in cross-ancestry combined meta-analysis by estrogen receptor status.
Full size table

Conditional analysis revealed three additional independent signals significant at p < 10−4 at the 1p13.3 locus (rs116363925, rs114351980, and 1:109969874:C:T), two independent signals at 15q24 (rs113939578, rs12917507), and one each at 5q31.1 (5:132149322:G:GGCCGCCGCC) and 15q26.3 (rs117793215) for overall breast cancer risk. Another independent SNP at 1q41 that was associated with ER-negative breast cancer was rs5780828 (Table 3).

Table 3 Conditional regression analysis of top SNPs and others in the loci.
Full size table

Concerning pleiotropy, none of the SNPs identified above have been reported in previous GWAS associations at genome-wide significance with cancers. Associations with mosquito bite size and asthma had been reported for rs2522057 and SNPs in LD with this lead SNP. For the 15q24 region, associations with cardiovascular phenotypes have been previously reported for rs1869959 while body height, glomerular filtration rate, and type 2 diabetes have been associated with rs60381548 and SNPs highly correlated with this lead SNP (Supplementary Table 4).

The eQTL analysis of breast tumors revealed significant associations in four loci: 1p13.3, 5q31, 15q24.1, and 15q24.2 (Supplementary Table 5A). There were significant associations (P < 10−6) between the protective allele of rs17024629 (T allele) at 1q13.3 and increased expression of GSTM1, GSTM2, and GSTM4, which are located 19 kb, 31 kb, and 51 kb downstream of the SNP. At 5q31, the top SNP rs2522057, located 15 kb downstream of IRF1, was most significantly associated with the gene’s expression levels. At 15q24.1, rs1869959, located 35 kb upstream of MPI and 12 kb upstream of ULK3, was significantly associated with the expression of these two genes. The other top SNP at the 15q24 locus, rs60381548, located intron of SIN3A, 30 kb downstream of PTPN9, 162 kb downstream of SNUPN, and 212 kb upstream of SNX33 was correlated with all four genes. The 1q41 locus revealed a significant association between rs67931591 and PTPN14. The SNP at 7q11.23 was significantly correlated with STAG3L2, a pseudogene. Previous published report on normal breast tissues from the GTEx revealed associations between rs2522057 and the SLC22A5 gene, and between rs1869959 and the ULK3 gene (Supplementary Table 5B).

Functional annotation analyses pointed out relationships with genomic functional biofeatures for rs2522057, rs17024629, rs1869959, and rs60381548 or SNPs in strong LD with these top SNPs in breast tissue-originated cell lines (Supplementary Tables 6, 7A, 7B). Active enhancer and promoter states were found for SNPs in strong LD with rs2522057 (rs2188962, rs4705950, rs4705950, rs72797306, rs11741255) using the 25-state chromatin model. Additional associations were found with histone modifications. These included: H3K4me1 and H3K27ac enhancer peaks for: rs2522057 and other SNPs in strong LD (rs2188962, rs17622378, rs12521868, rs146604341, rs11951091, rs6866614, rs4705950, rs72797303, rs2706396, rs2522052, rs2706403, rs2706336, rs72797306, rs2248116, rs11741255); those in strong LD with rs17024629 (rs538388, rs560674, rs568686, rs669426, rs3850616, rs17024628); a SNP in strong LD with rs1869959 (rs7180432); and for the top SNP rs60381548. H3K4me3 and H3K9ac promoter peaks were found for: rs2522057 and other SNPs in strong LD (rs12515180, rs11951091, rs72797306); SNPs in LD with rs17024629 (rs538388, rs669426, rs3850616); rs1869959 and other SNPs in strong LD (rs4886613, rs936230).

We evaluated the consistency of the association of the identified loci in Latinos, and found the effect and direction of the association were consistent in 8 out of 11 evaluated variants (Supplementary Table 8). However, none of these consistent variants was statistically significant at p < 0.05 in the Latino study of 2385 cases and 6416 controls.

Discussion

We found seven variants associated with breast cancer risk among women of African ancestry that may contribute to better prediction of breast cancer risk and provide further insights into mechanisms of breast cancer carcinogenesis. Although the discovery of the loci is largely driven by effects in European ancestry populations, observation of risk loci in multiple ancestral populations lends credence to the chances of those variants being causal. We designed our current approach of cross-ancestry meta-analysis to uncover genetic variants shared across ancestry.

The SNPs identified in this study lie in regions that are close to genes that have been previously implicated in cancer. Interestingly we found three variants located within the introns of genes. One of the variants, rs67931591 was found in KCNK2 (also known as TREK1), which encodes the protein potassium channel subfamily K member 2, a member of the two-pore-domain background potassium channel family. Potassium channels are known to play a role in cancer and studies using TCGA data have shown associations with DNA methylation in the KCNK genes and triple negative breast cancer. Additionally, overexpression of KCNK5, KCNK9, and KCNK12 and under-expression of KCNK6 and KCNK15 were associated with triple negative breast cancer17. Other studies investigated expression of KCNK2 gene as potential prognostic markers. For example, Innamaa et al.18 found increased KCNK2 expression in human ovaries and a role in cell proliferation and apoptosis for KCNK2 modulators in ovarian cancer cell lines. Li et al.19 found differential expression of KCNK2, KCNK15 and KCNK17 in liver cancer cells compared to healthy tissue. KCNK2 has also been reported in amplified regions in a genome-wide scan of chromosomal alterations in esophageal squamous cell carcinoma20.

We found two independent SNPs at the 15q24 locus at about 582 kb apart (rs1869959 at 15q24.1 in the SCAMP2 intron and rs60381548 at 15q24.2 in the SIN3A gene). The SIN3A gene was associated with rs60381548 in the eQTL analysis of breast tumor in the present study. Switch-independent 3 family A (SIN3A) is a transcriptional regulator, that along with its paralog and corepressor play important roles in normal breast development, cancer and metastasis21,22,23. Furthermore, SIN3A mediates STAT3 transcriptional repressor activity24 and along with genes involved in histone modification such as HDAC and Lysine specific demethylase (LSD), inhibits several cancer genes including CASP7, TGFB2, CDKN1A, HIF1A, TERT and MDM225. Studies have shown key roles for SIN3A in breast cancer including sensitivity to chemotherapy25 and breast cancer progression26,27.

The other SNP at 15q24, rs1869959, is located in the intron of the SCAMP2 gene that codes for secretory carrier associated membrane protein 2 that functions as carriers to the cell surface in post-golgi recycling pathways28. The recent GTEx project pilot study found significant associations between the SNP and SCAMP2 in esophageal mucosa, ULK3 in breast mammary tissue, adipose, whole blood, and lung tissue29,30. We also found that rs1869959 was associated the expression of ULK3 in breast tumor. ULK3 is a serine threonine kinase that activates GLI2, a key component of the Hedgehog signaling pathway, and implicated in many cancers31,32.

Similarly, the C5orf56 gene harboring the rs2522057 SNP returned no interesting associations with cancer. However, nearby genes in the 5q31 locus included RAD50, that codes for a DNA repair protein, a part of the MRE11-RAD50-NBS1 complex33. Other nearby genes include SLC22A5 solute carrier family 22 member 5 encoding the OCTN2 (organic cation transporter protein), and IRF1 that encodes interferon regulatory factor 1. SLC22A5 is an estrogen-dependent gene whose expression is associated with ER status in breast cancer cell lines and tissue specimens34. Significantly decreased levels of SLC22A5 have been reported in colorectal cancer tissues compared to normal tissues in eQTL studies35. Moreover, eQTL studies report associations between rs2522057 and gene expression in several tissues including breast mammary tissue, lymphocytes, esophageal mucosa, lung, skeletal muscle, skin, thyroid and whole blood29,30. We found an association between rs2522057 and IRF1 expression levels in the eQTL analysis of breast cancer in this study lending support to the likelihood of involvement of the IRF1 gene in the mechanism of the SNP on breast cancer carcinogenesis. Additionally, IRF1 has been shown to have tumor suppressor functions in breast cancer through its inhibition of NF-kB36 and CASP8 activation and induction of apoptosis37.

The majority of GWAS-identified SNPs were located in non-coding regions of the genome, and three loci in the present study were found in intergenic regions. The closest gene to rs17024629 is AMPD2 (high adenosine monophosphate deaminase 2) and has recently been shown to predict worse outcomes in undifferentiated pleomorphic sarcoma38. Earlier studies39 found high expression levels of AMPD2 in hepatocellular carcinoma, though the levels did not differ substantially from those in the non-tumorous organ. It is noteworthy that our eQTL analysis did not find a significant association with AMPD2 expression. The carcinogen metabolism genes, GSTM1, GSTM2, and GSTM4 are also located in this region and our eQTL analysis of breast tumor revealed highly significant associations between rs17024629 and these genes. The GSTM1 null genotype has been associated with several cancers including cancers of the colorectum, oral cavity, lung, cervix, and stomach40,41,42,43,44,45,46,47. In eQTL studies, GSTM4 was significantly associated with gene expression in several tissues including the aorta, lungs, tibia nerve and whole blood29,30.

The rs1637365 SNP at the 7q11.23 locus is near the CASTOR2 gene (cytosolic arginine sensor for MTORC1 protein, also known as GATSL1, GATS-like protein 1). The CASTOR proteins are arginine sensors that function as negative regulators of the TORC1 signaling pathway, an often dysregulated pathway in human cancer, through the GATOR complex, inhibiting mTORC148,49. The rs181337095 SNP is located 6 kb 5′ of RP11-168G16.2, an antisense DNA.

A potential limitation of this study is the different genotyping platforms used by the different consortia. However, stringent QC measures pre- and post-imputation were carried out. Additionally, the meta-analysis did not reveal significant heterogeneity across studies. Secondly, the sample size for ER-negative breast cancer cases was relatively small, thus reducing the precision of the estimates and providing less power for detecting risk loci. The third limitation is related to the additional SNPs identified at the same loci with the index SNPs from the conditional regression analysis. The regression procedures were based on a liberal p value cutoff of 10−4, and the chance that some of the identified SNPs could be spurious findings cannot be ruled out. Another noteworthy point is that identification of genetic variants in GWAS is just the first step of the discovery of true causal variants and genes associated with breast cancer. Further studies are needed, including in vitro and in vivo functional studies to elucidate the mechanisms by which identified putative causal variants are acting and identify the targeted genes, Finally, although the direction and strength of the associations were consistent between African and European populations, and mostly consistent with Latino populations, we could not find statistically significant replication of the identified variants, which are likely due to the modest sample sizes of the Latino study.

Our study found six loci that could provide further insights into pathways for breast cancer carcinogenesis. The genetic variants that shared across ancestry populations makes them possible causal variants. Functional studies on these loci are desirable to identify causal variants and elucidate the mechanisms of breast cancer carcinogenesis. In addition, future studies can evaluate these variants for breast cancer risk prediction, particularly in African ancestry populations.

Methods

Study population

Data for this study were obtained from four consortia of African ancestry populations (ROOT, AMBER, AABC, and BCAC-African ancestry)16 and the Ghana Breast Health Study (GBHS)50,51, with a combined sample size of 19434 participants including 9241 cases and 10193 controls (Supplementary Table 1). Estimates from these studies were meta-analysed to generate pooled African ancestry estimates of breast cancer risk. Additionally, we used summary estimates (odds ratios, ORs) of breast cancer from European ancestry BCAC datasets (GWAS, iCOGs and OncoArray) with a combined sample size of 228,951 (122,977 cases and 105,974 controls)6.

Genotyping and quality control

Genotyping and quality control (QC) procedures have been described in detail for the three consortia16 and the BCAC European ancestry data6. The AABC was genotyped using the Illumina Human 1M-Duo BeadChip. After QC, a total of 3007 cases (1518 ER-positive, 987 ER-negative) and 2720 controls remained in the analysis52. Genotyping in the ROOT consortium was done using Illumina HumanOmni 2.5-8v1 array and 1657 cases (374 ER-positive, 403 ER-negative) and 2029 controls passed QC. In the BCAC-African ancestry consortium, genotyping was done using the Illumina OncoArray (260K GWAS backbone) and after removing overlapped samples between OncoArray with AABC, AMBER and ROOT and samples failed in QC, a total of 2271 cases (1130 ER-positive, 613 ER-negative) and 1406 controls remained for analysis. The Illumina MEGA array was used for genotyping in the AMBER consortium, and 1407 cases (952 ER-positive, 385 ER-negative) and 2408 controls remained in analysis passed QC. In the GBHS, Illumina Global Screening Array was used for genotyping, and 899 cases (296 ER-positive, 277 ER-negative) and 1630 controls were included in analysis. Imputation for all studies was done using the cosmopolitan reference panel in the 1000 Genomes Project (Phase 3 release).

In addition, we examined the association between the identified SNPs of interest and breast cancer risk in a GWAS of Latinos (2385 cases and 6416 controls). Details of the genotyping, QC and data analysis have been published53.

Data analysis

GWAS

In the ROOT and AABC GWAS studies, genotyped SNPs were analyzed and imputed with imputation score >0.3 and minor allele frequency >0.01 to account for uncertainty in imputation. Unconditional logistic regression was used to examine the association of each SNP and breast cancer risk adjusting for age, study site and eigenvectors from Principal Components Analysis (PCA). In the ROOT GWAS, the first four eigenvectors were used to control for population stratification as only the first 4 eigenvectors were associated with case status. The AABC GWAS adjusted for the first 10 eigenvectors from the PCA. OR and 95% confidence intervals (CI) were calculated from the multivariable logistic regressions. All tests of statistical significance were two sided. Using similar methods, separate analyses were conducted to compare ER-positive and ER-negative breast cancers with controls. The AMBER consortium estimated ORs and P values using unconditional logistic regression, adjusting for 10-year age group, sample type (saliva, blood, other), study (Black Women’s Health Study (BWHS) versus others) and PCs that associated with breast cancer at P < 0.1. The GBHS estimated per-allele ORs and 95% CI for each SNP on allele counts (dosages) using unconditional logistic regression adjusting for the first ten principal components, self-reported ethnicity and age. In the Oncoarray African ancestry samples, a total of 27 million SNPs with MAF ≥ 0.1% and imputation quality score ≥0.3 were included in the analysis. PCs were estimated using EIGENSTRAT. ORs and P value of each SNP were estimated using unconditional logistic regression, adjusting for age, study (Women of African Ancestry Breast Cancer Study—WAABCS versus other) and the first ten PCs.

The BCAC European study used a two-stage imputation approach, using SHAPEIT2 for phasing and IMPUTE version 2 for imputation. The first ten principal components and country were adjusted for in the logistic regression, and per-allele ORs and standard errors were computed6.

Meta-analysis

Regression coefficient estimates from the five contributing African ancestry studies were combined in a fixed effects meta-analysis using METAL54. Variants associated with breast cancer at P < 0.05 from the African ancestry meta-analysis were then combined in another fixed effects meta-analysis with the coefficients from the BCAC European ancestry data. Heterogeneity in both meta-analyses was assessed using the I2 statistic. SNPs that were significant genome-wide (P < 5 × 10−8) in the cross-ancestry meta-analysis, and >500 kb away from the 180 loci known to be associated with breast cancer risk were identified5,6. Conditional analysis below confirmed the identified loci. All analyses were done separately for ER-positive, ER-negative, and overall breast cancer risk.

Regression analysis conditional on index SNPs

In order to identify independent SNPs in the identified loci, conditional analysis was done in each of the regions, including all variants in the flanking ±500 kb region of the lead SNP. The 15q24 region had two SNPs about 582 kb apart that were both genome-wide significant (see results for details). Hence, all variants in the region extending from 500 kb upstream of the proximal SNP and 500 kb downstream of the other SNP were included in the conditional analysis for this region. We used the GCTA software with the –COJO option55, that utilizes summary statistics and population-specific linkage disequilibrium (LD) from 1000 Genomes Project, for the computation of conditional beta coefficients. SNPs significant at p < 10−4 after adjusting for lead SNP were considered as independent signals. The p < 10−4 cutoff was derived by applying a factor of 3000 (the ratio of the 3 billion base pairs genome-wide to the 1 million base pairs in each region in the conditional analysis) to the GWAS significance of 5 × 10−8. This procedure was repeated until no additional independent signals were significant. In addition to the conditional analysis involving the lead SNP and one other candidate SNP, we also determined joint ORs including all independent loci in the same model. Separate analyses were done for African and European ancestry data, and the estimates from the conditional analysis were combined in a meta-analysis.

Functional annotation

The functional annotations of the SNPs were determined using HaploReg v4.156. Using data from ENCODE57 and the Roadmap Epigenomics Consortium58, we examined the chromatin states including core 15-state model and 25-state model using 12 imputed marks, H3K4me1 and H3K27ac (enhancers), and H3K4me3 and H3K9ac (promoters) for each identified SNP and other SNPs in strong LD with these lead SNPs (>0.8). We also assessed evolutionary conserved regions, DNase hypersensitivity sites, and variant effect on regulatory motifs, proteins bound and eQTL hits from previous studies.

eQTL analysis

We carried out a cis-eQTL analysis to understand possible target genes in the six loci. All genes within ±1 MB around each index SNP were evaluated and gene expression in breast tumors from TCGA breast cancer patients (African ancestry, n = 164 and European ancestry, n = 778) were used in the analysis. A linear regression model estimated additive effects for each SNP, adjusting for age, ancestry, copy number variation, batch effect, and molecular subtype. Separate analyses were done for African and European ancestry samples and the estimates were meta-analysed to obtain overall estimates. Bonferroni significance levels were applied to determine statistical significance. We also checked associations between the identified loci and gene expression in several tissues, including normal breast, that had been published from previous eQTL analyses on the Haploreg website.

Allelic pleiotropy

We assessed the GWAS catalog (www.ebi.ac.uk) for previously reported associations for the identified lead SNPs and all other SNPs in LD with r2 > 0.4 and phenotypes.

Ethical approval

Informed consent was obtained from all subjects included in the analysis. The relevant ethical review boards at all participating institutions approved study protocols.