Evolutionary histories of breast cancer and related clones

Revolutionized sequencing technologies have enabled a series of recent studies on clones evolving in apparently normal tissues, which are documented by sensitively detecting clone-defining somatic mutations^{1,2,3,4,5,6,7}. In view of cancer development, a key observation through these studies is that clonal outgrowth in normal or non-cancer tissues is common, often pervasive and frequently driven by common cancer mutations¹. This immediately points to an important implication to the early history of cancer that one or more of those positively selected non-cancer clones should be destined for subsequent cancer development^1,8. Here, among key questions that studies on normal tissues cannot answer are when cancer arises from these non-cancer clones by acquiring what additional mutations, while other clones partially sharing common mutations are still normal or precancer, and what the difference in mutation profile is between cancer clones and those non-cancer relatives. Phylogenetic analyses using multisampling of cancer specimens have been used to infer the life history of cancer in terms of driver events. However, analyses of cancer tissue alone frequently obscure the order of early driver events that are often assigned together to a long major trunk in the phylogenetic tree^9,10,11. Moreover, they do not help map the timing at which cancer clones emerged or track the fate of other related non-cancer clones. To answer these issues, analyses of both cancer and non-cancer lesions are absolutely needed, although these are frequently hampered by the fact that at the time of cancer diagnosis or surgery, genetically related non-cancer clones are probably swept out by rapidly expanded cancer clones^8,12,13.

Breast cancer is one of the most prevalent cancers among women. Annually, 2.3 million women are diagnosed with breast cancer worldwide, with 685,000 deaths reported in 2020, and the incidence is still increasing in many countries^14,15. Whereas so far there have been no large studies of clonal expansion in normal breast tissues, several reports have shown that, in some patients, cancer is accompanied by satellite benign breast lesions (BBLs), such as proliferative lesions with and without atypia, wherein both lesions share common genetic alterations and hence a common ancestor^{16,17,18,19,20}. The phylogenetic analysis of such cases containing both cancer and associated BBLs could provide a unique opportunity to decipher the entire life history and evolutionary dynamics of breast cancer, which are instrumental in understanding when and where the breast cancer ancestor is born and evolves to acquire the cancer phenotype, and also to establish new strategies to predict or even prevent breast cancer development^21,22.

In the current study, we performed whole-genome sequencing (WGS) of multiple microscale samples obtained by laser-capture microdissection (LCM) from both cancer and clonally related BBLs, together with apparently normal lobules. Then, on the basis of the rate of mutation accumulation estimated from WGS of single-cell-derived organoids established from mammary epithelia, we reconstructed the phylogenetic trees including both cancer and non-cancer clones to infer the entire history of breast cancer.

To estimate the rate of mutation accumulation in normal mammary epithelial cells with ageing, we established 71 single-cell-derived organoids from epithelial cell adhesion molecule (EpCAM)-positive cells isolated from histologically normal mammary tissues in patients with breast cancer or from breast milk provided by healthy volunteers who were breastfeeding (Methods, Extended Data Fig. 1a, Supplementary Note 1 and Supplementary Fig. 1). Excluding six organoids for which polyclonal origins were suspected and another with a germline CDH1 mutation, we finally evaluated somatic mutations for 64 organoids established from six pre- (n = 20) and nine postmenopausal (n = 32) patients with breast cancer and 12 from six healthy women (Fig. 1a and Extended Data Figs. 1 and 2). The mutation distribution in 15 organoids showed a minor subpeak at a variant allele frequency (VAF) less than 0.25 corresponding to subclonal mutations acquired during cell culture, which were eliminated assuming a Gaussian mixture model (Methods). We identified a total of 58,385 single nucleotide variants (SNVs) and 3,955 small insertions and deletions (indels) as clonal somatic variants, including four PIK3CA mutations (Fig. 1a). When fit to the known Catalogue Of Somatic Mutations In Cancer (COSMIC) single base substitution (SBS) signatures, most SNVs were assigned to three clock-like signatures²³, SBS1 (9.9%), SBS5 (80.7%) and SBS40 (9.4%). SBS1 is characterized by the prominence of C>T transitions at CpG dinucleotides resulting from the spontaneous deamination of 5-methyl-cytosine²⁴, whereas SBS5 and SBS40 are ‘flat’ signatures of unknown aetiology^24,25, which are difficult to separate from each other and, hence, designated collectively as ‘SBS5/40’ in the subsequent analyses. According to the linear regression model, the number of SNVs significantly depended on age at sample collection, years after menopause, parity and the presence of a driver mutation. SNVs were accumulated at 19.5 mutations per genome per year before menopause, which was reduced to 8.1 mutations per genome per year after menopause, while the mutation number was reduced by 54.8 per delivery (Fig. 1b,c). The mutation rate was also affected by PIK3CA mutations, which increased the number of SNVs by 210.4, although this needs to be validated using additional PIK3CA-mutant clones, because the number of driver-mutated samples was still small (n = 4). The accumulation rate of indels was also reduced by 45%, from 1.3 mutations per genome per year before to 0.72 mutations per genome per year after menopause, while the mutation number was spared by 3.9 per delivery (Fig. 1d,e).

The evolutionary history of breast cancer was then investigated using phylogenetic analysis using those surgical specimens that contained both pathologically confirmed cancer and multifocal (three or more) non-cancer proliferative lesions (greater than or equal to 3 mm in diameter). After reviewing pathology reports of 156 patients with breast cancer, we found five patients for whom formalin-fixed paraffin-embedded (FFPE) specimens fulfilling the above-mentioned criteria had been preserved. We obtained multiple microscale samples (2.9 mm (0.5–10) in diameter on average) from both cancer and non-cancer lesions (13.8 per patient) using LCM (Extended Data Fig. 3). We analysed a total of 69 LCM samples for somatic mutations and copy-number alterations (CNAs) using WGS, based on which we reconstructed phylogenetic trees (Methods, Supplementary Notes 1 and 2 and Supplementary Figs. 2 and 4). These samples comprised histologically normal lobules (n = 6), non-proliferative (n = 1) and proliferative lesions (n = 33), classic-type lobular carcinoma in situ (LCIS) (n = 1), ductal carcinoma in situ (DCIS) (n = 20) and invasive ductal carcinoma (IDC) (n = 8). All patients were in their 40s, premenopause and carried Luminal A-like IDC (n = 3) (KU539, KU779 and KU957) or oestrogen receptor (ER) (+) human epidermal growth factor receptor 2 (HER2) (−) DCIS (n = 2) (KU582 and KU873) (Supplementary Table 1). As was the case with normal organoids, the mutations in non-cancer lesions were dominated by SBS1 and SBS5/40 mutations (Fig. 2a and Extended Data Figs. 4 and 5). A minority of lesions had varying contributions from APOBEC-related signatures (that is, SBS2 and SBS13, ref. ²³) in peripheral branches, which supports the role of APOBEC-induced mutagenesis during relatively late phases of carcinogenesis (Extended Data Fig. 4c).

In all five cases, the phylogenetic trees comprised just one or two large clades, in which frequently acquiring unique driver alterations of their own, multiple progenies derived from a single common ancestor evolved to give rise to both cancer and BBLs. These progenies expanded over years and by the time of cancer diagnosis, occupied a large area in the affected breast, replenishing both cancer and non-cancer tissues as well as normal ones (Fig. 2a and Extended Data Fig. 5). Conspicuously, in four of the five cases, the most recent common ancestors (MRCAs) in each clade harboured der(1;16) (Fig. 2 and Extended Data Fig. 5a–c). Moreover, in KU779 and KU539, there were two independent der(1;16)(+) clones having distinct breakpoints, suggesting a strong selective advantage conferred by der(1;16) (Fig. 2a,d and Extended Data Fig. 5a). der(1;16) is a recurrent abnormal chromosome highly characteristic of Luminal A breast cancer, particularly those with invasive lobular histology (ILC)²⁶ and generated by an unbalanced translocation that fuses chromosome 1q and 16p arms near centromere sequences, that is, der(1;16)(q10;p10), leading to 1q gain and 16q loss in common (Fig. 2b,c). The widespread distribution of der(1;16) clones was further investigated using fluorescence in situ hybridization (FISH) by detecting unbalanced chromosomal copy numbers of 1q and 16q arms. Showing varying histologies, including non-cancer proliferative lesions, LCIS, DCIS and IDC, the der(1;16)(+) clones were widely distributed within the mammary gland, spanning 62 mm (range 35–90 mm) regions on average (Fig. 2a and Extended Data Figs. 3c and 5a–c). der(1;16) signals were also detected in histologically normal lobules (Fig. 2a and Extended Data Figs. 3c and 5b).

