High-Resolution Bisulfite-Sequencing of Peripheral Blood DNA Methylation in Early-Onset and Familial Risk Breast Cancer Patients

Breast cancer is the leading type of cancer in women worldwide, representing nearly 25% of all cases (1). Approximately 15% of invasive breast cancer can be traced to a hereditary origin, with patients having at least one first-degree female relative with the disease. The landmark discoveries of the first two major breast cancer susceptibility genes, Breast cancer 1 (BRCA1; ref. 2) and Breast cancer 2 (BRCA2; ref. 3), over 20 years ago revealed the power of genetic testing; over 60% of women with heterozygous germline mutations in either of these genes ultimately develop breast cancer by the age of 70. However, there have been diminishing returns with each subsequent tumor suppressor identified. Whereas mutations in BRCA1 and BRCA2 together account for around 15% of hereditary cases, deleterious genetic alterations in other genes such as Cadherin-1 (CDH1), PTEN, and Partner and localizer of BRCA2 (PALB2) account for less than 7% of hereditary cases combined (4). After taking into account the cumulative and putative effects of SNPs on this disease, 50% of hereditary breast cancer may still remain unexplained from a genetic perspective (4).

Tumor suppressor inactivation by hypermethylation has been well described, mainly as a “second hit” in tumors (5). Because sporadic breast cancer tumors with BRCA1 hypermethylation remain histologically and molecularly similar to hereditary breast cancer tumors with BRCA1 mutations (6, 7), we investigated whether somatic epigenetic alterations could exist more broadly in familial breast cancer. This hypothesis follows multiple epigenome-wide association studies that have used methylation arrays to suggest that heritable methylated marks can contribute to hereditary cancers, including breast cancer (8–11). Controversy has surrounded the potential for transgenerational epigenetic inheritance of methylated alleles that contribute to disease ever since the concept has been postulated (12, 13). While the bar for irrefutable proof is beyond the experimental conditions that can be performed in human studies, observational support for this phenomenon has come from heritable hypermethylated MutL homolog 1 (MLH1) alleles that lead to Lynch syndrome (14, 15), a familial predisposition to colorectal cancer (16). More recently, a study of 49 families looking for BRCA1 epimutations has identified two instances of a dominantly inherited, constitutively methylated allele at the BRCA1 promoter that leads to breast cancer (17).

Despite establishing the contribution of heritable BRCA1 epimutations to breast cancer, the prominence of this occurrence across the genome remains unknown. Here we utilize an unbiased approach to investigate breast cancer predisposition by DNA methylation. Using a cohort of 99 patients with familial and early-onset breast or ovarian cancer, we assessed putative germline DNA methylation from whole blood using bisulfite sequencing at known gene promoters and CpG islands. Although we did not observe promoter epialterations at BRCA1 or BRCA2 in patients with cancer, we did find other instances of allelic methylation with the potential to contribute to disease in 9% of our cohort. This includes allelic hypermethylation at known tumor suppressor loci such as PTEN as well as allelic hypomethylation at known proto-oncogenic loci such as Testis expressed 14 (TEX14). We also identify sites of recurrent allelic methylation specific to patients with breast cancer, relative to a control cohort of unaffected women. These results, which are independent from chemotherapy treatment, provide proof of principle that allelic DNA methylation is also frequently present in individuals with cancer predisposition. Further studiesare warranted to determine whether systemic DNA methylation may able to help explain hereditary breast cancer with unknown origins.

The cohort was selected from patients who visited Stanford Cancer Genetics Clinic between 2004 and 2013. Out of a total of 333 patients, 248 were affected. Samples were classified by the likelihood of hereditary breast cancer predisposition based on the age of first cancer diagnosis, number of cancer diagnoses, type of cancer, and family history of breast or ovarian cancer (Supplementary Fig. S1). Ninety-nine affected patients were selected for this study, of which 87 had history of breast cancer only, 6 were patients with ovarian cancer and 4 had both breast and ovarian cancers, and 2 had other cancers with history of breast and ovarian cancers. Fifty-nine percent of the cohort had family history of breast or ovarian cancers. Earlier cancer onset, multiple cancers in the same individual, and greater family history of cancer were indicators of hereditary disease.

