Identification of 31 loci for mammographic density phenotypes and their associations with breast cancer risk

Mammographic density (MD) phenotypes are strongly associated with breast cancer risk and highly heritable. In this GWAS meta-analysis of 24,192 women, we identify 31 MD loci at P < 5 × 10−8, tripling the number known to 46. Seventeen identified MD loci also are associated with breast cancer risk in an independent meta-analysis (P < 0.05). Mendelian randomization analyses show that genetic estimates of dense area (DA), nondense area (NDA), and percent density (PD) are all significantly associated with breast cancer risk (P < 0.05). Pathway analyses reveal distinct biological processes involving DA, NDA and PD loci. These findings provide additional insights into the genetic basis of MD phenotypes and their associations with breast cancer risk.

P ercent density (PD), the percentage of the breast area that appears radiodense or light on a mammogram, is one of the strongest risk factors for breast cancer, but the biological basis for this association is poorly understood 1,2 . Women with ≥75% density on a mammogram have a 4 to 5-fold increased risk of breast cancer compared to women with little or no dense tissue, independent of other known risk factors 2,3 . PD is a composite of two phenotypes: the dense area (DA) reflecting the amount of fibroglandular tissue in the breast, and the nondense area (NDA) consisting of predominantly fatty tissues that appear radiotranslucent or dark on a mammogram 4 . Recent studies have shown that NDA is associated with decreased breast cancer risk independently of DA, suggesting that breast adipose tissues play an important role in normal mammary gland growth and function 5,6 . PD, DA, and NDA each have heritability estimates of over 50% in twin studies [7][8][9] . However, only 15 independent genome-wide significant loci with P < 5 × 10 −8 have been identified to date, together explaining less than 1-3% of the total variance of mammographic density (MD) phenotypes [10][11][12][13][14] .
In this genome-wide association study (GWAS) meta-analysis of 24,192 women screened with full-field digital mammography (FFDM), we identify 31 MD loci, of which 17 also are associated with breast cancer in an independent study of over 200,000 breast cancer cases and controls 15 . These findings triple the total number of independent genome-wide significant MD loci now mapped to 46, enabling the first genetic pathway analyses and Mendelian randomization analyses to evaluate the causal nature of the association of MD phenotypes with breast cancer risk. Pathway analyses reveal distinct biological processes involving DA, NDA, or PD loci. Mendelian randomization analyses show that genetic estimates of DA, NDA, and PD are all significantly associated with breast cancer risk. These findings provide additional insights into the genetic basis of MD phenotypes and their relationship with breast cancer risk.