In line with previous reports^9,10, all cancer clones found in the large DCIS and IDC lesions in KU779 (nos. 7a–e and 11a–e) were branched off from a long common trunk with short branches, which otherwise obscured the past history of these cancer clones (Fig. 2a). However, many normal and BBL samples collected at the same time enabled the analysis of the early history of these cancer clones, although the very recent history was still unclear. In addition, we were able to estimate the approximate timing of early branch points in the tree using the mutation rate measured for normal mammary epithelial cells. In particular, the timing of the acquisition of der(1;16) was more accurately pinpointed than that of other driver events, by maximizing the posterior probability of the observed numbers of duplicated and unduplicated mutations on 1q arm in der(1;16)(+) MRCA (Extended Data Fig. 6a–d and Methods). On an average, der(1;16) in six clones was estimated to be acquired at 10.6 (range 5.8−16.9) years of age (Fig. 2a and Extended Data Fig. 5a–c). We also estimated the average timing at which MRCA emerged as 26.5 (range 18.1–34.4) years of age, assuming a constant mutation rate until the emergence of the MRCA. For example, two distinct der(1;16) detected in a 48-year-old woman (KU779) were estimated to occur in two mammary cells at the ages of 5.8 and 10.0 years, respectively (Fig. 2a). These ancestor cells then gave rise to the MRCAs at the ages of 18.1 and 22.3 years, respectively, from which a number of non-cancer progenies evolved, followed by the appearance of cancer founders at least 10 years after the initial acquisition of der(1;16). In the remaining case (KU582) lacking der(1;16), the MRCA carrying an AKT1 mutation emerged by 4.4 years of age. It was 9.0 years later when the most recent common non-cancer ancestor of cancer clones was confirmed (Extended Data Fig. 5d). Targeted-capture sequencing of additional LCM samples for mutations in the main trunk confirmed that AKT1-mutated lobules derived from the MRCA widely expanded over the 65 mm area, showing a variety of histologies, including proliferative and non-proliferative as well as DCIS lesions (Extended Data Figs. 3c and, 5d). Of particular interest is the observation that multiple cancer founders (nos. 7a–e, 10, 14 and 15 in Fig. 2a, for example) independently evolved from non-cancer ancestors (blue arrows in Fig. 2a and Extended Data Fig. 5a,c,d) (KU779, KU539, KU957 and KU582). This suggests that a cancer population can be initiated at different time points by multiple independent cancer founders originating from common ‘non-cancer’ ancestors (Extended Data Fig. 6e).

Because these cancer and non-cancer lesions shared many driver alterations, in addition to der(1;16) and an AKT1 mutation, and had an identical germline background, WGS of these lesions provides a unique opportunity to investigate the critical genetic events that discriminate among IDC, DCIS, non-cancer proliferative lesions and histologically normal lobules. In some cases, additional driver alterations were observed only in high grade lesions but not in non-cancer clones in the most recent clade, whereas the opposite was true of other lobules (Fig. 2a and Extended Data Fig. 5). In KU873, for example, a DCIS clone (in no. 3) had two driver mutations in GATA3 and CBFB in addition to der(1;16) and a PTEN mutation, whereas no additional driver mutations were found in the remaining seven non-cancer lesions except for no. 8 carrying a GATA3 mutation (Extended Data Fig. 5b). Conversely, in KU779, a DCIS lobule (no. 14) had no known driver alterations other than der(1;16), whereas lobules in the same clade (nos. 3, 4, 8, 9 and 13) were still proliferative lesions without progression to DCIS or IDC, although they acquired an additional PIK3CA (E726K) with or without a RUNX1 (P267fs) and another PIK3CA (H1047R) mutation, or a PIK3CA (E545K) mutation (Fig. 2a). In KU539, the DCIS clone in no. 5 had two driver mutations in PIK3CA and CBFB in addition to a GATA3 mutation and der(1;16), although the invasive cancer clone (in no. 4) in the same clade had no additional driver alterations (Extended Data Fig. 5a). Similarly, no additional driver alterations were found in DCIS clones (in nos. 1 and 2) in KU957, whereas proliferative but non-cancer lesions (nos. 4 and 7) carried an additional copy of chromosomes 10 and 15q, respectively, with no. 4 also carrying an IDH1 mutation (Extended Data Fig. 5c). Overall, the most common additional mutations acquired during evolution were those affecting PIK3CA and GATA3 (Extended Data Fig. 5e), which showed no clear correlations with cancer clones, although they are among the most frequent mutations in breast cancer^27,28. Thus, no consistent patterns of additional driver mutations were observed between cancer and non-cancer clones.

The unexpected enrichment of der(1;16) in the five index cases suggested that the widespread expansion of satellite lesions of varying histology was a common feature of der(1;16)(+) breast cancer. To confirm this, we screened another set of 33 specimens of Luminal A-like invasive cancer (n = 28) or its putative precursor lesion (ER(+)HER2(−) DCIS) (n = 5) for der(1;16) using FISH and identified an additional eight der(1;16)(+) specimens, two from premenopausal and six from postmenopausal patients (Fig. 3a and Extended Data Figs. 7 and 8). As was the case with der(1;16)(+) clones in the index specimens (Fig. 2a and Extended Data Fig. 5a–c), which were all from premenopausal patients, the two der(1;16)(+) clones in premenopausal patients showed a macroscopic expansion over an area greater than 20 mm in diameter (Fig. 3b and Extended Data Fig. 7), supporting the above-mentioned hypothesis. By contrast, most of the remaining der(1;16)(+) clones from six postmenopausal patients were found in cancer lesions, rarely involving non-cancer lesions and, if ever, the surrounding der(1;16)(+) non-cancer lesions were confined within small lobules less than 10 mm in diameter (Fig. 3a,b and Extended Data Fig. 8a–e). To exclude the possibility that this was due to the late acquisition of der(1;16), we estimated the timing of der(1;16) acquisition in five of the six postmenopausal patients on the basis of phylogenetic analysis. Of interest, the mean age of the acquisition of der(1;16) in the five postmenopausal patients was estimated as 11.7 years (0–18.7), which is comparable to the 10.6 years (5.8–16.9) (P = 0.54) in premenopausal patients (Fig. 3c,d). Thus, we speculate that there should have been a larger expansion of der(1;16)(+) clones, including non-cancer lesions, before menopause, which, however, regressed after menopause in the face of reduced oestrogen levels.

At the onset of puberty, the mammary glands start to proliferate rapidly and undergo a remarkable expansion until the end of adolescence^29,30. Thus, to see whether the large expansion of der(1;16)(+) clones could be explained only by the physiological enlargement of the growing mammary glands in this period, we evaluated the extent to which der(1;16)(−) non-cancer clones can expand after puberty. For this purpose, we microdissected multiple non-cancer lobules of the surgically resected and fresh-frozen specimens from three newly recruited premenopausal patients with breast cancer (Supplementary Table 1 and Methods). On average, 2.1 (1–6) LCM samples were collected from each of the 10 × 10 mm² areas aligned over 40–50 mm lengths from three consecutive tissue slices separated from adjacent slices by 10 mm (Fig. 4a and Extended Data Fig. 9). In total, we collected 77 LCM samples, comprising 66 histologically normal lobules, eight proliferative and three classic-type LCIS lesions, which were subjected to WGS (Methods, Supplementary Note 1 and Supplementary Fig. 3). Six of the eight proliferative lobules harboured one or more breast cancer driver mutations (Fig. 4b,c). Driver mutations were also found in histologically normal lobules, although less common (12 of 66, P = 2.0 × 10⁻³), in which PIK3CA and PIK3R1 were significantly mutated (dN/dS > 1.0) (Fig. 4c and Supplementary Table 10). The presence of driver mutations was associated with a higher clonality as suggested by a significantly larger median VAF of mutations in driver-mutated versus unmutated samples (0.33 versus 0.25, P = 1.8 × 10⁻³), supporting the role of driver mutations in positive selection (Fig. 4d). One patient had a known pathogenic germline BRCA2 variant that, however, did not seem to influence the clonality in normal lobules (Extended Data Fig. 9a,d). Among 48 clones that were estimated to have emerged by 1 year of age, 13 (27.1%) involved two or more lobules, of which three clones affected the lobules more than 10 mm apart (Fig. 4e,f). By contrast, only five (9.3%) of 54 clones that emerged after 13 years affected two or more lobules, all of which were still confined within an area less than or equal to 10 mm in diameter. These observations indicate that most of the der(1;16)(−) clones that emerge after puberty stay within a single lobule or, if not, are confined to adjacent lobules and rarely expand to a larger area as observed for those carrying der(1;16). Thus, the large expansion of all der(1;16)(+) clones is not explained only by the physiological development of the breasts during puberty but suggests the driver role of der(1;16). The role of der(1;16) in clonal expansion is further highlighted by the unexpected detection of der(1;16) in three LCIS lesions in KU1215, which was shared with an ILC lesion located in another quadrant and shown to have expanded over the region spanning more than 70 mm in diameter in the subsequent FISH analysis (Extended Data Fig. 9b,c).

