Evolutionary histories of breast cancer and related clones – Nature

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.