Results
GWAS of MD phenotypes. This GWAS meta-analysis comprised a total of 24,192 non-Hispanic white women with MD phenotypes measured centrally using Cumulus software 16 . The first study included 20,311 women screened using Hologic FFDM machines, and the second study included an independent sample of 3881 women screened using General Electric (GE) FFDM machines. Women in the GE cohort were 2.7 years younger, had lower BMI, and were less likely to be postmenopausal compared with women in the Hologic cohort (Supplementary Table 1). On average, DA was 1.1 cm 2 higher, NDA 31.2 cm 2 lower, and PD 4.0 percentage points higher in the GE cohort compared with the Hologic cohort ( Supplementary Fig. 1). PD, computed by DA divided by the total breast area (DA + NDA), was strongly correlated with DA (R, 0.8) and NDA (R, −0.8), and DA was moderately negatively correlated with NDA (R, −0.35) in both cohorts, as expected.
In the Hologic study, 37 SNPs were associated with MD phenotypes at P < 5 × 10 −8 . In the GE study, 3 of these SNPs could not be confirmed and were excluded, while 18 additional SNPs with P < 5 × 10 −5 in the Hologic study reached genomewide significance in the combined meta-analysis. In total, 52 SNPs at 40 independent chromosomal regions (loci) were associated with DA, NDA, and/or PD in the same directions in both studies, and met the conventional genome-wide significance threshold of P < 5 × 10 −8 (Fig. 1). The genomic inflation factors for the GWAS meta-analyses of DA, NDA, and PD were 1.06, 1.08, and 1.07, respectively, indicating that there was little evidence of uncontrolled population substructure ( Supplementary  Fig. 2).
We identified 39 previously unreported MD SNPs at 31 independent loci, including 16 SNPs for DA, 13 for NDA, and 12 for PD with P < 5 × 10 −8 in the GWAS meta-analysis, adjusting for ln(BMI), age at mammography, and principal components of European ancestry (Table 1, Supplementary Table 5). A single SNP showed the strongest association with both DA and PD at two loci (2q35, 6p22.3). Of the 16 DA loci, one (3p25.2) was also significantly associated with both NDA and PD, five (2p24.1, 2q35, 4q28.1, 6p22. 3, 19q13.33) were also associated with PD but not NDA, and 10 were associated with DA only. Of the 13 NDA loci, one (3p25.2) was also significantly associated with both DA and PD, three (2p23.3, 10q21.2(1), 11q24.3) were also associated with PD but not DA, and nine were associated with NDA only. Sensitivity analyses of the NDA SNPs showed that all 13 remained significant after additional adjustment for BMI using three polynomial terms (BMI, BMI 2 , and BMI 3 ) in addition to ln (BMI); and only 2 of the 13 NDA SNPs (rs6718628 and rs4132228) were associated with BMI at P < 0.05 in models unadjusted for BMI (Supplementary Table 13). Of the 12 PD loci, nine were also significantly associated with DA and/or NDA, and three (8p12, 13q13.3, 15q26.1) were associated with PD only. These three SNPs had opposite directions of association with DA and NDA that were not genome-wide significant when analyzed separately, but reached statistical significance for the PD composite measure of DA and NDA. Analyses stratified by menopausal status showed that MD SNP effects were similar in premenopausal and postmenopausal women (Supplementary Table 15).
Associations with breast cancer risk. We evaluated whether the newly identified MD loci were also associated with breast cancer risk using data from an independent sample of 122,977 cases and 105,974 controls of European ancestry from the Breast Cancer Association Consortium (BCAC) and Discovery, Biology and Risk of Inherited Variants in Breast Cancer Consortium (DRIVE) 15 . We found that 24 MD SNPs at 17 loci were associated with breast cancer risk at P < 0.05 (Table 2, Supplementary Table 6). Of these SNPs, 15 MD SNPs at 10 loci were in linkage disequilibrium with previously reported breast cancer susceptibility alleles 15         19q13.33, 20q13.13) were associated with both MD and breast cancer risk.
To explore the extent to which DA, NDA, and PD may be associated with breast cancer through shared underlying genetic factors and biologic pathways 4 , we performed Mendelian randomization analyses using the weighted median method 18 and summary statistics for all prior and new MD loci from this study and BCAC/DRIVE 15 . For each standard deviation (SD) increment in DA, NDA, and PD, the estimated odds ratios (95% confidence interval) for breast cancer were: 1.45 (1.30-1.61; P < 0.001), 0.84 (0.73-0.98; P = 0.029), and 1.68 (1.44-1.96; P < 0.001), respectively. Sensitivity analyses using mode-based estimates that are also robust to violations of the instrumental variable assumptions, but less powerful than the weighted median method showed similar results (Supplementary Table 14) 19 . These estimates were remarkably similar to a meta-analysis of 13 observational studies that reported adjusted odds ratios of breast cancer associated with each SD increment of DA, NDA, and PD of: 1.  6 . These findings support a biological basis for the positive association of DA and PD with breast cancer risk, and inverse association of NDA with breast cancer risk, in observational studies.
To complement the Mendelian randomization analyses, we estimated the genetic correlation between each MD phenotype and breast cancer based on all SNPs genome-wide. LD Score regression 20,21 of 779,828 SNPs using summary statistics from this study and the BCAC/DRIVE breast cancer GWAS 15 yielded estimates of 0.27 (P = 5.5 × 10 −6 ), −0.14 (P = 0.014), and 0.27 (P = 7.7 × 10 −10 ) for the genetic correlations of DA, NDA, and PD, respectively, with breast cancer. The significant positive genetic correlations of DA and PD with breast cancer, and significant inverse association of NDA with breast cancer were consistent with the Mendelian randomization results, as well as evidence from observational studies 6 , supporting the shared genetic bases of all three MD phenotypes and breast cancer.
Functional analyses. To identify potentially functional variants at the 31 MD loci, we examined whether any of the lead SNPs were associated with: protein-coding variants (Supplementary Table 2); gene expression levels in mammary tissue, primary fibroblast cells, subcutaneous fat, visceral fat, or whole blood (Supplementary Table 3); or promoter and enhancer regions in mammary epithelial cells or mammary fibroblasts (Supplementary Table 4). Regional association plots showed the nearby genes and linkage disequilibrium patterns in Europeans for each of the newly identified loci for DA, NDA, and PD (Supplementary Figs. [3][4][5]. MD SNPs at three loci (2p23.3, 15q25.2, 19q13.33) were strongly correlated (r 2 ≥ 0.80) with nonsynonymous mutations, and all of these SNPs were also significant (FDR < 0.05) expression quantitative trait loci (eQTLs) ( Supplementary  Tables 2 and 3). SNPs on 2p23.3 were associated with NDA and PD. The lead PD SNP rs11676272 encoded the S107P missense mutation in ADCY3 (adenylate cyclase 3), which catalyzes the formation of the secondary messenger cyclic adenosine monophosphate involved in signal transduction and metabolic processes. The lead SNPs for NDA rs6718628 and PD rs11676272 were also significant eQTLs associated with increased ADCY3 expression in subcutaneous and visceral fat, and whole blood. On 15q25.2, the lead NDA SNP rs1812707 was strongly correlated (r 2 = 0.88) with the V661L missense mutation variant in ADAMTSL3, involved in protein glycosylation and catabolism, and was also a significant eQTL in fibroblasts and subcutaneous fat for the GOLGA6L5P pseudogene, about which little is known.
On 19q13.33, the lead SNPs for DA rs492602 (r 2 = 0.99) and PD rs1704773 (r 2 = 0.89) were tightly linked with the W154X nonsense mutation in FUT2 (fucosyltransferase 2) resulting in a truncated protein. FUT2 is involved in the production of histoblood group antigens, and exhibits the non-secretor phenotype (lack of antigens in epithelial mucosa and exocrine secretions) when inactivating mutations are present. The lead DA and PD SNPs were also significantly associated with decreased expression of FUT2 in mammary tissue and fibroblasts, and increased expression in fibroblasts of the nearby MAMSTR gene, encoding a transcriptional regulator. This regulatory activity may be mediated in part by strong correlation (r 2 ≥ 0.80) with an intronic variant of MAMSTR within a promoter and enhancer-like region in mammary epithelial cells and mammary fibroblasts (Supplementary Table 4).
The lead MD SNPs at 13 loci were significantly associated (FDR < 0.05) with gene expression levels in normal human mammary tissue, primary fibroblast cells, subcutaneous fat, visceral fat, or whole blood, which were the tissues most closely related to cell types in the breast available in GTEx 22 MD SNPs at 13 loci, or strongly correlated (r 2 ≥ 0.80) variants nearby, were located within regions with promoter or enhancer activity in normal human mammary epithelial cells or mammary fibroblasts using data from ENCODE 24,25 (Supplementary  Table 4). On 10q21.2, the new NDA SNP in an intron of ZNF365 was perfectly linked (r 2 = 1.0) with an enhancer-region variant in mammary fibroblasts. ZNF365 is involved in regulating neuronal growth and DNA repair. On 14q24.1 the lead DA SNP was an intronic variant in RAD51B located within an enhancerlike region in mammary fibroblasts. The RAD51 protein family is essential for DNA repair by homologous recombination, and interacts with the major breast and ovarian cancer susceptibility genes BRCA1 and BRCA2. On 18q21.33 the lead DA SNP was an intronic variant in BCL2 located within an enhancer-like region in mammary fibroblasts. BCL2 suppresses apoptosis and constitutive expression is thought to cause follicular lymphoma.
Enrichment of MD loci in fibroblast regulatory regions. To identify cell types through which DA, NDA, and PD loci influence their respective phenotypes, we tested for the enrichment of all independent prior or new loci for each MD phenotype in the regulatory regions of 125 diverse human cell and tissue types using the Uncovering Enrichment through Simulation (UES) method 26 (Supplementary Table 9). Regulatory regions were defined by DNase I hypersensitive sites sequencing (DNase-seq) experiments available from the ENCODE and Roadmap Epigenomic consortia 27 . We found that NDA loci (n = 17) were significantly enriched in the open chromatin regions of fibroblast cell lines derived from three different normal human tissues at the Bonferroni threshold of P < 0.0004 accounting for 125 tests. NDA loci were observed within regulatory regions more often than expected under the null distribution in fibroblasts isolated from the lung (3.6-fold; P < 0.0001), skin (3.8-fold; P = 0.0001), and heart (3.2-fold; P = 0.0003). Suggestive enrichment of NDA loci at P < 0.01 was also found in the regulatory regions of 13 additional fibroblast cell lines from different tissues, including mammary fibroblasts (2.7-fold; P = 0.0066). The regulatory activity of NDA loci in mammary fibroblasts was further supported by our functional analyses showing that 6 of 13 new NDA loci were associated with promoter or enhancer regions in mammary fibroblasts (Supplementary Table 4). DA (n = 28) and PD (n = 20) loci were not significantly enriched in the regulatory regions of any normal human cell types. However, suggestive evidence of enrichment (P < 0.01) was found specifically in 6 fibroblast cell lines for DA loci, and 2 fibroblast cell lines for PD loci, and not in any other normal human cell type (Supplementary Table 9).
Biological pathways. We identified distinct biological pathways significantly enriched (FDR < 0.05) for DA, NDA, and PD loci using DAVID 28,29 to perform gene set enrichment analyses of all prior and new genome-wide significant loci. We found that NDA loci (n = 17) were significantly enriched by over ninefold for genes involved in mammary gland development (Supplementary  Table 11). Other NDA pathways were related to metabolism, cell differentiation, and reproduction. DA loci (n = 28) were significantly enriched for genes involved in reproduction, apoptosis, metabolism, and signaling (Supplementary Table 10). PD loci (n = 20) were significantly enriched for genes involved in anatomical structure development and metabolism (Supplementary Table 12). Additional pathways implicated by combining all MD loci (n = 46) were involved in regulation of gene expression, nucleic acid binding and metabolism, and cell proliferation (Supplementary Data 1).
Heritability. The proportion of phenotypic variance explained by all genotyped SNPs estimated using GCTA 30 10 . SNP-based heritability estimates represent the upper bound on the total proportion of phenotypic variance that could be explained by GWAS of common variants, and the lower bound for narrowsense heritability estimated from twin studies because they do not account for rare variants that are not in linkage disequilibrium with the genotyped SNPs 32 . Heritability estimates from twin studies can also be influenced by nonadditive genetic effects and shared environmental effects 32 . Altogether, the newly identified and previously known MD loci explained 12.3%, 9.1%, and 8.7% of the SNP-based heritability for DA, NDA, and PD, respectively, compared with 3.3%, 2.1%, and 3.2% explained by previously known MD loci.

Discussion
High MD is one of the strongest and most common risk factors for breast cancer, and has been estimated to account for up to one-third of all breast cancers 4 . PD, DA, and NDA are all highly heritable and significantly associated with breast cancer risk in observational studies 6 . However, the biological bases for how these breast tissue phenotypes are related to breast cancer development are poorly understood. In this GWAS meta-analysis, we identified 31 loci for MD phenotypes, tripling the total number of genome-wide significant loci from 15 previously to 46 presently. Seventeen of the MD loci also were associated with breast cancer risk in an independent large meta-analysis, identifying potential new breast cancer susceptibility alleles at seven loci. Mendelian randomization and genetic correlation analyses provided further evidence of the shared genetic etiology of all three MD phenotypes with breast cancer.
Among the seven loci newly associated with both MD and breast cancer risk, one locus on 3p25.2 was significantly (P < 5 × 10 −8 ) associated with all three MD phenotypes. Intronic variants in PPARG (peroxisome proliferator activated receptor gamma) were associated with DA, PD, and breast cancer risk in the same direction, and with NDA and breast cancer risk in the opposite direction, as expected. PPARG is a member of the nuclear receptor family of ligand-activated transcription factors, and a regulator of adipocyte differentiation 33 . PPARG has been shown to inhibit transcription of aromatase, the rate-limiting enzyme in estrogen biosynthesis, in primary breast adipocytes 34 . Higher aromatase expression has been observed in dense breast tissue 35 , and treatment with aromatase inhibitors lowers breast cancer risk although associations with MD changes have been less consistent 36 . Estrogen may influence breast cancer risk as well as breast tissue composition through its proliferative effects on mammary cells 37 . PPARG is therefore plausibly involved in the development of DA, NDA, and PD, as well as breast cancer.
Five additional loci were newly associated with both DA and breast cancer risk in the same direction, identifying LMNB1, MKX, PRKG1, MRPL17, and SMIM25 as candidate genes for both DA and breast cancer risk. On 5q23.2, rs6885843 was positively associated with DA and breast cancer risk, and located 44 kb upstream of LMNB1 (lamin B1) involved in autosomal dominant adult-onset leukodystrophy 33 . On 10p12.1, rs2642278 was positively associated with DA and breast cancer risk, and an eQTL for MKX (mohawk homeobox), which plays a role in cell adhesion 33 . On 10q21.1, rs1892368 is an intronic variant in PRKG1 (protein kinase cGMP-dependent 1) that was inversely associated with both DA and breast cancer risk (P = 4.8 × 10 −7 ). PRKG1 is a key mediator of the nitric oxide/cGMP signaling pathway important in many signal transduction processes 33 . On 11p15.4, rs11040963 was inversely associated with both DA and breast cancer risk (P = 4.5 × 10 −4 ), and located 8 kb upstream of MRPL17 (mitochondrial ribosomal protein L17) involved in protein synthesis in mitochondria 33 . On 20q13.13, rs17196752 was positively associated with DA and breast cancer risk, and uncorrelated (r 2 = 0.0001) with the closest known breast cancer risk allele rs6122906 15 . rs17196752 is an intronic eQTL within an enhancerlike region in mammary fibroblasts that down-regulated the mammary tissue expression of SMIM25 (Small Integral Membrane Protein 25). While little is known about SMIM25 function, rs17196752 has been associated with white blood cell traits that have been linked to cancer and other systemic diseases 38 .
Finally, on 19q13.33, SNPs for DA (rs492602) and PD (rs1704773) were tightly linked with a FUT2 protein-truncating mutation associated with the nonsecretor phenotype for histoblood group antigens. Both SNPs also were associated with increased expression of the MAMSTR transcriptional regulator in fibroblasts, and a MAMSTR regulatory region in mammary epithelial cells and mammary fibroblasts. rs492602 has been associated with serum lipid levels 39 but not breast cancer, and rs1704773 was associated with PD and breast cancer (P < 0.05) in opposite directions, indicating that FUT2 and MAMSTR are candidate genes for DA and PD but may not be directly associated with breast cancer risk.
Among the ten new MD loci associated with known breast cancer susceptibility alleles, all but one were associated with MD phenotypes and breast cancer risk in consistent directions. SNPs on 2p24.1 were associated with both DA and PD in directions opposite to their known associations with breast cancer risk 15 . The DA SNP rs11684853 was tightly linked (r 2 = 0.99) with variants in an enhancer element with strong activity in mammary epithelial cells and mammary fibroblasts. The PD SNP rs34331777 was about 100 kb away from the tumor suppressor gene OSR1 (Odd-skipped related 1) encoding a zinc-finger transcription factor that acts on the p53 and Wnt/β-catenin signaling pathways 40 , and the MIR4757 microRNA that could influence the translation of multiple target mRNAs with different effects on MD and breast cancer risk.
There are several potential biological mechanisms through which higher DA and PD may increase breast cancer risk. Dense areas of the breast contain more epithelial cells, fibroblasts, and collagen than nondense areas that contain more adipocytes 41 . Increased collagen alignment and stiffness of the extracellular matrix (ECM), associated with dense breast tissue, have been shown to induce malignant phenotypes in normal mammary epithelial cells 42 . Fibroblasts produce collagen and other ECM components, as well as proteases involved in ECM remodeling. In addition, fibroblast signaling is a critical determinant of normal mammary epithelial and adipocyte cell development and differentiation 43 , and cancer-associated fibroblasts can stimulate breast tumor progression 44 . Our finding that loci for all three MD phenotypes were enriched in the regulatory regions of fibroblasts more than in any other cell type is consistent with a key role of fibroblasts in regulating the stromal environment and normal breast tissue composition, as well as aberrant growth in breast tumors.
Mammary epithelial cells from high density tissues also have been shown to have greater DNA damage response signaling and shorter telomeres compared with mammary epithelial cells from low density tissues 45 . DNA damage may increase DA and PD, and decrease NDA, by repressing CD36 expression in mammary fibroblasts, which induces increased ECM deposition and decreased lipid storage in nonmalignant breast tissue 44 . CD36 is a widely expressed glycoprotein receptor that binds to a broad range of ligands, including ECM proteins and lipids, and modulates adipocyte differentiation, lipid metabolism, angiogenesis, apoptosis, cell-ECM interactions, and immune signaling 44,45 . Importantly, CD36 expression is primarily controlled by the PPARG transcription factor 45 , which we found to be associated with all three MD phenotypes and breast cancer risk. The etiologic role of DNA repair and apoptosis genes in dense breast tissue is further supported by our findings that variants in the RAD51B (14q24.1) DNA repair gene, and FAF1 (1p32.3) and BCL2 (18q21.33) apoptosis genes are significantly associated with DA, and that DA loci are enriched for genes in the apoptosis pathway.
NDA has been inversely associated with breast cancer risk independently of DA in observational studies 5,6 , but the etiologic nature of this association and underlying mechanisms are uncertain. Our findings demonstrating a significant inverse association of genetically estimated NDA and breast cancer risk in both Mendelian randomization and genetic correlation analyses provide strong evidence that this association is caused by shared underlying genetic and biological pathways. A limitation of this study is that relatively few women were diagnosed with breast cancer following the mammogram, precluding mediation analyses of the extent to which SNP associations with breast cancer risk are explained by MD phenotypes. A limitation of Mendelian randomization studies generally is the potential for bias due to horizontal pleiotropy, although the weighted median method is a relatively robust and statistically powerful approach 18,19 . While NDA and BMI effects are difficult to disentangle, the adipose tissues within the breast may play a more direct role in breast cancer etiology than distant adipose tissues. Mammary adipocytes secrete adipokines that modulate the stromal environment, and constitute a local source of lipids and metabolites that influence mammary epithelial cell growth and function 46 . In vivo models have shown that the mammary adipose environment is critical for mammary gland growth and development 46 .
The finding that NDA loci, but not DA or PD loci, are significantly enriched for genes involved in mammary gland development supports a key role of the nondense fatty tissues in breast health. Three of the four new NDA loci enriched for genes involved in mammary gland development were associated with NDA and breast cancer risk 15 in opposite directions, implicating GLI2, NCOA1, and NTN4 as candidate genes for both NDA and breast cancer. GLI2 at 2q14.2 encodes a zinc finger transcription factor that mediates hedgehog signaling 33 . NCOA1 expression was upregulated in mammary tissue by the NDA SNP at 2p23.3, and encodes a transcriptional coactivator for steroid and nuclear hormone receptors that stimulates transcriptional activity in a hormone-dependent fashion 33 . NTN4 expression was upregulated in mammary tissue by the NDA SNP at 12q22, and encodes a member of the netrin protein family involved in neuronal growth, angiogenesis, and tumorigenesis 33 . Netrin 4 has been implicated in controlling epithelial cell branching morphogenesis in the breast 47 . These findings provide insights into the genetic basis for the inverse association of NDA with breast cancer risk.
In summary, this GWAS of 24,192 women from two independent population-based cohorts screened using Hologic or GE digital mammography, and MD phenotypes measured centrally using Cumulus, identified 31 MD loci, and new candidate genes for MD and breast cancer risk. The study findings support the etiologic role of NDA as well as DA and PD as modifiable risk factors for breast cancer that provide potential for intervention. With this study, all 46 genome-wide significant loci identified to date explain 12.3%, 9.1%, and 8.7% of the SNP-based heritability for DA, NDA, and PD. Future studies are needed to discover additional MD loci and to elucidate the different roles that the fatty and dense breast tissue components play in breast health and cancer risk.

Methods
Study design. We conducted a GWAS meta-analysis within the Research Program on Genes, Environment and Health (RPGEH) administered by Kaiser Permanente Northern California (KPNC) Division of Research 48,49 . RPGEH is populationbased and participants were not selected based on any disease phenotype. All participants completed a health survey, and over 100,000 individuals provided a DNA sample that was genotyped genome-wide; this sample constitutes the Genetic Epidemiology Research on Adult Health and Aging (GERA) study 49 . Written informed consent was obtained from all participants. Institutional Review Board approvals for this study were obtained from KPNC, Stanford University, and the Icahn School of Medicine at Mount Sinai.
Hologic study: The first GWAS included 20,311 non-Hispanic white women who underwent bilateral screening mammography at age 39-80 years during 2004-2013 at one of 36 KPNC clinics using Hologic FFDM machines. Processed (for presentation) images were retrieved from the KPNC imaging archive. Hologic images were downsampled from a pixel size of 70 microns to 200 microns. The resulting image resolution exceeded that of the computer monitors used to view the images, and was therefore unlikely to influence the density measurements. A median filter with a 3-pixel radius was applied to the downsampled Hologic images to reduce digital noise and make the images appear more like screen film mammograms, which we found to improve the reproducibility of density measurements 50 . Hologic images were randomly assembled into 23 batches of up to 1100 images each, including randomly selected replicates for quality control.
GE study: The second GWAS included an independent sample of 3881 non-Hispanic white women who underwent bilateral screening mammography at age 38-77 years during 2004-2013 at one of 11 KPNC clinics using GE FFDM machines. GE images processed using Tissue Equalization software were retrieved from the KPNC imaging archive, and downsampled from a pixel size of 94 microns to 200 microns. We found that denoising of the GE images did not improve the reproducibility of the density measurements and therefore did not denoise the downsampled images 51 . GE images were randomly assembled into 6 batches of up to 700 images each, including randomly selected replicates for quality control.
Density assessments. Women were excluded if their mammograms contained breast implants (3.6%), did not contain the entire breast (1%), or were unreadable/ unavailable (2.6%) 51 . Women with a history of bilateral breast cancer (0.06%) were also excluded. For women with a history of unilateral breast cancer, we measured MD phenotypes using the contralateral (unaffected) breast image from the closest prediagnostic exam following the RPGEH survey when available, or prior to the survey otherwise 50 . For women with no history of breast cancer, we used the left breast image from the closest post-survey exam, except for a random 10% subset of women for whom the right breast image was used to blind the reader to cancer history. All density measurements were performed using a single craniocaudal view.
MD phenotypes were measured centrally using Cumulus6 16 (provided by M.J.Y.), a computer-assisted method that requires the reader to select the pixel intensity threshold for distinguishing the dense and nondense portions of the breast image. The reader must also define the pectoral muscle boundary, whereas Cumulus6 detects the outer edge of the breast automatically for most FFDM images. Cumulus computes the PD as the DA divided by the total breast area. The NDA is equal to the total area minus the DA. All Cumulus measurements were performed by a single radiological technologist (R.Y.L.) trained by the software developer (M.J.Y.) and a breast imaging specialist (J.A.L.), who was certified in the Cumulus method and provided the gold standard measurements used for training and longitudinal evaluation.
In the first study, 23 batches of up to 1100 Hologic images were read consecutively over a period of eight months. Each batch contained 10% quality control images, including random replicates used to assess reader reproducibility and images with gold standard measurements by JAL used for calibration and periodic retraining. The batch-adjusted Pearson R for PD, DA, and NDA were: 0.952, 0.925, and 0.996 in the Hologic cohort. Following the completion of all Hologic density assessments, the reader underwent a training period to attain high reproducibility on GE images. Then, in the second study, 6 batches of up to 700 GE images were read consecutively over a period of 3 months. The batch-adjusted Pearson R for PD, DA, and NDA were: 0.961, 0.941, and 0.995 in the GE cohort. When multiple measurements were obtained per image, the mean values were used in subsequent analyses.
Genotyping, quality control, and imputation. Over 650,000 SNPs were genotyped at the UCSF Institute for Human Genetics, Genomics Core Facility using a custom Affymetrix array optimized for individuals of European ancestry 48,52 . This array is estimated to report on 93% of common variants with minor allele frequency (MAF) > 0.05, and 73% of less common variants (MAF between 0.01 and 0.05) at r 2 > 0.80, based on the 1000 Genomes Project (http://1000genomes.org) European population. Arrays were processed using the Affymetrix Axiom reagent kit 1.0 (96.7%) or 2.0 (3.3%). Genotype quality control procedures have been described previously 49 . Samples were excluded if the genotyping call rate was <0.97, or if there was evidence of trisomy, monosomy, male or ambiguous sex (PLINK v1.07 53 X chromosome F-statistic >0.2), or excess heterozygosity (PLINK 53 F-statistic < −0.03). Among first-degree female relatives (457 pairs and 14 trios), only the youngest woman was retained for analysis. Principal components of ancestry were computed for the genotyped SNPs using EIGENSOFT4.2 54 , and women were excluded if their principal components were not consistent with European ancestry 48 . Over 30 million variants were imputed from the 1000 Genomes Project reference panel using IMPUTE2.2.2 55-57 , after pre-phasing the genotyped SNPs using SHAPEIT v2.r644 58 . After excluding variants with MAF < 0.01 or imputation r 2 ≤ 0.3 59 , there remained 9,906,178 variants available for analysis.
GWAS meta-analysis. The MD phenotypes for each study cohort were transformed separately to attain standard normal distributions with mean 0 and variance 1, to facilitate estimation of the combined meta-analytic effects and enable interpretation of effect sizes in SD units. The distributions of MD phenotypes differed between women in the two cohorts and required different transformations ( Supplementary Fig. 1). The optimal power transformations determined using the R boxcox package for DA, NDA, and PD, respectively, were: fifth-root, cube-root, and cube-root for the Hologic cohort; and cube-root, cube-root, and square-root for the GE cohort. For computational efficiency, each phenotype was pre-adjusted for image batch using linear regression, and the residuals were used in GWAS analyses.
Separate GWAS analyses of each standardized MD phenotype in the Hologic cohort (n = 20,311) and GE cohort (n = 3881) were performed with PLINK 53 using linear regression models of each SNP as an additive dosage effect 60 , adjusted for ln(BMI), age at mammography, genotyping reagent kit, and the first ten principal components of European ancestry 48,49 . BMI was determined from electronic health records for the patient visit closest to the date of mammography. There was a linear relationship of age ( Supplementary Fig. 6) and ln(BMI) with the normalized density phenotypes, except in the extreme tails of the BMI distribution where the data were sparse (Supplementary Fig. 7). GWAS meta-analyses were conducted using an inverse-variance weighted fixed-effects model implemented in METAL 61 . To be considered statistically significant, we required SNPs to: meet the conventional genome-wide significance threshold of P < 5 × 10 −8 in the metaanalysis of the Hologic and GE studies combined; and have the same direction of association in both studies. A single genotyped SNP rs3819405 on 6p22.3 that was significantly associated with DA and PD had low levels of LD with all nearby SNPs (Supplementary Figs. 3H and 5F); rs3819405 had a high call rate of 99.89%, and the MAF of 0.33 was similar to the MAF of 0.34 among individuals of European ancestry in the 1000 Genomes Project.
Conditionally independent SNPs within the same chromosomal region were identified by conditional analyses adjusting for the lead SNP. To be considered statistically significant, both SNPs were required to have: meta-analytic P < 5 × 10 −8 in the conditional analysis; the same direction of associations in both studies; and low linkage disequilibrium (LD; r 2 < 0.10). Novel loci were identified by conditional analyses adjusting for the nearest known SNP for any MD phenotype. To be considered novel, SNPs were required to meet the genome-wide significance threshold of P < 5 × 10 −8 in conditional analyses adjusting for the nearest known MD SNP, and to have low LD (r 2 < 0.10) with previously reported genome-wide significant SNPs for any MD phenotype on the same chromosome.
Quantile-quantile plots and genomic inflation factors 62 were used to assess the presence of inflated significance levels due to uncontrolled population substructure. LocusZoom v1.3 63 plots of the 400 kb region centered around each novel lead SNP were used to visualize the GWAS meta-analysis significance levels, linkage disequilibrium with the lead SNP, local recombination rates from HapMap, and nearby genes.
Associations with breast cancer. We evaluated associations of newly identified MD SNPs with breast cancer risk in 122,977 cases and 105,974 controls of European ancestry from the BCAC and Discovery, Biology and Risk of Inherited Variants in Breast Cancer Consortium (DRIVE) 15 . Associations of MD SNPs with risk of estrogen receptor (ER)-negative breast cancer were evaluated in a subset of 21,468 cases and 100,594 controls 15,64 . Summary statistics were obtained from: http://bcac.ccge.medschl.cam.ac.uk/bcacdata/. MD SNPs were considered to be potentially novel breast cancer loci if they were associated with breast cancer with P < 0.05 and P > 5 × 10 −8 in the BCAC/DRIVE data, and were uncorrelated (r 2 ≤ 0.01) with previously reported genome-wide significant breast cancer SNP in the NHGRI-EBI GWAS catalog 65 (https://www.ebi.ac.uk/gwas/).

Mendelian randomization and genetic correlation of MD and breast cancer.
Mendelian randomization analyses were performed to evaluate the potential causal associations of MD phenotypes with breast cancer risk. We used the weighted median method to estimate the causal effect because it is more efficient and robust to violations of instrumental variable assumptions than other Mendelian randomization methods for summary statistics, and provides consistent estimates even when up to half of the information comes from invalid instrumental variables 18 . We considered the first reported SNP at all independent prior and new genomewide significant loci (Table 1 and Supplementary Table 7) for DA (n = 28), NDA (n = 17), and PD (n = 20), and estimated their associations with the relevant MD phenotype in the GWAS meta-analysis. Summary statistics for SNP associations with breast cancer in 122,977 cases and 105,974 controls of European ancestry were obtained from BCAC/DRIVE 15 .
LD Score regression was performed to estimate the genetic correlation between MD phenotypes and breast cancer from GWAS summary statistics using the LDSC v1.0.1 software 20,21 . We used the LD scores for the European ancestry population from the 1000 Genomes Project provided by the software developers 20,21 . A total of 779,828 SNPs were included that had available LD scores and summary statistics from this MD GWAS and the BCAC/DRIVE breast cancer GWAS 15 .
Regulatory function of MD SNPs. We evaluated whether MD SNPs were associated with gene expression levels in: human mammary tissue (n = 251), primary fibroblast cells (n = 300), subcutaneous fat (n = 385), visceral fat (n = 313), and whole blood cells (n = 369) using data from the Genome-Tissue Expression (GTEx) project version 7 (https://gtexportal.org/) 22,23 . Significant cis-eQTLs within 1 Mb of the gene transcription start site were identified by computing q-values for SNP-gene pairs involving one of the lead MD SNPs using the R qvalue package, and controlling for a false discovery rate (FDR) of 0.05. We also assessed whether lead MD SNPs or nearby proxies (r 2 ≥ 0.80 in Europeans from the 1000 Genomes Project) were located within promoter or enhancer regions in human primary mammary epithelial cells and human primary mammary fibroblasts using data from the ENCODE and Roadmap Epigenomic consortia (https://www. encodeproject.org/) 24,25 . Promoter-like regions were identified by combining DNase hypersensitivity and histone modification H3K4me3 signals in the same cell type 24,25 . Enhancer-like regions were identified based on DNase hypersensitivity and histone modification H3K27ac signals in mammary epithelial cells, and DNAse hypersensitivity only in mammary fibroblasts 24,25 .
Tissue enrichment of MD SNPs in regulatory regions. We tested whether previously identified and new loci for each MD phenotype were enriched in regulatory regions in 125 diverse human cell and tissue types using the UES method (https://github.com/robertkleinlab/uesEnrichment) 26 . Open chromatin regions were defined using DNase I hypersensitive sites sequencing data from the ENCODE and Roadmap Epigenomic consortia 27 . The empirical P value for the observed enrichment of independent MD SNPs in the regulatory regions of each cell line was computed by generating 10,000 sets of randomly selected SNPs that were matched to MD SNPs by the distance to the nearest transcription start site and number of correlated SNPs. We applied a Bonferroni correction for the 125 cell lines tested to determine the significance threshold of P < 0.0004. Pathway analysis. We performed gene set enrichment analyses to identify biological pathways implicated by prior and new MD loci using DAVID v6.8 28,29 . For each MD locus, we included the nearest flanking protein-coding genes within 500 kb of the lead SNP, as well as target genes whose expression levels were associated with the lead MD SNP (cis-eQTL) in mammary tissue, primary fibroblast cells, subcutaneous fat, visceral fat, or whole blood (Table 1 and Supplementary Tables 3,  7, and 8). We estimated the fold enrichment of MD genes among gene sets or pathways with FDR < 0.05 from the Gene Ontology, KEGG, Reactome, and Biocarta databases, and identified the responsible genes using DAVID 28,29 .
Heritability. We estimated the proportion of phenotypic variance explained by the additive genetic effects of all genotyped SNPs using GCTA v1.02 31 . Array-based heritability was estimated separately for the Hologic and GE study cohorts, and the resulting estimates were then combined using inverse-variance weighting. We also estimated the proportion of phenotypic variance explained by all independent prior and new genome-wide significant MD loci. For each MD phenotype, the residual variance of linear regression models were estimated using 100-fold cross-validation, where model 1 included all independent genome-wide significant loci in addition to the adjustment variables in the GWAS model, and model 2 included only the non-SNP covariates. In both models, the SNP effects and covariates were nested within the Hologic or GE studies to account for differences across studies. The proportion of variance explained by the genome-wide significant loci for each MD phenotype was computed by 1 − V1/V2, where V1 and V2 represent the estimated residual variances of models 1 and 2 respectively.
Reporting summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Data availability
Genotype data of RPGEH GERA participants are available from the database of Genotypes and Phenotypes (dbGaP) under accession phs000674.v3.p3. This includes individuals who consented to having their data shared with dbGaP. The complete GERA data are available upon application to the KP Research Bank (https://researchbank. kaiserpermanente.org/our-research/for-researchers). Breast cancer summary statistics are available at http://bcac.ccge.medschl.cam.ac.uk/bcacdata/. All remaining relevant data are available in the article, Supplementary Information, or from the corresponding author upon reasonable request.