Finally, we characterized clinical and pathological features of der(1;16)(+) breast cancers and also investigated its role in breast cancer pathogenesis, using published data of 610 breast cancer cases from The Cancer Genome Atlas (TCGA) (Extended Data Fig. 10). According to the copy-number measurement based on exome sequencing², der(1;16) was detected in 119 (19.5%) cases as a concomitant whole-arm 1q gain and 16q loss (Extended Data Fig. 10b). In accordance with previous reports^26,31, most of the der(1;16)(+) cases (86.6%) were classified as Luminal A-type cancers (Extended Data Fig. 10c), in which der(1;16) was more enriched in the ILC tumours compared to that in the IDC tumours (49.6 versus 12.7%, P = 3.2 × 10⁻¹⁶) (Extended Data Fig. 10e). Whereas Luminal A tumours generally have a favourable prognosis, der(1;16)(+) tumours were associated with a significantly longer overall survival than der(1;16)(−) Luminal A tumours (a median of 33.1 months, compared with 28.3 months, P = 1.0 × 10⁻³) (Extended Data Fig. 10h). In multivariable analysis, der(1;16) remained to be a significant predictor of prolonged overall survival, even after adjustment for age and stage (P = 2.4 × 10⁻³, Supplementary Table 16). These pathological subtypes and mutation profiles were largely recapitulated in the der(1;16)(+) tumours that were analysed (Extended Data Fig. 8f). After excluding PIK3CA, which was frequently mutated in both der(1;16)(+) and der(1;16)(−) cases, the most frequent mutational targets in der(1;16)(+) Luminal A tumours included CDH1, and GATA3 and CBFB in ILC and IDC tumours, respectively. MAP2K4 and TP53 mutations were less common compared to those in der(1;16)(−) tumours (Extended Data Fig. 10i,j and Supplementary Table 17). CDH1 is a putative target of 16q deletion^26,32 and was mutated in 86 (14.1%) of the 610 TCGA cancer cases. Among these, most cases (94%) had biallelic alteration, in which most were associated with der(1;16) (n = 46), followed by other 16q loss (n = 21) and 16q UPD (n = 14). Heterozygous mutation was found in only five cases. Except for one frameshift change, all were variants of unknown significance that were rarely registered in the COSMIC database. These findings indicate that CDH1 is a bona fide recessive tumour suppressor gene, and haploid loss of CDH1 alone may not be sufficient for positive selection or clonal expansion during breast cancer development. Regarding this, we analysed the effect of the allelic imbalance caused by der(1;16) on gene expression in 323 Luminal A cancer cases, including 103 der(1;16)(+) cases from the TCGA. The mean expression levels of 1q and 16q genes were substantially increased or decreased in der(1;16)(+) cases compared to those in cases without 1q gain or 16q loss, respectively (Extended Data Fig. 10k). In particular, the mean expression levels of CDH1 in der(1;16)(+) cases without CDH1 mutations were almost halved compared with those in cases without 16q loss. However, the effect of haploid gain and loss was not confined to CDH1 but also seen in other known or putative oncogenes and tumour suppressor genes^27,28, although the effect was highly variable across genes. Thus, the exact molecular targets of der(1;16) are still elusive.

Through phylogenetic analyses, we successfully traced the evolution of breast cancer and precursor lesions, from the acquisition of initial driver alterations to the development of clinically diagnosed disease. The absolute timing and the order of early driver events were more accurately estimated than in previous studies^{9,10,11,33,34} by analysing both cancer and non-cancer lesions and by using the rate of mutation accumulation measured for normal mammary epithelium. As demonstrated in a recent study on myeloproliferative neoplasms³⁵, the first driver events occurred long before the cancer diagnosis, around puberty or late adolescence or, in one case, as early as in early infancy. However, unlike the case with the myeloproliferative neoplasms study, discrimination between cancer and non-cancer clones along the evolutionary tree was enabled to some time point after the acquisition of initial driver events. In most cases, the appearance of the MRCA of cancer and non-cancer clones was no earlier than the patient’s 20s to early 30s. Thus, it seems to still take more than 10 years from the acquisition of the initial driver alterations (at 6–17 years of age) before the initial cancer founders appear. Expanding along the mammary ducts, these clones ultimately occupied an unexpectedly large area in the breast by the time of cancer diagnosis. Of note, the cancer clones often evolved multifocally from clonally related but still ‘non-cancer’ ancestors. Being distinct from the classical linear model for the evolution of a single cancer founder, such a branching pattern of evolution of multiple cancer founders from within a non-cancer population might be more common than expected during cancer development (Extended Data Fig. 6e). Another finding of interest is the lack of consistent correlations between histologies and the number and/or type of driver events. Although we cannot exclude the possibility of the presence of undetected driver mutations and structural variations, this may suggest the role of epigenetic changes³⁶ and/or locally defined microenvironments in cancer development.

It should be noted that such a unique pattern of cancer evolution could be biased by the selection of specimens harbouring multiple satellite BBLs for LCM, which was highly enriched for der(1;16). The analysis of additional cases with der(1;16) confirmed that the presence of persistent non-cancer clones in a large area is an intrinsic feature of der(1;16)(+) breast cancer at least in premenopausal cases. The parallel evolution of multiple independent der(1;16) clones in two cases supports the strong driver role of der(1;16) in puberty or late adolescence. Accounting for 20% of all breast cancers and one- and two-thirds of Luminal A and invasive lobular breast cancers, respectively, der(1;16) defines a major subtype of breast cancers. However, it is still open to question whether or not this pattern of cancer evolution is also common in other breast cancer subtypes. It was observed at least in an AKT1-mutated case (KU582). Mutations affecting PIK3CA and PIK3R1 are among the most frequent targets of somatic mutations in breast cancer^27,28 and also common in apparently normal mammary lobules (10 out of 66 lobules) (Fig. 4a and Extended Data Fig. 9a,b). However, none of the clones carrying these mutations showed a widespread expansion. Further investigations are needed to clarify this.

Another finding of interest is the mutational profile of mammary epithelium, which is distinct from that in other tissues^1,37,38. Throughout a woman’s lifespan, the mammary epithelium undergoes dynamic changes that are synchronized with menstrual cycles as well as pregnancy, delivery and breastfeeding thereafter^{29,30,39,40,41,42}, which is also reflected by its unique mutation profile. A significantly reduced mutation rate after menopause might be explained by the reduced cell turnover associated with the cessation of menstrual cycles and/or reduced oestrogen levels. In agreement with this is the negative effect on accumulation of mutations imposed by pregnancy or delivery, during which menstrual cycles are spared. The findings on mutation profile are also in agreement with the epidemiological findings that late menopause and low parity correlate with an elevated risk of breast cancer development^43,44,45. However, the reduction of 50 SNVs per parity seems substantially larger than expected from the typical period (1.1–1.5 years) when menstrual cycles are interrupted by pregnancy and breastfeeding (roughly 30 SNVs)⁴⁶. This raises the possibility that after the effacement of markedly proliferated mammary glands after delivery or breastfeeding, the mammary epithelium might be reconstructed by newly recruited, ‘dormant’ stem cells, in which the SNV burden had been spared. Such dormant stem cells have been proposed to explain the disappearance of clones carrying tobacco signatures in bronchial epithelium after cessation of tobacco smoking⁴⁷.