All of the patients received clinical multigene panel testing. Seventy-six patients had variants of uncertain significance or benign alterations, and 23 patients had pathogenic or likely pathogenic mutation in one of these genes: Ataxia-telangiectasia mutated (ATM), BRCA1, BRCA2, Checkpoint kinase 2 (CHEK2), MLH1, MutS homolog 2 (MSH2), MutY DNA glycosylase (MUTYH), PALB2, PMS1 homolog 2 (PMS2), or RAD51 homolog C (RAD51C). Thirteen of these patients were specific to BRCA1 or BRCA2 mutations. Patients were classified in three different categories according their exposure to chemotherapy treatment: active treatment (including patients who had completed their treatment within two months of collection), prior treatment, and never treated with chemotherapy (Table 1).

The control cohort was recruited from nonsmoking female Stanford employees (graduate students, laboratory personnel, nurses, physicians, research scientists, medical assistants, residents, and fellows) on a volunteer basis as described previously (18). Healthy controls were women between the ages of 30 and 65 years and were unaffected by breast cancer and had no personal or first-degree relative with a history of any type of cancer except for basal cell carcinoma. Only data about age, race, family history, and menopausal status were available for the healthy controls (Table 2). There were no significant statistical differences between the cohorts based on age of disease onset and healthy controls (P > 0.05, Student t test), race (χ² P > 0.05), or menopausal status (χ² P > 0.05, combining premenopause and treatment-induced menopause). However, there was a significant statistical difference between the age of collection between both groups (P = 7.7 × 10⁻⁵, Student t test) and menopausal status when excluding treatment-induced menopause (χ² P = 0.023).

This study was carried out with informed consent from all the participating individuals. The study was approved by the Institutional Review Board (IRB) of Stanford University (Stanford, CA) and was conducted in accordance with the approved methods and guidelines of the U.S. Common Rule.

Genomic DNA was extracted from 3 mL peripheral whole blood using the Gentra Puregene Blood Kit (Qiagen Sciences) according to the manufacturer's protocol.

Sequencing libraries were constructed using the MethylSeq^XT kit (Agilent Technologies) according to manufacturer's instructions using 1 μg of gDNA. This target enrichment system captures 84 Mb of the genome focused on known cancer and developmental differentially methylated regions (DMR) and CpG islands and provided a cost-effective means to sequence at a high depth-of-coverage (averaging 71× depth across all samples), providing highly quantitative methylation data at base-pair resolution (Supplementary Fig. S2A and S2B).

Briefly, sheared DNA was repaired, A-tailed, and ligated to adaptors prior to overnight hybridization with biotinylated baits. Target enrichment was performed using Streptavidin T1 magnetic beads (Life Technologies), and the purified DNA was bisulfite converted and desulphonated using the EZ DNA Methylation Gold kit (Zymo Technologies). Inserts were PCR amplified and indexed using 8-bp methylation barcodes. All DNA purification steps were performed using AMPureXP Magnetic Beads (Beckman Coulter) and the libraries were sequenced on an Illumina HiSeq4000 at the Stanford Center for Genomics and Personalized Medicine.

Bisulfite HTS data were analyzed on the Stanford's SCG compute cluster. FASTQ files were trimmed and mapped to hg19 using BSMAP (v. 2.89; ref. 19). PCR duplicates were removed using Picard Tools. Reads mapping to the reference sense strand were filtered out to minimize off-target, spurious results to the (+)-stranded target enrichment baits. Methylation ratios, CpG read coverage, and estimated bisulfite conversion were all calculated using MOABS (v 1.3.0; ref. 20). Publicly available whole-genome bisulfite sequencing (WGBS) data for three breast tumors, BT089, BT126, and BT198, were obtained as raw FASTQ files from ArrayExpress under accession E-MTAB-2014 and processed in the same pipeline, with reads accepted from both strands.