In summary, we revealed mutational processes in the mammary epithelium and the entire life history of breast cancer, highlighting a unique role of der(1;16) in the major subset of Luminal A breast cancer. Our findings may contribute to the understanding of breast carcinogenesis and the development of new strategies for prediction, early diagnosis and even prevention of breast cancer.

Data reporting

No statistical methods were used to determine the sample size. The experiments were not randomized. Pathologists were blinded to the genetic alterations in each sample during histopathological evaluation.

Participants and materials

We enroled 207 female patients with breast cancer who underwent surgery at the Kyoto University Hospital and eight healthy breastfeeding women who delivered at the Kyoto University Hospital or Adachi Hospital. Written informed consent was obtained from all the participants. The study was reviewed and approved by the ethics committees of the Kyoto University and Adachi Hospital, and all the participants consented to publication of data. The characteristics of the participants are summarized in Supplementary Table 1. Invasive cancer lesions in FFPE surgical specimens were immunostained for ER (no dilution), PR (no dilution), HER2 (no dilution) and Ki-67 (1:100 dilution), for which histological grade was also evaluated according to the modified Scarff–Bloom–Richardson grading system⁴⁸, to surrogate cancer subtypes. The cut-off for ER and PR-positivity was set at greater than or equal to 1%. For HER2 status, immunohistochemistry scores of 0 and 1+ were considered negative, whereas 3+ was considered positive. Tumours with scores of 2+ were further evaluated by means of dual colour in situ hybridization, wherein an HER2/CEP17 ratio greater than or equal to 2.2 was considered positive. The Ki-67 labelling index was determined in the hotspots. Surrogate subtype classification was defined as follows: Luminal A-like, ER(+) and HER2(−) and histological grade 1 or 2 and Ki-67 less than or equal to 15%; Luminal B-like, ER(+) and HER2(−) and histological grade 3 or Ki-67 greater than 15%; Luminal HER2, ER(+) and HER2(+); HER2-enriched, ER(−) and HER2(+); triple-negative breast cancer (TNBC), ER(−) and HER2(−). Only the nuclear grade was evaluated for DCIS classification according to the WHO Classification of Tumours of the Breast⁴⁹; ER, PR, HER2 and Ki-67 were not routinely evaluated. BBLs were also evaluated according to the World Health Organization (WHO) classification and classified as follows: non-proliferative lesions (fibroadenoma without atypia, columnar cell change and apocrine metaplasia), proliferative lesions without atypia (usual ductal hyperplasia, columnar cell hyperplasia, sclerosing adenosis, radial scar and papilloma) and proliferative lesions with atypia (flat epithelial atypia, atypical ductal hyperplasia and atypical lobular hyperplasia). Classic-type LCIS was also classified as BBL per the recent clinical practice wherein LCIS and atypical lobular hyperplasia were grouped as ‘lobular neoplasia’, which is considered to be a risk factor for cancer development but is not an obligate precursor of invasive cancer^22,50. When a lesion consisted of two or more differentially classified epithelia, the most severe diagnostic category was assigned. Breast cancers were classified by referring to the medical records, whereas BBLs and cancer lesions subjected to genetic analysis were independently reviewed by three experienced pathologists blinded to the genetic alterations. In case a unanimous agreement was not reached for the lesions, the most experienced pathologist reviewed them again and determined the consensus diagnosis.

To establish normal epithelial cell-derived organoids, fresh normal breast tissue and matched blood samples were obtained from 15 patients with breast cancer who underwent total mastectomy, and excess breast milk and oral mucosa were obtained from eight healthy breastfeeding women. For the analysis of cancer-related clonal evolution, 156 patients with cancer who underwent surgery without any preoperative treatment were screened by searching terms related to proliferative lesions in the pathological reports, based on which five archival FFPE surgical specimens and matched blood DNA or FFPE normal lymph nodes were provided by the Kyoto Breast Cancer Research Network (KBCRN) Breast Oncology Research Network (BORN)-BioBank. We also obtained an additional 33 FFPE surgical specimens and five matched FFPE normal lymph node or skin samples from the BORN-Biobank, for which tumour cores in the tissue microarray had already been evaluated for surrogate subtype classification. To investigate the structures of der(1;16)(−) non-cancer clones, fresh-frozen tissues and matched blood samples were obtained from three breast cancer patients who had undergone total mastectomy; FFPE surgical specimens were also used for one of three patients. The sample information is summarized in Supplementary Table 2.

External datasets