Differential methylation analysis was performed using Metilene (v. 0.2-6; ref. 21), using the provided input script and filtering for at least 5 reads per CpG. DMRs were filtered for an absolute difference of 0.1 and at least 3 CpGs. Quantifying manually observed DMRs was performed by metilene analysis of a subset of the genome containing the region of interest comparing the patient against the healthy control cohort. P values reported are products of Mann–Whitney U tests. Unbiased DMR screening in individuals was performed by iterating through patients compared with the healthy control cohort. Candidate regions were identified where the healthy control cohort averaged >90% methylation or <10% methylation and the affected patient averaged 50% ± 10% methylation. Regions were filtered statistically for FDR values <0.05. Candidate allelic methylation regions in the healthy controls were identified similarly by iterating across each healthy individual compared with the remaining 48 control cohort. As a final quality control check, potential DMRs were inspected visually to filter out regions that were subject to artificial read saturation introduced by the target enrichment platform.

CpG methylation status of individual reads was assessed using the tanghulu function of cgmaptools (0.0.1; ref. 22) using samtools (v. 1.4). Kernel densities were calculated and plotted in R (v 3.3.1) using the density function.

Methylation boxplots were calculated by taking the mean of methylation ratios across a window starting 1.5 kb upstream of the most 5′-TSS and ending 500 bp downstream of the most 3′-TSS (if any other proximal, alternative TSS exist). Boxplots were displayed in R using the ggplot2 package. Statistical significance was determined with ANOVA, using the aov and TukeyHSD functions. FDR values were obtained using the p.adjust function (Bonferroni).

Methylation heatmaps were constructed by creating TDF files using IGVTools (v. 2.3.3). Genomic regions were visualized using the IGV Browser (v. 2.3.90) and exported.

Enrichment of cancer predisposition genes was determined using the hypergeometric test using the phyper function in R. The significance of the overlap of genes (114) and gene-associated epialterations (184) was determined against the genes included in the MethylSeq^XT panel (22,160).

Principal component analysis was performed in R using the prcomp function on CpGs that had a SD > 0.2 across all samples (n = 7,587). CpGs were required to have at least 5 reads to be included. PCA results were displayed using the ggplot2 and ggfortify packages.

TCGA level 3 data was downloaded from the Broad Institute's Firebrowse GDAC (v. 20160128; ref. 23) and analyzed in R. Log₂-transformed RSEM values are presented with a Welch P value.

All genomic positions reported are referenced to hg19.

Raw sequencing data files are available through dbGaP under data access control guidelines in accordance with informed consent under the study accession phs001699.

We selected a cohort of 99 patients with early-onset breast or ovarian cancer who visited the Stanford Cancer Genetics Clinic between 2004 and 2013. Selection was designed to capture samples with the greatest likelihood of hereditary breast cancer predisposition (see Materials and Methods). These patients participated in clinical genetic screening across a panel of breast cancer genes and could be divided into two subgroups: 76 patients without any known deleterious mutations in cancer-related genes (affected noncarriers), and 23 patients with known or likely pathogenic mutations in cancer-related genes such as ATM, BRCA1, BRCA2, CHEK2, MLH1, MSH2, MUTYH, PALB2, PMS2, and RAD51C (affected carriers). We also included 6 healthy BRCA1 or BRCA2 carriers who have never had cancer (unaffected carriers), and 49 healthy individuals who never had cancer served (except for basal cell cancer) as healthy controls (Table 1).

The most common cancer type was invasive ductal breast cancer with ER⁺, PR⁺, and Her2-negative receptor status (16%). The mean age of the cancer diagnosis was 41 years, varying between 18 and 61 years. Thirty-three percent of the patients had 2 or more cancers. Patients were mostly Caucasian (61%) with 2 children (34%). Their BMI was normal at the time of the blood collection (mean of 25.7), and 72% were either premenopausal or treatment had induced the menopause. Patients had varying exposures to chemotherapy treatment: 39 patients had previously completed chemotherapy and 20 patients were undergoing active chemotherapy at the time of collection (Tables 1 and 2).

Using genomic DNA extracted from whole-blood, we performed high-depth bisulfite sequencing across selected regions of the genome focused on known cancer and developmental differentially methylated regions (DMR) and CpG islands. As opposed to dramatic shifts in global methylation profiles in studies that have compared tumor tissue to normal tissue (24), in the blood DNA we observed only subtle differences in global methylation between affected patients and unaffected patients or healthy controls. Similarly, patients with known mutations in cancer-related genes did not have major differences when compared with patients without cancer. Global peripheral blood methylation ratios remained similar across the cohort (Supplementary Fig. S2C). Principal component analysis did not clearly cluster these subgroups, and only 5.615% and 1.677% of the variance could be accounted for by principal components 1 and 2, respectively (Supplementary Fig. S2D). Despite interindividual differences in methylation ratios across the genome, we did not identify any statistically significant DMRs that passed quality control between healthy controls and affected patients or unaffected carriers.

Previous studies employing probe-based and array-based strategies have observed low-level hypermethylation at the BRCA1 promoter (9, 10, 25), suggesting that germline DNA methylation at this locus could contribute to breast cancer predisposition. Coupled with the fact that inactivating mutations in either BRCA1 or BRCA2 are found in 15% of hereditary breast and ovarian cancer, representing the greatest risk factors for disease (4, 26–28), and that promoter hypermethylation at BRCA1 in breast tumors is associated with poor outcome (29), we examined whether similar promoter hypermethylation existed in our cohort. We examined these promoters using approximately 2 kb windows encompassing the transcriptional start site (TSS, 1.5-kb upstream and 500-bp downstream). This led to the identification of only modest changes (<2%) in average methylation state across these regulatory regions (Fig. 1A). Because one of the main advantages of our study is the increased resolution afforded by bisulfite sequencing, we looked directly at the methylation states for individual CpGs surrounding the TSS of BRCA1 and BRCA2. At both of these loci, methylation ratios remained similar across all groups (Fig. 1B and C).

We extended this approach to examine approximate 2 kb windows encompassing the TSS of other genes known to be frequently inactivated in hereditary breast and ovarian cancer (30). This panel of 24 additional genes includes BRCA1-A complex subunit (ABRAXAS1), ATM, BRCA1-associated RING domain 1 (BARD1), Bloom syndrome protein (BLM), BRCA1-interacting protein 1 (BRIP1), CDH1, CHEK2, FA complementation group C (FANCC), FA complementation group M (FANCM), MLH1, Meiotic recombination 11 homolog (MRE11), MSH2, Nibrin (NBN), Neurofibromin (NF1), PALB2, PMS2, PTEN, RAD51 homolog B (RAD51B), RAD51C, RAD51 homolog D (RAD51D), RecQ protein-like (RECQL), RAD50 interactor 1 (RINT1), Serine/threonine kinase 11 (STK11), and Tumor protein p53 (TP53). Akin to our observations at BRCA1 and BRCA2, the average methylation ratios at these gene promoters were similar across the cohort, with few significant differences between subgroups (Supplementary Fig. S3).

Despite the statistically significant changes in methylation ratios between healthy controls and an affected subgroup at ABRAXAS1, MSH2, STK11, and TP53, a closer inspection at base-pair resolution at these loci indicated that these results were likely due to noise at the terminal end of target enrichment regions and low read coverage, and not reflective of true biological differences (Supplementary Fig. S4A–S4D). Similarly, the average hypermethylation observed at the RAD51C promoter in both affected cancer cohorts when compared with healthy controls was not driven by any obvious changes at the individual CpG or sample level. Instead, we identified a significant DMR (chr17:56,769,048-56,769,761; P = 4.4 × 10⁻¹⁴) upstream of the 5′UTR in patient 2417 that averaged 56.9% methylation while all other samples displayed nearly complete methylation (Fig. 2A, arrow). All healthy controls averaged 90.7% methylation within this region. As a tumor suppressor involved in the DNA damage response pathway, canonical gene activation of RAD51C by hypomethylation seemed to be a counter-intuitive mechanism for predisposition and oncogenesis. However, RAD51C shares a bidirectional promoter with TEX14, and we found that this hypomethylated region completely encompassed the TSS and first exon of TEX14. TEX14 cooperates with SCY1-Like 1 (SCYL1) and Polo-like kinase 1 (PLK1) to degrade the RE-1 silencing transcription factor (REST) tumor suppressor and is frequently overexpressed in triple-negative breast cancer (31). We observed that the majority of sequencing reads within this locus were either completely methylated or entirely unmethylated (Fig. 2B), suggesting the presence of a heterozygous unmethylated allele with the potential for aberrant, overexpression of TEX14. Notably, patient 2417 presented clinically with triple-negative breast cancer, and the time of collection was 8 years posttreatment with chemotherapy.