Whole-exome sequencing (WES) data (.bam files) of paired tumour and germline control samples from female patients with invasive breast cancer (n = 661) were downloaded from TCGA data portal (https://portal.gdc.cancer.gov/). The .bam files were converted to the fastq format using biobambam⁵¹ (v.0.0.191) and processed using the Genomon2 pipeline (v.2.6) for mutation calling (below) and our in-house pipeline ‘CNACS’ for copy-number analysis as described by Yokoyama et al.². RNA-sequencing (RNA-seq) data in the transcripts per million (TPM) format were also downloaded from the TCGA data portal. Clinicopathological information of these samples was also downloaded from TCGA data portal, whereas the information about PAM50 messenger RNA subtypes was extracted from the study by Ciriello et al.³²; if data were lacking, information was extracted from TCGA Network²⁷. Samples included in this study are summarized in Supplementary Table 2.

Organoid culture

Single normal epithelial cell-derived organoids were established according to the protocols described by Lim et al.⁵², Wong et al.⁵³ and Dekkers et al.⁵⁴, with some modifications. First, single-cell suspensions were obtained from normal breast tissues of patients with breast cancer or the breast milk of healthy breastfeeding women. Fresh mammary tissue from the contralateral quadrant of the cancer-containing one was obtained from surgical specimens, which were confirmed to be pathologically normal by three pathologists reviewing the haematoxylin and eosin (HE)-stained fresh-frozen and FFPE sections. The tissue was minced manually and digested for 8–10 h at 37 °C with 150 U ml⁻¹ collagenase I (Thermo Fisher Scientific (Thermo)), 50 U ml⁻¹ hyaluronidase (Merck) and 100 U ml⁻¹ DNase I (Roche) in Dulbecco’s modified Eagle’s medium (DMEM)/F-12 supplemented with 5% fetal bovine serum (FBS), 0.5 mM glutamine (Thermo), 5 μg ml⁻¹ insulin (Merck), 10 ng ml⁻¹ epidermal growth factor (EGF) (PeproTech) and 500 ng ml⁻¹ hydrocortisone (Merck)⁵². The resulting cell suspension was sequentially digested with 0.25% trypsin and 1 mM EDTA (1 min, 37 °C), and 5 mg ml⁻¹ dispase (Thermo; 1 min, 37 °C). A single-cell suspension was obtained by means of filtration through a 40 μm cell strainer (Corning) after removing red blood cells using RBC Lysis Solution (QIAGEN). Excess breast milk was stored at 4 °C and transported to the laboratory within 24 h. Milk was diluted 1:1 with PBS and centrifuged at 1,000g for 10 min (ref. ⁵³). The supernatant, including milk fat, was discarded and the cell pellet was washed 3–5 times with PBS; a single-cell suspension was then obtained by means of filtration through a 40 μm cell strainer.

Next, mammary epithelial cells were isolated from the single-cell suspension using CD326 EpCAM MicroBeads (Miltenyi Biotec (Miltenyi)) (1:5 dilution) and the MACS Cell Separation System (Miltenyi) according to the manufacturer’s instructions. For the cell suspension obtained from the breast milk of healthy participants without available paired oral mucosa, the negative selection was performed using CD45 MicroBeads (Miltenyi) (1:5 dilution) before EpCAM-positive epithelial cell isolation to use CD45-positive leukocytes as paired control samples for WGS.

The isolated epithelial cells were cultured using two different methods. In culture method 1, the cells were resuspended in Matrigel (Corning) at a density of 5,000–20,000 cells per ml for tissue-derived cells and 20,000–80,000 cells per ml for milk-derived cells; 100 μl suspension was seeded into each well of six-well plates, which was filled with 2 ml DMEM/F-12 supplemented with 1× B27 supplement (Thermo), 0.5 mM glutamine, 5 μg ml⁻¹ insulin, 10 ng ml⁻¹ EGF, 500 ng ml⁻¹ hydrocortisone, 20 ng ml⁻¹ cholera toxin (Merck), 100 U ml⁻¹ penicillin and 100 μg ml⁻¹ streptomycin⁵². The cells were cultured for 13–37 days (median 17 days) until the organoids grew to equal to or more than 80 μm in diameter and were collected. In culture method 2, the epithelial cells that had been stocked frozen in CELLBANKER 1 (Takara Bio (Takara)) were suspended in 25 μl Matrigel and plated on 48-well plates, wherein each well was filled with the culture medium based on the previous protocols⁵⁴ with modifications: 250 μl DMEM/F-12 supplemented with 10 mM HEPES, 1× B27 supplement, 2 mM glutamine, 10 nM Gastrin I (Merck), 1 mM N-acetylcysteine (FUJIFILM Wako Pure Chemical (FUJIFILM)), 1 μg ml⁻¹ R-Spondin 1 (R&D Systems), 5 nM Neuregulin 1 (PeproTech), 20 ng ml⁻¹ FGF-10 (PeproTech), 5 ng ml⁻¹ EGF, 100 ng ml⁻¹ Noggin (PeproTech), 500 nM A83-01 (Tocris Bioscience), 5 μM Y-27632 (FUJIFILM), 0.5 μg ml⁻¹ hydrocortisone (Selleck Chemicals), 100 nM β-Oestradiol (Cayman Chemical Company), 10 μM Forskolin (Merck), 100 U ml⁻¹ penicillin, 100 μg ml⁻¹ streptomycin and 10% Afamin-Wnt-3A serum-free conditioned medium⁵⁵. After 15–82 days (median 20 days) of primary culture, single-cell-derived organoids were established by seeding dissociated cells at a limiting dilution, treated with TrypLE Express (Thermo) and then cultured for 13–88 days (median 15 days). In total, we established 71 single-cell-derived organoids from 15 patients with breast cancer and eight healthy participants.

LCM

FFPE surgical specimens from patients with breast cancer were prepared for LCM using the protocol described by Uehiro et al.⁵⁶. Individual 10 μm-thick FFPE specimens were placed on PEN membrane (4 µm)-coated glass slides (Leica Microsystems (Leica)) and immunohistochemically stained with pan-cytokeratin antibody cocktails (AE1/AE3) no. 412811 (NICHIREI) (no dilution). HistoZyme (Diagnostic BioSystems) was used for antigen retrieval, and the VECTOR Red Alkaline Phosphatase Substrate Kit (VECTOR Laboratories) was used for visualization. For the analysis of der(1;16)(−) non-cancer lobules, surgical specimens of premenopausal patients with breast cancer who underwent total mastectomy were sliced into 10-mm-wide slices and 10 × 10 mm² tissues were consecutively obtained from the slices of the contralateral quadrant of the cancer-containing one and immediately frozen. Frozen tissues were sectioned using a cryostat CM1950 (Leica). Ten micrometre-thick sections were placed on PEN membrane-coated slides and stained with Mayer’s Hematoxylin solution (FUJIFILM) and Eosin Y (FUJIFILM). LCM of the stained FFPE and frozen tissue slides was performed using the LMD7000 or LMD7 system (Leica). The pathologists reviewed the HE-stained slides, and CK5- and E-cadherin-stained slides (1:100 and 1:50 dilutions, respectively) if needed, for each of the 10–15 sequentially sectioned 10-μm-thick LCM slides and diagnosed each dissected lesion.

WGS

The DNA extracted from each single organoid established in culture method 1 was divided into two aliquots, each of which was independently subjected to whole-genome amplification (WGA) with the REPLI-g Single Cell Kit (QIAGEN) and analysed by means of WGS and subsequent validation sequencing³. DNA from the fresh organoids successfully expanded in culture method 2, peripheral blood, oral mucosa and leukocytes derived from breast milk was extracted using the QuickGene DNA whole blood kit (Kurabo Industries), Gentra Puregene Kit (QIAGEN), QIAamp DNA Blood Mini Kit (QIAGEN) or QIAamp DNA Micro Kit (QIAGEN). DNA from FFPE and fresh-frozen LCM samples was extracted using the GeneRead DNA FFPE Kit (QIAGEN) and Maxwell 16 Cell LEV DNA Purification Kit (Promega), respectively.

WGS libraries were prepared as follows: 100 ng of the WGA DNA extracted from single organoids in culture method 1 was used to prepare a library using the TruSeq Nano DNA Library Prep Kit (Illumina) or Lotus DNA Library Prep Kit (Integrated DNA Technologies (IDT)); 5–50 ng of DNA extracted from fresh organoids in culture method 2 was used to prepare a library using the xGen Prism DNA Library Prep Kit (IDT); 10–200 ng of DNA extracted from FFPE LCM samples was used to prepare a library using the SMARTer ThruPLEX DNA-seq Kit (Takara) or xGen Prism DNA Library Prep Kit; 2.5–30 ng of DNA extracted from fresh-frozen LCM samples was used to prepare a library using the Lotus DNA Library Prep Kit. These libraries were sequenced on a NovaSeq 6000 system (Illumina) or DNBSEQ-G400RS (MGI Tech) in 100–150-basepair (bp) paired-end mode, according to the manufacturer’s instructions. In total, paired WGA samples from 65 organoids, six fresh organoid samples, 84 FFPE LCM samples, 79 fresh-frozen LCM samples and 36 matched germline controls were used. The target coverage was 35× for FFPE samples and 30× for other types of sample, and the actual average coverage was as follows: 35× (19–59×) in WGA organoid samples, 45× (40–61×) in fresh organoid samples, 46× (24–100×) in FFPE LCM samples, 34× (18–63×) in frozen LCM samples and 41× (28–73×) in germline samples.

Raw sequence data were processed into .bam files using Genomon2, as previously described². In brief, sequencing reads were aligned to the human reference genome (GRCh37) using the Burrows–Wheeler Aligner⁵⁷ (v.0.7.8) with default parameter settings. The PCR duplicates were eliminated using biobambam. For organoid samples, mouse-derived sequencing reads resulting from Matrigel contamination were removed using Xenome⁵⁸ (v.1.0.0), and only the reads classified as ‘human-mapped’ were processed into .bam files. Mutation calling for WGA organoid samples was performed after merging each pair of two .bam files of WGA samples derived from a single organoid using Samtools⁵⁹ (v.1.10); each mutation call was then reviewed back to the original two .bam files using GenomonMutationFilter² (v.0.2.1); mutations detected in both the WGA samples with two or more variant reads each were considered true somatic mutations, whereas mutations detected in only one sample were eliminated as WGA-related errors. The mutations listed in our in-house mouse-derived variant list were eliminated as artefacts due to Matrigel contamination.

In the entire WGS analysis, mutation calling was performed by means of paired analysis using a ‘three-caller combination’ to improve the sensitivity and true positive rate, wherein mutations were called by three different callers (Genomon2, Mutect2 (ref. ⁶⁰) (GATK4, ref. ⁶¹ (v.4.1.2)), and Strelka2 (ref. ⁶²) (v.2.9.3)) independently; the mutations detected by two or three callers were considered ‘high-confidence’ mutations. As described previously², Genomon2 first discards the low-quality, unreliable reads and variants with mapping quality of less than 20 and/or base call quality of 15 or lower. Next, the variants that did not meet the following criteria were further excluded as sequencing errors: (1) a sufficient depth (six or more) in both samples and the matched controls; (2) VAFs greater than 0.1, 0.15 and 0.12 in WGA and/or fresh organoid, FFPE LCM and fresh-frozen LCM samples, respectively; (3) variant reads three, five and four or more in WGA and/or fresh organoid, FFPE LCM and fresh-frozen LCM samples, respectively; (4) VAFs less than 0.07 and variant reads one or less in matched controls; (5) a strand ratio not equal to 0 or 1; (6) Fisher’s P < 10^−1.5 for WGA and/or fresh organoid samples and <10^−1.3 for FFPE and fresh-frozen LCM samples; (7) EBCall⁶³ P < 10⁻⁵, <10⁻⁴, <10⁻⁵ and <10⁻³ for WGA organoid, fresh organoid, FFPE LCM and fresh-frozen LCM samples, respectively, which were evaluated with a ‘control panel’ consisting of WGS data of 39, 20, 42 and 14 blood or normal tissue samples of unrelated participants prepared using the corresponding library preparation kits; for the analysis of WGA organoid and FFPE LCM samples, which were expected to contain a lot of sample type-specific artefacts, a ‘control panel’ consisting of an increasing number of normal samples of corresponding sample type was used to eliminate artefacts. The .bam files were further edited by means of Samtools to remove sequencing reads with mapping quality below 20 and duplicated reads and then analysed using Mutect2 and Strelka2. Mutation calling using Mutect2 was performed as follows: initially, variants were called by Mutect2 using panel of normals made from ‘control panel’ data, to filter out sequencing noise; the raw output of Mutect2 was subsequently processed by means of FilterMutectCalls and FilterByOrientationBias in the default settings to filter out the remaining sequencing errors. We excluded the variants that were not supported by (1) a sufficient read depth (eight or more) in both samples and the matched controls; (2) VAFs greater than 0.1, 0.15 and 0.12 in WGA and/or fresh organoid, FFPE LCM, and fresh-frozen LCM samples, respectively; (3) allelic depths for the alternative alleles four or more in WGA and/or fresh organoid and fresh-frozen LCM samples and five or more in FFPE LCM samples; (4) VAFs less than 0.07 and allelic depths for the alternative alleles two or less in matched controls; (5) Phred-scaled quality for the possibility of sequencing errors (SEQQ) greater than 40 (for SNVs); (6) Phred-scaled quality of strand bias artefact (STRANDQ) greater than 70 for WGA and/or fresh organoid samples and greater than 50 for FFPE and fresh-frozen LCM samples (for SNVs). Mutation calling using Strelka2 was performed in the default settings, and the variants that were not supported by (1) a sufficient read depth (ten or more) in both samples and the matched controls; (2) VAFs greater than 0.1, 0.15 and 0.12 in WGA and/or fresh organoid, FFPE LCM and fresh-frozen LCM samples, respectively; (3) the alternative alleles five or more in WGA and/or fresh organoid and fresh-frozen LCM samples and six or more in FFPE LCM samples; (4) VAFs less than 0.07 and the alternative alleles two or less in matched controls; (5) a somatic Empirical Variant Score (SomaticEVS) greater than 17 for SNVs, greater than 16 for indels in WGA organoid samples and greater than 6 for indels in fresh organoid, FFPE LCM and fresh-frozen LCM samples, were excluded. The variants identified by each caller were annotated using ANNOVAR⁶⁴; the variants that were listed in the 1000 Genomes Project dataset or gnomAD database with a minor allele frequency of more than or equal to 0.001 and variants within segmental duplications reported in the GenomicSuperDups database or repetitive sequences reported in the University of California, Santa Cruz (UCSC) Genome Browser⁶⁵ were further excluded, except for the driver mutations defined below, to achieve a high true positive rate. After that, each variant was reviewed in .bam files of each sample and the matched control as well as ‘control panel’ samples using GenomonMutationFilter: variants that were not supported by (1) VAFs greater than 0.1, 0.15 and 0.12 in WGA and/or fresh organoid, FFPE LCM and fresh-frozen LCM samples, respectively; (2) the variant reads of three, four and five or more in fresh organoid samples, WGA organoid and fresh-frozen LCM samples and FFPE LCM samples, respectively; (3) variant read one or less in matched controls; (4) average VAFs less than or equal to 0.05 and total variant read of ten or less in ‘control panel’ samples; (5) VAFs in samples 15 times or more higher than average VAFs in ‘control panel’ samples, were further excluded as sequencing artefacts. Variants clustered within a short length (150 bp) were subjected to visual inspection using the Integrative Genomics Viewer⁶⁶ to further eliminate sequencing errors. Somatic driver mutations were defined as loss-of-function mutations in tumour suppressor genes or mutations reported in the COSMIC database with ten or more mutated tumours in all types of cancer or five or more mutated tumours in breast cancers.

For the evaluation of germline variants, germline samples were analysed using Genomon2 in the unpaired mode. The variants in breast cancer susceptibility genes that were registered as ‘pathogenic’ or ‘likely-pathogenic’ in the ClinVar database were considered pathogenic variants.

Copy-number analysis

CNAs were analysed using Control-FREEC⁶⁷ (v.11.0) with the contaminationAdjustment option, except for the WGA organoid samples in which detection of small copy-number changes was difficult due to WGA-related artefacts. The copy-number gains and losses, and the uniparental disomies were visually confirmed using normalized copy numbers and beta allele frequency plots. CNA results are summarized in Supplementary Table 3.

Estimation of mutation rate in normal epithelial cells

The mutation accumulation rate was estimated using clonal mutations detected in each single-cell-derived organoid, which was thought to have been inherited from the original single cell. We excluded six organoids in which polyclonal origins were suspected (median VAF less than 0.4); furthermore, we excluded an organoid from a participant carrying a CDH1 germline pathogenic variant from the analysis. Clonal mutations were defined by the following method: (1) the Gaussian mixture model was adapted to VAF distribution of SNVs in each organoid (n = 64) using the R package mclust⁶⁸ (v.5.4.7) to determine the optimal number of mixture components, one or two, by which the organoids consisting of two components with VAF peaks greater than or equal to 0.4 and less than 0.25 were classified as ‘bimodal’, and those consisting of one or two components with VAF peaks greater than or equal to 0.4 were classified as ‘unimodal’ and (2) all the mutations in ‘unimodal’ organoids (n = 49) were considered clonal mutations; by contrast, mutations with VAF values lower than the intersection point of two components in ‘bimodal’ organoids (n = 15) were eliminated as subclonal mutations acquired during cell culture.

A linear regression model without intercept was fitted to estimate the effects of age and known breast cancer risk factors on the number of clonal SNVs and indels in 61 organoids, wherein the reciprocal of the number of analysed organoids per case was weighted for each organoid. Initially, we tested the effects of years before menopause (age at sample collection in premenopausal women and age at menopause in postmenopausal women), years after menopause (zero in premenopausal women and the difference between age at sample collection and age at menopause in postmenopausal women), the presence of driver mutation, presence of breast cancer, known breast cancer risk factors (body mass index, parity, alcohol consumption per day, smoking pack-years and number of first- and second-degree relatives with breast cancers) and the factors that might affect the sensitivity of mutation call (sample type (WGA or fresh organoids) and the fraction of genomic region with more than or equal to 15× coverage) on the number of mutations simultaneously. Then, we removed the least significant variable one by one by comparing models with and without each variable until the removal of any variable led to a significant difference (P < 0.05) by analysis of variance. Finally, we defined the final model with the least significant coefficients. Years before and after menopause were both extracted as significant coefficients on mutation number in the final models for SNVs and indels, which showed significant improvement compared to the models that did not incorporate the age at menopause (by analysis of variance). All statistical analyses are summarized in Supplementary Table 4.

Phylogenetic analysis

In the multi-sampling analysis, phylogenetic trees were reconstructed using somatic mutation data, copy-number information and mutant cell fraction (MCF) across all LCM samples in each subject. First, the Gaussian mixture model was adapted to VAF distribution of mutations within diploid regions to decompose major and minor clones. Then, MCF was calculated by doubling the mean value of the largest Gaussian distribution.

Before reconstructing a phylogenetic tree using MEGA⁶⁹ (v.11.0.11), we made an input matrix comprising all mutations detected across all samples, combined with their mutation status for all samples, which was determined according to the depth, number of supportive reads and copy-number status, as follows: (1) for samples with two or more supportive reads, the mutation status was assigned as ‘mutation_present’; (2) for samples with only one supportive read, the mutation status was assigned as ‘mutation_unknown’; (3) for samples with no supportive reads, (a) when chromosomal loss or other loss of heterozygosity was present at the mutation locus, the mutation status was assigned as ‘mutation_unknown’; (b) when no supportive read for the mutation was well expected (P > 0.05) according to the binomial distribution determined by sequencing depth, total copy number (TCN) and MCF, the mutation status was also assigned as ‘mutation_unknown’ and (c) otherwise, the mutation status was assigned as ‘mutation_absent’. On the basis of these criteria, we confirmed high accuracy for mutation status assignment as ‘mutation_present’ or ‘mutation_absent’ by validation sequencing (below) for 780 randomly selected mutations in two cases, wherein mutation status in two and nine samples, respectively, was evaluated for each mutation (accuracy 99.4%, 3,055 out of 3,072 mutation statuses); the results of validation sequencing are summarized in Supplementary Table 5. Then, a maximum parsimony tree was established using MEGA with 1,000 bootstrap replicates.

Subsequently, we used the R package treemut³⁵ (v.1.1) to assign each mutation to a branch using the expectation maximization method based on the number of supportive reads and the sequencing depth for all potential mutations for all samples, as well as the tree information from MEGA (branching pattern and branch length). However, because treemut was originally developed for the analysis of monoclonal diploid samples and assumes that VAF = 0.5 for mutated loci and VAF = 0 for wild-type loci, we corrected variant read counts on the basis of adjusted VAF (aVAF), MCF, TCN and minor copy number (MCN) for bulk LCM samples, as if they consisted of a clonal population derived from a single cell. At first, mutant allele number (MAN) in a single mutant cell was calculated assuming an MCF of 100%, using the following formula:

Next, adjusted TCN (aTCN) was determined as follows:

if MAN ≤ 0.5, aTCN = 2,

if MAN > (TCN − MCN), aTCN = (TCN − MCN) × 2,

otherwise, aTCN = (the nearest whole number from MAN) × 2.

Then, adjusted VAF (aVAF) was calculated by dividing MAN by aTCN.

Finally, corrected variant read counts were calculated by multiplying depth by aVAF.

The resulting phylogenetic trees were drawn with the SNV number to infer the time course of clonal evolution by adapting the SNV accumulation rate estimated in single-cell-derived organoids. Cancer driver mutations and CNAs were annotated manually in the trees as follows: in the cases in which more than one sample in a clade carried the same hot spot mutations with VAFs greater than 0.15 and variant reads of five or more for FFPE LCM samples, or VAFs greater than 0.12 and variant reads of four or more for fresh-frozen LCM samples, the mutation was assigned to the shared branch if it was accompanied by other shared mutations in that clade or assigned to the private branches if it was the only mutation shared by the samples in that clade; in the cases in which more than one sample in a clade shared CNAs, losses or gains of the same paternal or maternal alleles as determined by single-nucleotide polymorphism analysis were considered the same events and assigned to the shared branch, whereas those of the different paternal or maternal alleles were considered different events and assigned to the private branches. Information about whether each CNA was assigned to a shared or private branch is summarized in Supplementary Table 6.

The method of phylogenetic analysis using MEGA and treemut was validated by reconstructing the representative two trees (KU539 and KU779) using PyClone-VI (ref. ⁷⁰) (v.0.1.0), wherein no major inconsistency was observed between the two methods in terms of the topology and assignment of driver mutations to corresponding branches (Supplementary Note 2).

Analysis of mutational signatures

Mutational signatures were extracted using the R package MutationalPatterns⁷¹ (v.3.4.0) in strict signature refitting mode. In the multi-sampling analysis, signature extraction was performed using SNVs on branches with 100 or more SNVs in the phylogenetic trees, wherein each branch was treated as an individual sample. Initially, SNVs were allocated to all known COSMIC SBS signatures (v.3.1). The overall contribution of each COSMIC SBS signature was calculated in organoid, FFPE LCM and fresh-frozen LCM samples. Then, SBS signatures with a 5% or more contribution in one or more sample type were selected to prevent overfitting. Finally, SNVs were re-allocated to five SBS signatures (SBS1, SBS2, SBS5, SBS13 and SBS40). Because two clock-like ‘flat’ signatures, SBS5 and SBS40, are difficult to separate^24,25, they were collectively designated as ‘SBS5/40’ in the subsequent analyses. To validate the signatures extracted using MutationalPatterns, we also analysed SBS signatures using two more algorithms, the SigProfiler Bioinformatic Tools MatrixGenerator⁷² (v.1.1.27) and Extractor⁷³ (v.1.1.1), and the R package HDP⁷⁴ (v.0.1.5). In the SigProfiler analysis, de novo signatures were extracted using 1,000 nonnegative matrix factorization iterations, followed by further decomposition to COSMIC SBS signatures (SBS1, SBS2, SBS5, SBS7a, SBS7b, SBS13 and SBS18). In the HDP analysis, SBS signatures were analysed using COSMIC SBS signatures as priors, wherein six prior SBS signatures (SBS1, SBS2, SBS13, SBS16, SBS18 and SBS45) and two new signatures were extracted. The new signatures were further deconvoluted to COSMIC SBS signatures using the R package deconstructSigs⁷⁵ (v.1.8.0), resulting in SBS1, SBS5 and SBS7a. Because SBS7a, SBS7b, SBS16, SBS18 and SBS45 were not extracted by MutationalPatterns and were inconsistent between SigProfiler and HDP in many parts, we considered these signatures with less confidence (Extended Data Fig. 4a). SBS1, SBS2, SBS5/40 and SBS13 were extracted by means of all three extraction algorithms, and the contributions in each sample showed a similar pattern (Extended Data Fig. 4b).

Estimation of the timing of MRCA emergence and der(1;16) acquisition

The MRCA in the clonally related samples was identified on the basis of the phylogenetic tree. For common MRCAs of cancer and non-cancer clones (Fig. 2 and Extended Data Fig. 5), we assumed that the mutation accumulation rate was constant until the timing of emergence of the MRCA (T_MRCA) and equal to that of normal cells. A point estimate of T_MRCA was calculated by dividing the number of SNVs in MRCA (N_MRCA) by the SNV accumulation rate before menopause in normal cells (R₀), and the 95% CI (confidence interval) was calculated using the R₀ values randomly sampled 1,000 times on the basis of a normal distribution estimated from the linear regression model in the single-cell-derived organoids.

For the timing of der(1;16) acquisition, that is, the timing of 1q gain (T_{1q_gain}), we first obtained the number of duplicated mutations in MRCA using the MAN for each SNV on 1q, which was calculated as described above; if the average number of MAN of the related samples was 1.5 or less, such a mutation was considered ‘unduplicated’ or otherwise considered ‘duplicated’; then, the number of ‘unduplicated’ and ‘duplicated’ SNVs in MRCA was counted as N1 and N2, respectively. T_{1q_gain} was estimated by maximizing the posterior probability of the observed number of N2 based on the simulation method wherein the value of T_{1q_gain} was increased by 0.1 years from 0 to T_MRCA or age at sampling and N2 was simulated 1,000,000 times using the R₀ value that was randomly sampled 1,000 times as in T_MRCA estimation, assuming that mutations occurred on each 1q allele at a constant rate (R_{0_1q} = R₀ × (the size of 1q gained)/(the size of the genome) × 0.5 (for haploid)) according to the Poisson distribution. For the clones with two 1q gains (KU957, TMA114 and TMA125 (Extended Data Figs. 5c and 8a,b)), the timing of der(1;16) genome doubling (T_gd) was also simulated, wherein T_{1q_gain} and T_gd were increased by 0.1 years from 0 to T_gd, and from T_{1q_gain} to T_MRCA or age at sampling, respectively. The SNVs with the average MAN of more than 2.5 were considered ‘triplicated’ (N3). The values of T_{1q_gain} and T_gd with the maximum probability for the observed N2 and N3 were considered to be point estimates of T_{1q_gain} and T_gd, respectively. The details of the simulation are shown in Extended Data Fig. 6a–d, and the results are summarized in Supplementary Table 7.

FISH analysis

Unstained, 4 μm-thick FFPE sections were subjected to hybridization with bacterial artificial chromosome clone-derived probes for 1p.31.3, 1q23.3 and 16q23.2 to detect all types of der(1;16)(+) clone with concurrent whole-arm 1q gain and 16q loss, wherein the 1p.31.3 signals were used as controls. For KU779, in which two independent der(1;16)(+) MRCAs were identified with WGS, FFPE sections were also subjected to hybridization with Vysis DNA probes for CEP6 (D6Z1) and CEP16 (D16Z3) (Abbott Laboratories) to distinguish the two clones by D16Z3 signal counts, wherein the D6Z1 signals were used as controls (Fig. 2a,d). The probes used are listed in Supplementary Table 8. The hybridized slides were stained with 4,6-diamidino-2-phenylindole and examined under a BZ-X800 fluorescence microscope (KEYENCE). The FISH-probed signals were counted in 100 or more nuclei or in three or more microscopic fields per lesion.

Measurement of the clonal expansion areas

In the der(1;16)(+) breast cancer cases, lesions originating from der(1;16)(+) MRCAs were identified by performing FISH on FFPE surgical specimens. For five of 141 der(1;16)(+) lesions identified by means of FISH, we confirmed that they actually carried both der(1;16) and more than 90% of shared mutations in the MRCAs by means of targeted-capture sequencing using Next-Generation Sequencing (NGS) Discovery Pools (IDT) and the xGen CNV Backbone Hyb Panel (IDT) (below) (Supplementary Table 9). In the case KU539 with two independent der(1;16)(+) MRCAs (Extended Data Fig. 5a), two clones could not be distinguished using FISH, unlike that in the case KU779 (Fig. 2a,d), but could only be distinguished by evaluating the shared mutations in the MRCAs by means of sequencing, wherein the der(1;16)(+) lesions evaluated by means of FISH were considered as carrying ‘undetermined’ type of der(1;16) and were eliminated from the target lesions for measurement of the clonal expansion area. In the AKT1-mutated breast cancer case (Extended Data Fig. 5d), lesions originating from an AKT1-mutation(+) MRCA were identified by means of targeted-capture sequencing (below), wherein lesions carrying both AKT1 mutation and more than 50% of shared mutations in the MRCAs were defined as originating from the MRCA (Supplementary Table 9). In these cases, the distance between the most distant non-cancer lesions originating from the same MRCA was measured as the maximum diameter of the area over which the non-cancer clones had expanded. In der(1;16)(+) cancer cases that were not accompanied by the formation of any der(1;16)(+) non-cancer lesions, the diameter was considered to be zero. For the fresh-frozen LCM samples, we first reconstructed phylogenetic trees as described above. Next, the cell fraction of mutations allocated to private branches was estimated using PyClone-VI. The mutations in shared branches and those in private branches whose cellular prevalence was high enough to determine clonal structure on the basis of the Pigeonhole principle, were used to detect clones that existed at the age of 1 or 13 years: only when the sum of their cellular prevalence exceeded 1.0, the two mutations were thought to be in the same structures. These clones unrelated to cancer were identified only through WGS, which might lead to lower sensitivity in the detection of clonally related lesions because of the limited number of observations compared to those obtained through FISH analysis. We measured the distance between clonally related lesions and the nearest clonally unrelated lesions in these cases as the maximum diameter of the area of clonal expansion to prevent underestimation.

Analysis of significantly mutated genes

Significantly mutated genes in the pathologically normal lobules, non-cancer proliferative lesions and cancer lesions in multi-sampled cases were analysed based on dN/dS using the R package dndscv⁷⁶ (v.0.0.1.0) to investigate the difference in mutational processes during the evolution of cancer. Mutations detected in the shared branches were counted once in each group of lesions. Genes with dN/dS > 1 and q < 0.1 were considered to be significantly mutated (Supplementary Table 10).

Targeted-capture sequencing of breast cancer-associated genes

For the ten breast cancer cases in which tumour cores in the tissue microarray were defined as der(1;16)(+) by FISH analysis, tumour lesions were macro-dissected from FFPE surgical specimens and analysed by means of targeted-capture sequencing using the xGen Predesigned Gene Capture Pools (IDT) designed for 189 genes associated with breast cancer (Supplementary Table 11) and the xGen CNV Backbone Hyb Panel. Driver mutations were called using Genomon2 as described above: the variants that were not supported by (1) a sufficient depth (eight or more); (2) VAFs greater than 0.02 and variant reads of five or more and (3) EBCall P < 10⁻⁴ were excluded as sequencing artefacts and the variants that met the criteria of driver mutations used in WGS analysis were identified. In these cases, no pathogenic variants of breast cancer susceptibility genes were identified. Copy-number changes were evaluated on the basis of the sequencing data using CNACS as described above. Eight out of nine tumours successfully evaluated for copy-number changes were confirmed to carry der(1;16), detected as concurrent whole-arm 1q gain and 16q loss. The results are summarized in Supplementary Table 12.

WES and transcriptomic analysis of TCGA cohort

Germline variants and somatic mutations in TCGA breast cancer WES cohort were evaluated using Genomon2, as described for the WGS samples. Somatic mutations were evaluated for 610 sporadic cases as follows: the variants that were not supported by (1) a sufficient depth (eight or more) in both samples and matched controls; (2) VAFs greater than 0.05 and variant reads of four or more in samples; (3) VAFs below 0.05 and variant reads of two or less in matched controls; (4) a strand ratio not equal to 0 or 1; (5) Fisher’s P < 10^−1.3 and (6) EBCall P < 10⁻⁴ were excluded as sequencing artefacts; the variants within repetitive sequences were further excluded; the variants clustered within a short length (150 bp) were subjected to visual inspection on the Integrative Genomics Viewer to further eliminate sequencing errors. Copy-number changes were evaluated using CNACS, as described above. Tumours with concurrent flat whole-arm 1q gain and 16q loss were considered der(1;16) positive (Extended Data Fig. 10b). If there were any apparent breakpoints in the middle of the 1q or 16q arm, the tumours were considered der(1;16) negative.

To investigate the effect of allelic imbalance of der(1;16), that is, trisomy 1q and monosomy 16q, on gene expression, TPM values of genes on 1q or 16q in der(1;16) positive tumours were compared with those in tumours without +1q or −16q, respectively. Tumours with mutations in CDH1, CBFB and CTCF were excluded from the analysis of genes on 16q.

Estimation of sensitivity and evaluation for true positive rate in WGS mutation calling

The sensitivity of WGS mutation calling was estimated by calculating the fraction of unique heterozygous germline polymorphisms in a sample detected by means of paired analysis using another participant’s germline sample to mimic a matched control (Supplementary Note 1). The estimated sensitivity for each sample type is summarized in Supplementary Table 13.

The true positive rate of WGS mutation calling was evaluated using targeted-capture sequencing for randomly selected mutations in WGA organoid and fresh-frozen LCM samples, and randomly selected private mutations and all of the shared mutations assigned to the MRCAs in five cases for FFPE LCM samples. Libraries were prepared using the xGen Prism DNA Library Prep Kit for WGA organoid and FFPE LCM samples; extra WGS libraries were reused for fresh-frozen LCM samples because of the small amount of DNA extracted, followed by target capture using NGS Discovery Pools or xGen Custom Hyb Panel-Accel (IDT), and then sequenced on a DNBSEQ-G400RS. The sequence data were processed into .bam files through Genomon2, as described above, and each mutation was reviewed using GenomonMutationFilter, wherein mutations with sequencing depths of 100 or more in both test and germline control samples were evaluated. A mutation was considered to be validated when (1) the VAF in the test sample was five or more times higher than that in the corresponding germline control sample; (2) the VAFs in the test sample greater than or equal to 0.05 and (3) for WGA organoid samples only, (1) and (2) was achieved in both of the two paired WGA samples. The true positive rates were 96.6% for organoid WGA samples (134 out of 139 SNVs and 9 out of 9 indels), 97.3% for private mutations in FFPE LCM samples (209 out of 215 SNVs and 7 out of 7 indels) and 99.7% for fresh-frozen LCM samples (308 out of 309 SNVs and 23 out of 23 indels). As for the validation of the MRCA mutations in FFPE LCM samples, we successfully evaluated 10,190 of 10,629 mutations (95.9%); the true positive rate was 99.3% (9,657 out of 9,728 mutations) for SNVs and 98.7% (456 out of 462 mutations) for indels. The results of validation sequencing and true positive rate are summarized in Supplementary Tables 14 and 15, respectively.

Statistical analysis

Statistical analyses were performed using the R software (v.3.6.3). All P values were calculated using a two-sided analysis. Fisher’s exact test or Mann–Whitney U-test was used for group comparisons. Survival analysis for TCGA dataset was performed using the R package survival⁷⁷ (v.3.2.11), wherein the overall survival time was determined from the date of diagnosis of breast cancer to the time of last follow-up or death, and survival curves were estimated using the Kaplan–Meier product-limit method. The differences in overall survival between the patients with der(1;16)-positive and der(1;16)-negative Luminal A breast cancer were tested for statistical significance using the log-rank test. We also performed multivariate analysis using a Cox proportional hazards regression model to evaluate the effect of der(1;16) on survival when adjusted for age and stage. The results of survival analysis are summarized in Supplementary Table 16. Multiple testing was corrected based on the Benjamini–Hochberg method to compare the frequency of driver mutations in TCGA Luminal A cancer cases with and without der(1;16) (Supplementary Table 17).

Reporting summary

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