We assessed for additional DNA methylation differences in the form of aberrant, allelic methylation in this panel of known genes in hereditary breast and ovarian cancer. We identified a similar DMR within the 5′-UTR of PTEN in patient 2282 (chr10:89,623,973-89,624,168; P = 6.0 × 10⁻⁶) that averaged 30.9% methylation, whereas most other samples were devoid of methylation (Fig. 1C, arrow); this patient also had a heterozygous PALB2 mutation. Healthy controls averaged 1.39% methylation within this region. Patient 2282 was diagnosed with ductal breast cancer at 39 years of age and was actively receiving chemotherapy at the time of collection. PTEN encodes for a tumor suppressor whose primary function is to regulate the PI3K-Serine/threonine Kinase 1 (AKT) pathway. Inactivating mutations in PTEN are frequently found in Cowden syndrome, an inherited disorder characterized by pathway hyperactivity that promotes cell growth and survival, leading to predisposition to several types of cancer, including breast (32). In addition to altered PTEN methylation in Cowden syndrome (33), promoter methylation at PTEN is commonly found in breast tumors, where it associates with lower gene expression (34). Similar to what we observed at the TEX14 promoter in patient 2417, we found that the majority of sequencing reads at this locus were either completely methylated or entirely unmethylated (Fig. 2D), indicating a heterozygous methylated allele with the potential to inhibit PTEN expression and drive a haploinsufficient phenotype. Seven additional patients with cancer exhibited promoter hypermethylation at the PTEN promoter, but to a lesser degree (patients 513, 2203, 2208, 2217, 2221, 2260, and 5111; Fig. 2C, arrowheads). When analyzed in aggregate, these samples share the same DMR as patient 2282 (chr10:89,623,973-89,624,168; P = 4.4 × 10⁻¹⁴), but average only 18.2% methylation and have methylated:unmethylated allele ratios that range from 0.12–0.43 (Supplementary Fig. S5) compared with 0.88 in patient 2282. Three of these patients were heterozygous carriers for BRCA1 or CHEK2 mutations (patients 513, 2203, 5111); there were no commonalities across these tumors based on type or hormone receptor status. Patients 513, 2217, 2221, and 2260 were never treated with chemotherapy by the time of collection. Samples from patients 2203, 2208, and 5111 were collected at 2 years, 6 years, and 8 years postchemotherapy, respectively, but each of these patients were not on chemotherapy at the time of collection.

We also observed an intragenic DMR within CDH1 in patient 2980 (chr16:68,856,954-68,857,037; P = 1.3 × 10⁻⁴) that suggested allelic methylation (Fig. 3A and B). Healthy controls averaged 94.9% methylation, whereas patient 2980 averaged 61.8% across the region. This patient is a BRCA2 carrier and had triple-negative invasive ductal breast cancer when she was 38 years old and was actively receiving chemotherapy at the time of collection. CDH1 encodes for E-cadherin, a cellular adhesion molecule that cooperates with the β-catenin/Wnt signaling pathway, and is a common breast cancer tumor suppressor (35). Germline CDH1 mutations confer a high risk for signet ring cell diffuse gastric cancer and lobular breast cancer (36). This region of allelic methylation is located upstream of an alternative start site encoding a transcriptional variant lacking the extracellular domain.

Our initial results indicate that the most informative epigenetic changes contributing to breast cancer predisposition may arise as allelic methylation. We next turned to an unbiased approach to identify potentially similar epialterations by searching for DMRs in individual affected patients compared with healthy controls. Potential sites for aberrant allelic gain or loss of methylation were established by searching for uniformly methylated or unmethylated regions across all healthy controls that then approach 50% methylation in individual patients with cancer, at an FDR < 0.05. Candidate epialterations were rare, as 68% of samples had one or no regions identified, and only 9% of samples had five or more regions identified in our screen (Supplementary Tables T1 and T2). We identified an enrichment of candidate epialterations in cancer predisposition genes in the cancer cohort: 4 instances of candidate epialterations in the cancer cohort overlapped with a curated list of 114 cancer predisposition genes (37), whereas no candidate epialterations in the healthy controls were found in the list (Fig. 4A; P = 2.65 × 10⁻³).

The unbiased nature of this screen allowed us to identify patients with other epialterations in regions of interest. Among these include patient 2203 who had a hypermethylated DMR that encompassed miRNA-596 (chr8:1,765,429-1,765,789; P = 1.6 × 10⁻⁹; Fig. 4B). This patient averaged 40.1% methylation across this region, whereas all healthy controls averaged just 4.01% methylation, and methylation read densities suggested this methylated region was allelic (Fig. 4C). DNA methylation can also regulate miRNA expression, and disruption of this process leads to functional consequences in tumorigenesis and metastasis (38). miR-596 has been previously shown to induce p53-mediated apoptosis (39) and has tumor-suppressive functions when expressed in melanoma and oral squamous cell carcinoma (40, 41). Patient 2203 was previously mentioned to have increased PTEN promoter methylation and is a BRCA1 carrier; her sample was collected 2 years postchemotherapy.

We also observed a DMR near the TSS of TatD DNase domain containing 1 (TATDN1 - chr8:125,550,351-125,551,216; P = 1.3 × 10⁻¹²) in patient 14246 that extended through the first exon of the neighboring gene NADH:ubiquinone oxidoreductase subunit B9 (NDUFB9) upon visual inspection (Fig. 5A). Patient 14246 averaged 43.4% methylation over this region, whereas healthy controls averaged 6.6% methylation. This patient had ductal breast cancer at 33 years of age and never received chemotherapy by the time of collection. TATDN1 amplifications and fusion products with Gasdermin B (GSDMB) have been reported in breast carcinomas (42, 43), and NDUFB9 knockdown has been shown to promote proliferation in breast cancer cells by disrupting mitochondrial metabolism (44). The hypermethylation in this individual appeared to be allelic based on the methylation status of sequencing reads strictly over the region encompassing TATDN1 (Fig. 5B), whereas the region overlapping NDUFB9 has a broader density of methylated reads, contraindicating allelic methylation (Fig. 5C). Allelic methylation overlapping TATDN1 was confirmed by the presence of a heterozygous SNP within the DMR (chr8:125,550,868; C>G; rs7009290). We identified 128 reads that covered this specific SNP, with all 64 methylated reads associated with the polymorphic guanine allele and all 64 unmethylated reads associated with the reference cytosine (Fig. 5D). Because only 14 bp separates the TSS of TATDN1 and NDUFB9, this methylated allele could function to affect transcription of either or both genes.

Our screen also revealed recurrent hypomethylation at the insulin-like growth factor–binding protein, acid-labile subunit (IGFALS) promoter. IGFALS cooperates with insulin growth factor (IGF)-binding proteins and drastically enhances the circulating half-lives of IGF-I or IGF-II when combined in a ternary complex (45). GWA studies have identified several SNPs within the IGF pathway that are associated with breast cancer (46, 47), and high IGF concentrations are associated with increased cancer risk, including breast (48, 49). We observed this allelic methylation (chr16:1,844,768-1,845,114, P = 5.0 × 10⁻¹⁴) in 4 patients with breast cancer who averaged 59.3% methylation (2981, 6305, 8192, and 16251; Fig. 6A, arrows), but not in healthy controls where methylation ratios averaged 93.1%. One of these patients (2981) was a carrier for a BRCA1 mutation. Both patients 2981 and 8192 were actively receiving chemotherapy at the time of collection, whereas patients 6305 and 16251 never received chemotherapy by the time of collection. This region encompassed the first exon of an alternative isoform of IGFALS and overlapped with a hypomethylated region in three breast tumors with publicly available WGBS data (50). Sequencing reads indicated allelic methylation in 4 patients in our cohort (Fig. 6B), whereas the hypomethylation observed in tumor WGBS was more extensive and complete in some areas. All 4 of these patients had a heterozygous SNP within the DMR (chr16:1,844,801; T>C; rs9924224) where unmethylated reads were associated with the polymorphic cytosine allele and methylated reads were associated with the reference thymine allele (Fig. 6C–F). In patient 2981, this effect was consistent over 37 reads (24 unmethylated reads overlapping the cytosine and 13 methylated reads overlapping the thymine). This effect was also consistent in patient 8192 over 25 reads (9 unmethylated reads with cytosine and 16 methylated reads with thymine). Patient 16251 had 27 total reads overlapping this region with 20 unmethylated reads overlapping the cytosine and 7 methylated reads overlapping the thymine. We observed more mixed methylation in patient 6305 where 4 methylated reads were associated with the polymorphic cytosine (25% of reads) and one unmethylated read associated with the thymine (5.3% of reads). Average IGFALS expression is also significantly increased in TCGA breast tumors compared with normal tissue (RSEM+1 5.60 vs. 3.78, P < 2.2 × 10⁻¹⁶; Fig. 6G).

DNA methylation can have profound effects on transcription and splicing (51, 52), and deregulation of both of these processes are hallmarks in the initiation of cancer (53, 54). Despite clear changes in breast tumor methylomes (10, 55), the role of methylation in cancer predisposition has remained less clear, with some proposing that epigenetic changes in normal body cells are indicative of global mechanistic disruption (56). Using high depth-of-coverage bisulfite sequencing across the genome, we provide evidence for epimutations in the form of allelic methylation at TEX14, PTEN, CDH1, IGFALS, miR-596, and TATDN1/NDUFB9, comprising 9% of our entire cohort (6.5% and 17% of affected noncarriers and affected carriers, respectively). Our methylation study is the first of its kind to definitively identify both hypermethylated and hypomethylated alleles at known tumor suppressor and proto-oncogene loci in peripheral blood, respectively. Bisulfite sequencing at base-pair resolution should remain the gold standard to identify allelic methylation changes.

Our findings did not detect altered DNA methylation at the BRCA1 promoter, whereas previous studies targeted at this locus have identified slightly elevated promoter hypermethylation in germline DNA (9, 10, 25). Although we identified statistically significant changes in methylation ratios across the TSS windows, the fact that they that do not hold up to further scrutiny at base-pair resolution suggests that probe-specific assays could be subject to significant noise or false positives. Our deep sequencing approach is also not well-suited to detect methylated mosaicism at very low frequencies (i.e., <2%), which would have unknown biological impact and could be attributed to numerous sources. Although the recent discovery of heritable BRCA1-methylated alleles in breast cancer was found in 2 of 49 families (17), the size of our cohort suggests that it could exist at much lower rates in the general population.

Because all patients with cancer in this cohort, regardless of carrier mutation status, were selected for early onset breast cancer, the cooccurrence of mutations in both the genome and epigenome of some patients could raise the possibility of epialleles as modifiers of breast cancer risk in the setting of known pathogenic high-risk mutations. Because of our sample size, our data remains insufficient for any predictive cancer-risk assessment of these epialleles in such confounded settings. It remains possible that the altered methylated regions in these patients represent secondary events, or acquired changes in systemic DNA methylation. This could be particularly true for the additional 7 patients who we identify to have increased PTEN promoter hypermethylation at levels that are inconsistent with allelic methylation. A potential mechanism behind this possibility remains unclear, as each of the three carrier cases have mutations in different genes (BRCA1, BRCA2, and PALB2). Also, because we identify candidate epialleles in patients with various chemotherapy regimens, including those who never had chemotherapy by the time of collection, we rule out the potential for therapy-related changes in immune cell composition or systemic therapy–related methylation changes.

Genetic control of methylation, or methylation QTL (meQTL), may exist for 15% of methylated regions (57, 58). Previous studies have also identified meQTL in blood that coincide with breast cancer GWAS SNPs in the NHGRI catalog (59), inviting speculation that these SNPs may reflect heritable epivariants such as the ones we identify here. Because cis-meQTL have been identified at distances from 2–10 kb and even up to 1 Mb from the methylated region itself, our target enrichment approach does not allow us to corroborate these previous findings due to the lack of sequencing coverage surrounding DMRs. Indeed, although we find heterozygous SNPs in patients that confirm the allelic nature of the epivariants observed at IGFALS and TATDN1, we did not observe heterozygous SNPs within other methylated alleles. It should be noted that other patients, including controls, were also heterozygous at the IGFALS and TATDN1 DMR SNPs without a change in methylation, suggesting that these are not putative meQTL.

Despite the inability to make predictive cancer-risk claims or to identify novel candidate predisposition genes, recurrent IGFALS methylation identifies an intriguing candidate that might not have been identified in a standard comparison between tumor and normal. Notably, the epialleles we observe overlap with the TSS of an alternative IGFALS isoform, and any potential functional roles for IGFALS in breast cancer will need to be studied in more detail with a more appropriate cohort and system.

Importantly, this study is based on a single-institution cohort, and thus faces several limitations that can only be addressed with multiple validation studies and larger cohorts. We did note a statistically significant difference of 6.8 mean years between the ages of collection between the patients with cancer and the healthy controls, and methylation has been previously reported to change with age (60). Although we are unable to correct for these differences in our study, we note that age-related global DNA hypomethylation can be separate from changes at discrete regions, of which we identify instances of both hypo- and hypermethylation. Also, as our cohort consists of patients spread throughout the country who had their blood collected over the span of 9 years, additional biospecimens from these individuals remain difficult or impossible to obtain. Thus, we cannot conclusively show the functional effects of these methylation alleles. This is evident in case of the intragenic hypomethylation we identify within CDH1, which could result in alternative or aberrant expression of the variant containing only cytoplasmic domains that could then affect Wnt-signaling (61). Alternatively, because gene-body methylation is associated with highly expressed genes and can function to maintain high fidelity transcription and proper splicing (52, 62), loss of methylation at this locus may result in lowered gene expression or hypomorphic transcripts.

In addition, without access to these individuals' first-degree relatives, we cannot construct pedigrees to determine if what we have identified are sporadic epimutations, or truly heritable epivariants. As we examined methylation in peripheral blood, these DMRs largely exist in white blood cells and could be attributable to acquired or environmental factors, such as immune responses to the cancer or other infection. Future studies involving familial cohorts and multiple tissue types, including tumors, will be needed to address whether these alleles are constitutional and heritable. Despite these limitations, our study suggests that more research is warranted to identify the impact of the epigenome in the predisposition and development of breast cancer.

J.J. Gruber reports receiving other commercial research support from Curis, Inc., and holds ownership interest (including patents) in Tocagen, Inc. and BD Biosciences. M.P. Snyder holds ownership interest (including patents) in Abcam, Epinomics, Personalis, SensOmics, Qbio, January, Akna, Filtricine, and Tailai, and is a consultant/advisory board member for Personalis, SensOmics, Qbio, January, Akna, Filtricine, Tailai, Genapsys, and Jupiter. No potential conflicts of interest were disclosed by the other authors.

Conception and design: J. Chen, M.K. Haanpää, J.J. Gruber, J.M. Ford, M.P. Snyder

Development of methodology: J. Chen, M.K. Haanpää

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): J. Chen, M.K. Haanpää, J.M. Ford

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): J. Chen, M.K. Haanpää, N. Jäger, J.M. Ford

Writing, review, and/or revision of the manuscript: J. Chen, M.K. Haanpää, J.J. Gruber, J.M. Ford, M.P. Snyder

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): J. Chen, M.K. Haanpää, J.J. Gruber, J.M. Ford

Study supervision: J.M. Ford, M.P. Snyder

We extend our thanks to the many women who generously participated in the study. We also thank all members of the Snyder and Ford labs for helpful discussions. This work used the Genome Sequencing Service Center by Stanford Center for Genomics and Personalized Medicine Sequencing Center, supported by the grant award NIH S10OD020141. M.K. Haanpää is supported by grants from Sigrid Juselius Foundation, Orion Research Foundation, and Päivikki ja Sakari Sohlberg Foundation. J.J. Gruber was supported by fellowships from the Jane Coffin Childs Memorial Fund for Medical Research, Stanford Cancer Institute, and Susan G. Komen Foundation, as well as funding from ASCO, the Conquer Cancer Foundation, and the Breast Cancer Research Foundation. N. Jäger was supported by an EMBO Long-Term Fellowship (ALTF 325-2014). J.M. Ford was supported by the BRCA Foundation and the Breast Cancer Research Foundation. M.P. Snyder is supported by grants from the NIH, including a Centers of Excellence in Genomic Science award (5P50HG00773504).

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.