Articles
https://doi.org/10.1038/s41588-019-0423-x

Quantitative evidence for early metastatic
seeding in colorectal cancer
Zheng Hu1,2,3, Jie Ding1,2,3,12,13, Zhicheng Ma1,2,3,13, Ruping Sun 1,2,3, Jose A. Seoane 1,2,3,
J. Scott Shaffer3, Carlos J. Suarez4, Anna S. Berghoff5,6,7, Chiara Cremolini8, Alfredo Falcone8,
Fotios Loupakis9, Peter Birner5,10, Matthias Preusser 5,6, Heinz-Josef Lenz11 and Christina Curtis

*

1,2,3

Both the timing and molecular determinants of metastasis are unknown, hindering treatment and prevention efforts. Herewe characterize the evolutionary dynamics of this lethal process by analyzing exome-sequencing data from 118 biopsies from
23 patients with colorectal cancer with metastases to the liver or brain. The data show that the genomic divergence between
the primary tumor and metastasis is low and that canonical driver genes were acquired early. Analysis within a spatial tumor
growth model and statistical inference framework indicates that early disseminated cells commonly (81%, 17 out of 21 evaluable patients) seed metastases while the carcinoma is clinically undetectable (typically, less than 0.01 cm3). We validated theassociation between early drivers and metastasis in an independent cohort of 2,751 colorectal cancers, demonstrating their utility as biomarkers of metastasis. This conceptual and analytical framework provides quantitative in vivo evidence that systemic
spread can occur early in colorectal cancer and illuminates strategies for patient stratification and therapeutic targeting of thecanonical drivers of tumorigenesis.

M

etastasis is the primary cause of cancer-related death inpatients with cancer, but the timing and molecular determinants of this process are largely uncharacterized1–3. In
particular, when and how metastatic competence is specified are ofclinical importance. The prevailing linear progression model posits that metastatic capacity is acquired late following the gradualaccumulation of somatic alterations, such that only a subset of cells
evolve the capacity to disseminate and seed metastases4–7. However,at odds with this model, gene-expression signatures from primarytumors are predictive of distant recurrence, indicating that metastatic cells constitute a dominant subpopulation in the primarytumor8,9. In addition, disseminated tumor cells have been identifiedin patients with early breast lesions10 and in mouse models of earlybreast and pancreatic cancers11–13. However, the timing of metastatic
dissemination has not been evaluated in human cancers due to thechallenge in obtaining paired primary tumors and distant metastases and the limitations of applying phylogenetic approaches to bulk
tissue samples.
Colorectal cancer (CRC) is the third most commonly diagnosed cancer and leading cause of cancer death14, as well as a suitable model for studying tumor progression given that the initiating
driver alterations are well-characterized4. The site and resectabilityof CRC metastases dictate treatment options and prognosis15,16; theliver is the most common site of metastasis, presumably becauseof venous drainage, and one-third of patients with metastatic CRC(mCRC) have liver-exclusive metastasis16. By contrast, brain metastasis is a rare (less than 4% of mCRCs), but devastating diagnosis
with limited therapeutic options and a median survival of three to

six months17. In CRC, metastasis is assumed to be seeded by genetically advanced cancer cells that have evolved through a series ofsequential clonal expansions4,18. However, CRC progression is notnecessarily linear. Rather, we described a Big Bang model of tumor
evolution, in which after transformation some CRCs grow as a single expansion populated by heterogeneous and effectively equally fit
subclones, and from which most detectable intratumor heterogeneity arises early19. These data suggest that some CRCs may be ‘born tobe bad’, wherein invasive and even metastatic potential is specified
early19,20. Effectively neutral evolution has since been reported inother primary tumors21–24, but the mode of evolution (effective neutrality versus subclonal selection) has not been evaluated in paired
primary tumors and metastases.
Although the metastatic process is largely occult, spatio-temporal patterns of genomic variation are embedded in the evolutionary
histories of paired primary tumors and metastases. Here we analyze
exome-sequencing data from 118 biopsies from 23 patients with
mCRC who had paired distant metastases to the liver or brain todelineate the timing and routes of metastasis and to define metastasis-competent clones (Fig. 1). The data show that primary tumor–
metastasis genomic divergence (PMGD) is low and that genomic
drivers were acquired early. Moreover, through simulation studies,
we establish that low PMGD in bulk-sample sequencing data is
indicative of early dissemination, contrary to current assumptions2.
Phylogeny reconstruction and analysis of the mutational cancer
cell fraction (CCF) revealed the early divergence of metastatic
lineages and their monoclonal origin. To overcome the limitationsof phylogenetic approaches—which cannot resolve the timing of

Department of Medicine, Division of Oncology, Stanford University School of Medicine, Stanford, CA, USA. 2Department of Genetics, Stanford University
School of Medicine, Stanford, CA, USA. 3Stanford Cancer Institute, Stanford University School of Medicine, Stanford, CA, USA. 4Department of Pathology,
Stanford University School of Medicine, Stanford, CA, USA. 5Comprehensive Cancer Center CNS Tumor Unit, Medical University of Vienna, Vienna, Austria.
6
Division of Oncology, Department of Medicine I, Medical University of Vienna, Vienna, Austria. 7Institute of Neurology, Medical University of Vienna,
Vienna, Austria. 8Department of Oncology, University Hospital of Pisa, Pisa, Italy. 9Unit of Medical Oncology 1, Department of Clinical and Experimental
Oncology, Istituto Oncologico Veneto, IRCCS, Padua, Italy. 10Department of Pathology, Medical University of Vienna, Vienna, Austria. 11Department ofMedical Oncology, Norris Comprehensive Cancer Center, Keck School of Medicine, University of Southern California, Los Angeles, CA, USA. 12Present
address: Veracyte Inc, South San Francisco, CA, USA. 13These authors contributed equally: Jie Ding, Zhicheng Ma. *e-mail: cncurtis@stanford.edu
1Nature Genetics | VOL 51 | JULY 2019 | 1113–1122 | www.nature.com/naturegenetics

1113

Articles

Nature Geneticsab

P-private subclonal
P1

P-private clonal
Brain metastasis (n = 10)(72 tumor biopsies; 6 MRS)

M

P2
P3
P4
M1

Shared clonal
N

P

M2
M3

M-private clonalM

M4

M-private subclonal
M

1.0
Metastasis CCF

Liver metastasis (n = 13)(46 tumor biopsies; 3 MRS;
Uchi, Kim, Leung, Lim)P

Primary CRC

M-privateclonalSharedclonal

0.8
0.6
0.4
0.2

P-privateclonal00 0.2 0.4 0.6 0.8 1.0
Primary CCF

c

Time

Normal cell

Primary tumor
founderdLimit of detection(~10 8 cells)

Early stage CRC(n = 1,813)

Metastatic CRC(n = 938)

Primary tumor
Early driver

Late driver
Dissemination —early or late
Metastasis

Fig. 1 | Study overview. a, The cohort of patients with mCRC includes 118 tumor biopsies from 23 patients. Paired CRCs with metastases to the brain andother sites (liver, lung or lymph nodes) from 10 patients and 72 tumor biopsies were whole-exome sequenced, including 6 cases with MRS of 3–5 regions,each from the primary CRC and metastases. Additionally, four publicly available cohorts with paired CRCs and liver metastases from 13 patients and46 tumor biopsies were reanalyzed within the same bioinformatics framework, including 3 cases with MRS. b, Tumor phylogenies were reconstructedfrom somatic alterations (sSNVs and indels). The mutational CCF was computed for each primary CRC and metastasis pair. c, Schematic illustration oftumor evolution starting from a normal cell that acquires mutations leading to malignant transformation, growth of the primary tumor and metastatic
dissemination, seeding and outgrowth. It is unknown whether dissemination occurs early from a dominant subclone when the size of the primary tumor
is below the limits of clinical detection (108 cells or 1 cm3) (early dissemination) or later from a minor subclone after the acquisition of additional driveralterations (late dissemination). To address this question, we developed a three-dimensional model of tumor growth and statistical inference framework
to time metastasis from patient genomic data. d, We further leveraged a large collection of metastatic (n = 938) and non-metastatic (n = 1,813) CRCs withtargeted sequencing data to evaluate the association between specific combinations of early driver genes (modules) identified in the mCRC cohort.

dissemination2,25–28—we developed a spatial computational modelof tumor progression and Bayesian statistical inference framework
to time dissemination in a patient-specific fashion. Furthermore,
we validated the association between combinations of early drivergenes and metastasis in an independent cohort of 2,751 CRCs,
demonstrating their utility as biomarkers of aggressive disease.
These results provide quantitative in vivo evidence for early metastatic seeding in mCRC with implications for systemic therapy andearlier detection.

Results

Overview of clinical cohorts. Patients with mCRC exhibit varied
progression paths, of which liver-exclusive metastasis and brain
metastasis represent extreme scenarios with distinct prognoses15,16.
We therefore characterized the genomic landscape, routes and timing of metastasis in mCRC by analyzing exome-sequencing datafrom 118 biopsies from 23 patients with paired distant metastases to the liver or brain (referred to as the mCRC cohort, Fig. 1a,
Supplementary Fig. 1, Supplementary Table 1 and Methods). To
investigate these patterns, we sequenced 72 samples from a uniquecohort of 10 patients with mCRC who had paired brain metastases
1114

and some of whom had additional metastases to the liver (n = 1),
lung (n = 1) and lymph nodes (n = 4). Five patients had brainexclusive distant metastasis (V402, V514, V855, V953 and V974),
which is estimated to occur in only 2–10% of patients with brain
metastasis16. For six patients, multi-region sequencing (MRS) of thepaired primary tumor and metastasis (P/M pairs) was performed(3–5 regions each), enabling the detailed reconstruction of tumor
phylogenies (Fig. 1b). Additionally, we included 46 tumor biopsies
from 13 patients with mCRC who had paired liver metastases after
excluding cases with low tumor cell purity (<0.4; Supplementary
Fig. 2) from four published datasets21,29–31, analyzed using the sameunified bioinformatics framework (Methods). No other sites ofmetastasis were reported for these patients and MRS was available for 3 P/M pairs (n = 2–9 regions each). As we have previously
shown, MRS enables more accurate estimation of the cancer cell
fraction (CCF) of sSNVs and discrimination between clonal andsubclonal mutations relative to single-sample sequencing23 (Fig. 1band Supplementary Fig. 3). Additionally, we leveraged an independent collection of 2,751 patients with CRC, including 938 patientswith metastatic disease (stage IV) and 1,813 patients with earlystage disease (stages I–III) for whom targeted sequencing data fromNature Genetics | VOL 51 | JULY 2019 | 1113–1122 | www.nature.com/naturegenetics

ArticlesaNumber ofsSNVs + indels

Nature Genetics
300
200

Not shared
Shared100

Number of affected patients

0

0

5

10

15

APC
TP53
KRAS
SYNE1
SMAD4
PIK3CAAPOB
TCF7L2
DNAH8
PTPRT
AMER1
ROBO2
TGIF1
FN1
ELF3
ATM
MUC2
BRAF
PIK3R1PTEN
Mutation
Nonsynonymous
Indel
Stopgain
Splicing
Copy number
1
2
3
4
5
6
7
8

Brain metastasis CRC (nn = 10 patients)

Lim21 P
Lim21 LI

Lim16 P
Lim16 LI

Lim12 P
Lim12 LI

Lim11 P
Lim11 LI

Lim8 P
Lim8 LI

Lim7 P
Lim7 LI

2.5

Fold enrichment in
pan-cancer drivers

5.0

10.0

P < 2.2 × 10–16

7.5
P = 0.167
5.0
2.5

Pr
pr ima
iv ry
a
M te
et
as
pr tas
iv is
at
e

ar

ed

0

0

Sh

V4
V4 6 B
0 M
V5 2 B
14 M
V5 BM
V5 59
5 L
V7 9 B I
5 M
V8 0 B
2 M
V8 4 B
5 M
V9 5 BM
V9 30
3 LU
V9 0 B
5 M
V9 3 B
74 M
U B
ch M
Ki i2 Lm I
K 1
Le im2 LIu
Le ng1 LI
un L
g I
Li 2 Lm I
Li 3 Lm I
Li 6 Lm I
Li 7 Lm
Li 8 I
m LI
Li 11
m L
Li 12 Im L
Li 16 I
m LI
21
LI

0

P = 0.367

pr ima
iv ry
a
M te
et
as
pr tas
iv is
at
e

0.25

7.5

ed

0.50

P < 2.2 × 10–16

10.0

Pr

Percentage

0.75

Paired primary CRC and metastasis (nn = 25)

Lim6 P
Lim6 LIc

Metastasis privatearPrimary tumor private

1.00

Sh

SharedLiver metastasis CRC (nn = 13 patients)

Fold enrichment in CRC driversb

Lim3 P
Lim3 LI

Leung2 P
Leung2 LI

Leung1 P
Leung1 LI

Kim2 P
Kim2 LI

Kim1 P
Kim1 LI

Uchi2 P
Uchi2 LI

V974 P
V974 BM

V953 P
V953 BM

V930 P
V930 LU
V930 BM

V855 P
V855 BM

V824 P
V824 LN
V824 BM

V750 P
V750 LN
V750 BM

V559 P
V559 LI
V559 BM

V514 P
V514 LN
V514 BM

V402 P
V402 BM

0
1
2
3
4
V46 P
V46 LNV46 BM

Ploidy1p36.33−13.2
3p26.3−14.3
4
5q13.3−14.3
7p22.3−12.1
8p23.3−23.2
8p21.3−21.2
13q
14q
15q
17p
18
20q
21q

Fig. 2 | The mutational landscape and patterns of genetic divergence in paired primary CRCs and metastases. a, Concordance among somatic alterations(sSNVs, indels and CNAs) in known CRC driver genes between paired primary CRCs and metastases. Stacked bar plots illustrate the total number ofsSNVs and indels in exonic regions with a lower cut-off of variant allele frequency = 0.1 in the corresponding site (primary tumor or metastasis). BM,
brain metastasis; LI, liver metastasis; LN, lymph node metastasis; LU, lung metastasis. b, The percentage of clonal sSNVs/indels that are shared, primarytumor-private or metastasis-private out of all clonal sSNVs/indels with CCF > 60% in any of paired primary tumors and distant metastases. c, Violin plots
illustrate the probability density of driver gene fold enrichment among shared, primary tumor-private and metastasis-private clonal non-silent sSNVs/
indels based on known CRC or pan-cancer drivers. The inset box corresponds to the 25th to 75th percentile (interquartile range); the horizontal line
indicates the median; and the vertical line includes data within 1.5× the interquartile range. A test statistic was computed based on n = 100 downsamplings
among patients (Methods). P value, two-sided Wilcoxon rank-sum test. P, primary tumor; LN, lymph node metastasis; LI, liver metastasis; LU, lungmetastasis; BM, brain metastasis.

the MSK-Impact32 and GENIE33 studies were available in order toevaluate the association between specific combinations of early
driver genes (modules) defined in the mCRC cohort and metastatic
propensity (Fig. 1d and Methods).
Genomic heterogeneity in CRCs and paired metastases. Highconcordance among putative driver genes was observed in themCRC cohort (Fig. 2a), consistent with previous studies21,29,34–37. For
instance, mutations in KRAS, TP53, SMAD4, TCF7L2, FN1, ELF3
Nature Genetics | VOL 51 | JULY 2019 | 1113–1122 | www.nature.com/naturegenetics

and ATM were completely concordant between P/M pairs (Fig. 2aand Supplementary Table 2). On average, 70% of high-frequency
somatic single-nucleotide variants (sSNVs) and small insertions
and deletions (indels) with CCF > 60% (Methods) in any primarytumor or metastasis were shared by both lesions (Fig. 2b). Amonggenes that were mutated in more than five patients, SYNE1 (fourout of six patients) and APOB (three out of five patients) tended
to be private to the primary tumor or metastasis and thus likely
arose after transformation. Although metastases usually had more1115

Articles
private high-frequency sSNVs than the primary tumor (P = 0.020,
Wilcoxon rank-sum test; Fig. 2b), they were not enriched for CRC
drivers (defined based on IntOGen38 and The Cancer Genome
Atlas (TCGA)39) or a published list of pan-cancer drivers40 (Fig. 2c,
Supplementary Table 3 and Methods). Similar results were obtained
when stratifying by brain or liver metastases (Supplementary Fig. 4).
These data reflect limited driver-gene heterogeneity between P/Mpairs and suggest that few additional private genomic drivers wererequired for metastasis when the primary CRC is already advanced.
Somatic copy-number alterations (CNAs) were also generally concordant, with chromosomes 7p22.3–12.1, 13 and 20q11–13 exhibiting
recurrent amplification and chromosomes 8p23.3–23.2, 8p21.3–
21.2 and 18 exhibiting recurrent deletion in P/M pairs41 (Fig. 2aand Supplementary Fig. 5). Several putative oncogenes, including
PIK3CA, GNAS, SRC, FXR1, MUC4, GPC6 and MECOM wererecurrently (≥4 patients) amplified in metastases relative to paired
primary tumors. Notably, HTR2A (encoding 5-hydroxytryptamine
receptor 2A)—which encodes a receptor for the neurotransmitter
serotonin (which also functions as a regulatory factor in the gastrointestinal tract42)—was amplified more frequently in brain (4 out of 10)
than liver (1 out of 13) metastases (Supplementary Fig. 5).
We defined the number of metastasis-private clonal sSNVs
as Lm (merged CCF > 60% in the metastasis samples and <1%in the primary tumor samples) and the number of primary tumorprivate clonal sSNVs as Lp (merged CCF > 60% in the primary and<1% in the metastasis), where a cut-off of 60% accurately distinguished clonal and subclonal sSNVs (Fig. 1b and SupplementaryFigs. 6, 7a,b). Therefore, we used a merged CCF value of 60% as
the cut-off to distinguish clonal and subclonal mutations throughout. Brain metastases exhibited higher Lm than liver metastases (median = 24.5 compared to 9.5, respectively, P = 0.01,
Wilcoxon rank-sum test), whereas no difference was found for Lp(median = 8.5 compared to 6.0, respectively, P = 0.70, Wilcoxon
rank-sum test; Supplementary Fig. 7c), potentially reflecting longer
progression times (and more cell divisions). Neither Lm (P = 0.68,
Wilcoxon rank-sum test) nor Lp (P = 0.95, Wilcoxon rank-sum test)
differed significantly in chemotherapy-naive versus treated cases
despite a slight shift in mutational spectra (A|T>C|G) after chemotherapy (Supplementary Fig. 8).
Gene ontology analysis showed enrichment for cellular adhesion
terms among both brain and liver metastasis-private non-silent
clonal mutations, but not primary tumor-private clonal or subclonal
mutations (Supplementary Table 4). Nervous system development
and neuronal differentiation terms were enriched among brain
and liver metastasis-private clonal mutations and primary tumorprivate mutations, consistent with hijacking of the enteric nervous
system in gastrointestinal malignancies43. By contrast, primarytumor-private non-silent clonal mutations were enriched for metabolic processes, DNA repair and damage, suggestive of more general deregulation and resource constraints during tumor expansion.
Phylogenetic reconstruction of metastatic CRC. The MRS data
revealed extensive intratumor heterogeneity both within tumors
and between P/M pairs (Fig. 3a,b, Supplementary Fig. 9 andSupplementary Table 2) and ample mutations for phylogeny reconstruction. We used the FST statistic44 to quantify intratumor heterogeneity within tumors (primary tumor or metastasis) in the mCRC
cohort based on subclonal sSNVs23 (Methods). Clonal mutations
present in all samples do not contribute to intratumor heterogeneity and were excluded from FST calculations. Both the primarytumor (median FST = 0.180, range = 0.150–0.430) and paired metastases (median FST = 0.178, range = 0.123–0.271) exhibited high FST
values, consistent with rapid genetic diversification (Supplementary
Fig. 10a). Proliferative indices based on Ki-67 staining were alsosimilar between paired CRCs and metastases (P = 0.765, Wilcoxon
signed-rank test, Supplementary Fig. 10b).
1116

Nature Genetics
Tumor phylogenies were reconstructed using sSNVs and indels
across multiple regions of each P/M pair using the maximum parsimony method45. Distant metastases corresponded to monophyleticclades in all but one (Kim1) case (eight out of nine cases with MRS;
Fig. 3c, Supplementary Fig. 9 and Methods), consistent with theunique origin of the metastatic lineage. Inspection of the phylogeny
for Kim1 indicated that the liver metastasis preceded the primarytumor, which is improbable and likely due to metastasis-specificloss of heterozygosity (LOH) spanning multiple mutations. In mostpatients, the metastatic lineage diverged before genetic diversification of the primary tumor (V402, V930, V953, V974 and Uchi2;
early divergence), whereas divergence occurred during diversification of the primary tumor in patients V750, V824 and Kim2
(late divergence). All brain metastases and most liver metastases
contained many private clonal sSNVs, but lacked shared subclonalsSNVs with the primary tumors, consistent with monoclonal seeding (Supplementary Figs. 11, 12), as demonstrated by simulation
studies (Supplementary Fig. 13). Two liver metastases (Lim6 andLim11) exhibited enrichment of shared subclonal mutations, butlacked metastasis-private clonal mutations, consistent with polyclonal seeding (Supplementary Figs. 12, 13). These data suggest that
distant metastases are often seeded by a single clone (a single cell or agroup of genetically similar cells). Notably, the phylogenetic tree forcase V930 indicates that the brain metastasis derived from the lungmetastasis, in line with the clinical history of the patient (Fig. 3).
Brain metastases and regional lymph node metastases formed
separate clades in the two cases in which they were profiled (V750
and V824), indicative of their independent clonal origin from theprimary tumor (Fig. 3c and Supplementary Fig. 9) and consistent
with polyguanine-repeat analysis46.
The finding that paired CRCs and metastases formed separate
phylogenetic clades in most patients suggests that metastatic dissemination may occur early during cancer development, such thatthe primary tumor has sufficient time to accumulate many uniqueclonal mutations after dissemination. However, phylogenetic divergence may occur much earlier than dissemination (Supplementary
Fig. 14) and phylogenetics cannot resolve the timing of dissemination2,25–28. As such, we next investigated the determinants of PMGD
and quantified the timing of metastasis.
The timing of dissemination and PMGD. To model the evolutionary dynamics of metastasis, we developed a three-dimensional
agent-based computational model to simulate the spatial growth,
progression and lineage relationships of realistically sized patient
tumors under varied parameters19,23 (Fig. 4a, Supplementary Fig. 15,
Supplementary Table 5 and Methods). We modeled the growth of a
primary CRC starting from a single founder cell and assumed thatthe metastasis was seeded by a random single cell from the periphery of the primary tumor, yielding primary and metastatic tumors
composed of approximately 109 cells (around 10 cm3). To accountfor distinct modes of tumor evolution, we simulate effective neutrality and stringent subclonal selection19,23, resulting in four evolutionary scenarios for P/M pairs: neutral/neutral (N/N), neutral/selection (N/S), selection/neutral (S/N) and selection/selection(S/S) (Fig. 4a, Supplementary Figs. 15, 16 and Methods). Using thissimulation framework, for which ground-truth values are known,
we evaluated the relationship between the number of metastasisprivate clonal sSNVs (Lm) and the primary CRC size at the time ofdissemination (Nd) in hundreds of virtual paired P/M tumors, forwhich size is a surrogate measure for time, as cell division rates are
unknown (Methods).
To define Lm, we first evaluated metastasis-private clonalsSNVs with relatively high-frequency sSNVs in the whole primarytumor (CCF > 1%). Therefore, any clonal sSNVs in the metastasis
will be private to the metastasis if CCF < 1% in the primary tumor.
We found that Lm is positively correlated with Nd under all four
Nature Genetics | VOL 51 | JULY 2019 | 1113–1122 | www.nature.com/naturegenetics

Articles
8

ye
ar
s,
4
V402

No evidence of disease

c

ERBB3

MED24

AMER1

ye
ar
s
(d , 8
ea m
th on
)
th
s

V402 P 3m

V402 P 1

7

on
th
s

b

ATRX

a

WIPF1
APOB
KRAS
APC
PTEN

Nature Genetics

ATRX

V402 P 2

P1
BM 1

x
,S
BM

ERBB3

BM 3

BM 4

PARK2

VAF ≥ 0.1
VAF < 0.1
P2

PARK2

APC, KRAS, TP53,
PIK3CA, TCF7L2, EYS

P3
P1
BM 4
BM 3
BM 2

Normal
133

20.0

84

65
PTPRT
RP1L1

No WBRT

BM 5

Early
divergence

BM 1

BM

FN1
MSR1

Normal20.0
5 3 57

162

SMAD4
TP53
SYNE1
KRAS
APC
PIK3R1
APC
APC
SYNE1
DNAH8
APOB
APOB
LIG1
PPP2R1A

x
,S
BM

x
P
LN , S
,Sx
x
L
LUI

D

Disease progression

1(d 0 m
ea on
th th
)
s

s
th
on
m
6
Stable disease
Chemotherapy, capecitabine
(5FU)

P2
AKAP9

P5
P4
LU 2

MSR1, FN1

61

V750

P3

P1
APC, KRAS,SYNE1, TCF7L2

56

V750 P 1
V750 P 2
V750 P 3
V750 P 4
V750 P 5
V750 LN 1
V750 LN 2
V750 LN 3
V750 BM 1
V750 BM 2
V750 BM 3
V750 BM 4
V750 BM 5

LU 3

PTPRT

Early
divergence

LU 1
BM 5
BM 2
BM 3
BM 4
BM 1

48

SYNE1

,S

,S
LU

D

V930 P 1
V930 P 2
V930 P 3
V930 P 4
V930 P 5
V930 LU 1
V930 LU 2
V930 LU 3
V930 BM 1
V930 BM 2
V930 BM 3
V930 BM 4
V930 BM 5

CYLC1
CONT1
RP1L1

hs
ar
ye

10

ye
8

Diseaseprogression

x

x

x
,S
LI

(d s, 6
ea m
th o
) nt

s
th
on
m
7
s,

ye
4
No
evidenceof disease

Chemotherapy
(5FU)

x

P,

ar

ye
2
No
evidenceof disease

Sx

V930

ar

ar

s,

s,

8

2

m

m

on

on

th

th

s

s

TCF7L2
APC
KRAS
APC
SYNE1

R

No WBRT

EYS
APC

TP53
TCF7L2APCKRASPIK3CA1(d 3 m
ea on
th th
)
s

on
th
s
m
5
x
,S
BM

No chemotherapy

ec
ur
rin
g

x
P
BM , S,S xx

D

No chemotherapy
y

144

216

V974 P 1
V974 P 2
V974 P 3
V974 BM 1
V974 BM 2
V974 BM 3
V974 BM 4
V974 BM 5

No gamma knife

MED24
Early
divergence

20.0

AKAP9

Sx
P,
x
D

V402 BM 4

BM to local
reccurence

BM 2

NormalV402 BM 3
67

No gamma knife

P4

V402 BM 1
V402 BM 2

V974

P3

V402 P 4

Chemotherapy

P2

APC, KRAS,
PTEN, APOB,
WIPF1, AMER1

APC, KRAS,
TP53, SMAD4,
SYNE1, PIK3R1

P3
PPP2R1A
LIG1Normal

APOB

LN 3
LN 2
LN 1
P4
P2
P1
P5

RP1L1,
CNOT1,
CYLC1

BM 5

BM 2

BM 3
20.0
116

18 21

115

75

287

BM 1
Late
divergenceFig. 3 | Within- and between-lesion heterogeneity in paired primary CRCs and metastases. a, Clinical and treatment history for four representative
patients who had CRC with brain metastases. Dx, diagnosis; Sx, surgical resection. b, Patterns of within- and between-lesion heterogeneity among sSNVs
and indels based on MRS of paired primary CRCs and metastases, for which canonical CRC driver genes are labeled. The number of mutations thatare shared or private among different lesions is indicated below as the corresponding colored horizontal bars: ubiquitously P/M shared (red), partially
P/M shared (green indicates M1; blue indicates M2), primary tumor-private (pink) or metastasis-private (yellow indicates M1 and gray indicates M2;
or cyan indicates M1 and M2). P corresponds to primary tumor. M1 and M2 correspond to different metastatic sites in the same patient when multiple
metastatic sites were sampled (V974: M1-LU, M2-BM; V750: M1-LN, M2-BM). c, Phylogeny reconstruction using maximum parsimony (PHYLIP) basedon mutational presence or absence, for which canonical CRC drivers genes are labeled. VAF, variant allele frequency. WBRT, whole-brain radiation therapy.
5FU, 5-fluorouracil.

evolutionary scenarios (Fig. 4b). The positive relationship between
Lm and Nd remains significant when accounting for variation in
mutation rate, cell birth and death rate, and selection intensity during tumor growth (Supplementary Fig. 17). We next evaluated Lmby simulating sequencing reads from variable numbers of primarytumor regions (n = 1, 10, 50 or 100) while considering the wholemetastasis as a bulk sample within our computational model. The
positive correlation between Lm and Nd was highly significant
under all sampling scenarios, pointing to the robustness of thisobservation (Supplementary Fig. 18). As expected, smaller Lm wasobserved when a greater number of primary tumor regions weresequenced because fewer mutations were private to the metastasis(Supplementary Fig. 18). Mathematical analysis of the special caseof neutral evolution and exponential growth further demonstrates
Nature Genetics | VOL 51 | JULY 2019 | 1113–1122 | www.nature.com/naturegenetics

the positive relationship between Lm and Nd (Supplementary Note,
Supplementary equation (6)). These data suggest that later dissemination results in more clonal mutations in the metastasis, manyof which are at low frequency in the primary tumor and are often
undetectable in bulk sequencing. Accordingly, later dissemination
will give rise to more metastasis-private clonal mutations in real
sequencing data, leading to higher PMGD. It should be noted that
if sampling of the primary tumor was exhaustive or if the metastasis-founder clone could be traced—neither of which are generally practical for studies of tumors in human patients—one would
expect very small Lm values and no correlation between Lm and Nd
since all mutations in the metastasis-founder cell that accumulated
during primary tumor growth would be captured. By contrast, thenumber of primary tumor-private clonal sSNVs (Lp) exhibited a1117

Articles

Nature GeneticsaEffectively neutral

Subclonal selection

Time
Normalcell

TimePrimary tumor
founder

Beneficial
subclones

Primary tumor
founder

Normalcell

Metastasis

b

N/N
160

N/S
160

r = 0.90

S/N160

r = 0.86

80

80

Lm

Metastasis

0

60

3

4

5

6

7

8

9

3

4

5

6

7

8

0
2

9
60

r = –0.01

r = 0.72

80

0
2

60

r = –0.03

r = 0.74

80

0
2

S/S
160

3

4

5

6

7

8

9

2
60

r = –0.38

40

40

40

20

20

20

20

4 5

6

7

8

9

6

7

8

9

5 6 7
log10(NNd)

89r = –0.36

Lp

40

3

0

0
2

H = Lm/(Lp + 1)

160

3

4

5

6

7

8

9
160

r = 0.85

80

0
2

3

4

5

6

7

8

9
160

r = 0.83

80

0
3

4

5 6 7
log10(N d)

8

9

3

4

5

6

7

8

9

3

4

5

6

7

8

9

log10(N d)

3

4

5

r = 0.69

80

0
2

2
160

r = 0.71

80

0
2

0
2

0
2

3

4

5 6 7
log10(N d)

8

9

2

3

4

Fig. 4 | Correlation between the Lp, Lm and H and primary carcinoma size at the time of dissemination. a, Schematic illustration of effectively neutral (N)evolution and stringent subclonal selection (S), two distinct evolutionary modes that can occur during the growth of the primary tumor or metastasis. It
is assumed that metastatic dissemination occurs during expansion of the primary CRC following malignant transformation, where Nd corresponds to thesize of the malignant clone (carcinoma) at the time of dissemination. b, The correlation between the timing of dissemination and Lm, Lp or H, based on thespatial simulation of tumor growth (n = 100 tumors for each scenario; Pearson’s r is reported). Lp and Lm correspond to the number of private clonal sSNVs(CCF > 60% in one site and CCF < 1% in the other site) in the whole primary carcinoma and metastasis, respectively, and H = Lm/(Lp + 1).

slightly negative correlation with Nd when CRCs grew under stringent selection (S/N or S/S), whereas under neutral evolution (N/Nor N/S) Lp ≈ 0, regardless of the timing of dissemination (Fig. 4b andSupplementary Fig. 17).
We defined early dissemination as Nd < 108 cells (around 1 cm3
in volume)—the size at which CRCs are generally clinically detectable—and late dissemination as Nd ≥ 108 cells. To establish intuition
for the relationship between PMGD and Nd, we defined H = Lm/(Lp + 1). In the simulation studies, H was positively correlated withNd (Fig. 4b and Supplementary Fig. 17), indicating that larger H
values are associated with later dissemination. Indeed, late dissemination typically results in large H (>20) (Fig. 4b). The observationthat most patients in the mCRC cohort exhibited small H values(median = 2.4, range = 0.5–23.5) suggests that early dissemination
may be relatively common. Although H is strongly associated withthe timing of dissemination, it does not capture all componentsof PMGD, including the mutation rate, as this is cancelled out inthe division of Lm over Lp. Additionally, variation in Lp due to differences in the mode of evolution and sampling bias contribute to1118

noise in H. To account for these sources of variability while estimating the timing of dissemination in individual patients, we turned toa powerful statistical inference framework grounded in population
genetics theory.
Quantitative evidence for early metastatic seeding in CRC. In
order to infer the timing of dissemination Nd, mutation rate u (percell division in exonic regions) and mode of tumor evolution in P/Mpairs, we developed SCIMET (spatial computational inference ofmetastatic timing), which couples our spatial (three-dimensional)
agent-based model of tumor evolution with a statistical inference
framework based on approximate Bayesian computation (ABC)47,48(Fig. 4a, Supplementary Figs. 15, 16, 19, Supplementary Tables 6,
7 and Methods). The use of ABC is well-established in population
genetics and has been utilized to infer the parameters of tumor
evolution19,49. As the patient genomic data were generally consistent with monoclonal seeding, we assumed that a single cell seeds
the metastasis (Lim6 and Lim11 were therefore excluded from thisanalysis). Evaluation of SCIMET on virtual tumors demonstrates
Nature Genetics | VOL 51 | JULY 2019 | 1113–1122 | www.nature.com/naturegenetics

Articles

Nature Geneticsa10

r = 0.63; P = 0.001

6

4

10 5

10 6

10 7

10 8

×

×

×

×

1

1

1

1

Probability

0.2
0

1
r = 0.75; P = 4.4 × 10–520

8~
log10(NNd)~
log10(NNd)

825

r = 0.42; P = 0.046

0.4

1 cm3

0.01

H = Lm /(Lp + 1)

b 10

0.6

10 9

10 4
×

1

0.0001

0.8

×

10 3
×
1

5

0
3.

6

1.

3

0.

00

00

0.

0.

Evolution P/Mmode (P/M) pair[3, 4][4, 6][4, 6][4, 6][4, 6][4, 6][4, 6][4, 7][5, 7][4, 8][7, 8]

0.

Kim1 LI
Leung2 LI
Lim3 LI
Lim7 LI
Lim8 LI
Lim12 LI
Lim16 LI
Uchi2 LI
Lim21 LI
Leung1 LI
Kim2 LI

6
01
5
0.
03
0.
06
0.
15

Liver metastasis

N/N
S/S
S/N
S/N
S/S
S/S
S/S
S/S
S/S
S/S
S/N

0.

V824 BM
V953 BM
V974 BM
V402 BM
V514 BM
V855 BM
V46 BM
V559 LI
V559 BM
V750 BM
V930 LU
V930 BM

3

Brain metastasis

N/N
S/N
S/N
S/N
S/S
S/N
S/S
S/S
S/S
S/N
S/S
S/S

log10-transformed(lower bound, upper bound)[3, 5][3, 6][3, 6][4, 6] Early[4, 6][4, 6][4, 7][4, 7][4, 7][6, 8] Late[6, 9][6, 9]

Nd

u

6

4

15
10
5
0

0

5

10

15

H = Lm /(Lp + 1)

20

25

0

25

50

75

100

Time elapsed from Dx of P to M
(months)

0

25

50

75

100

Time elapsed from Dx of P to M
(months)Fig. 5 | Patient-specific inference of the timing of metastasis in CRC. a, Heat map of the posterior probability distributions inferred by SCIMET forthe mutation rate u (per cell division in exonic regions) and Nd (timing of metastatic dissemination relative to primary carcinoma size) in individual∼
P/M pairs (n = 23) from 21 patients with mCRC. The median of the posterior distribution (Nd) is indicated by a white circle at the corresponding value.
For patients with more than one distant metastasis, each metastasis was analyzed independently. The mode of tumor evolution in each P/M pair wasdetermined based on model selection within the statistical inference framework (Methods). We define early dissemination as Nd (upper bound) < 108 cells(approximately 1 cm3 in volume) and use the third quartile of the posterior distribution as the upper bound to be conservative. Late dissemination isdefined as Nd (upper bound) ≥ 108. P/M pairs for which dissemination and seeding are inferred to have occurred early are denoted in blue, whereas those∼inferred to have disseminated late are denoted in magenta. b, Correlations between Nd based on SCIMET and the H metric as well as the time elapsed
from diagnosis of the primary tumor to diagnosis of the metastasis (n = 23). The Pearson’s r and P values are reported. Shading corresponds to the 95%
confidence interval of the linear regression.

the accurate recovery of the mutation rate and timing of dissemination (Supplementary Fig. 20).
The majority (90%) of CRCs and metastases (57%) exhibited
patterns consistent with subclonal selection (Fig. 5a). Inference of
patient-specific mutation rates using SCIMET showed an orderof magnitude variation across patients (inferred u or ũ = 0.06–
0.6, corresponding to 10−9–10−8 mutations per base pair per cell
division). Notably, in 83% (19 out of 23) P/M pairs from 17 out of
21 patients, dissemination was estimated to occur early when the
primary CRC was below the limits of clinical detection (inferred
Nd or N͠ d < 108 cells) and typically when the primary tumor wascomposed of fewer than 106 cells using conservative estimates(Fig. 5a and Methods). The N͠ d values were also significantly smaller
than the tumor size documented at the time of diagnosis in thisNature Genetics | VOL 51 | JULY 2019 | 1113–1122 | www.nature.com/naturegenetics

cohort (Supplementary Table 1). Of note, early dissemination
was common irrespective of the site of distant metastasis (8 out
of 10 brain and 10 out of 12 liver). Congruent results wereobtained when accounting for higher ratios of cell birth and death
rates in the primary CRC and metastasis (Supplementary Fig. 21),the collective dissemination of small clusters of cells (n = 10 cells;
Supplementary Fig. 22) or single-region sampling (Supplementary
Fig. 23). Among the four cases for which late dissemination wasinferred, three had MRS data, enabling comparison with their phylogenies. For two patients (V750 brain metastasis and Kim2 liver
metastasis), late dissemination was consistent with the tumor phylogeny (Fig. 3c and Supplementary Fig. 9). For patient V930, latedissemination was inferred for both the lung and brain metastases,
consistent with the large H values (brain, H = 23.5; lung, H = 11).
1119

Articles

Nature Geneticsa
APC/KRAS/
TP53/SMAD4
q = 7.6 × 10–17

0.8

82.7%

Fraction of patients

66.5%

Early-stage CRC (nn = 1,813)

APC/KRAS/TP53/
SMAD4 and X
q = 2.9 × 10–20
18.0%

0.15

0.6

APC, TP53, TCF7L2
q = 0.028

0.04
0.03
0.02
0.01
0

34

0.10

0.4

0.020
0.015
0.010
0.005
0

40

APC, KRAS, PTPRT
q = 0.00015

7

0.025
0.020
0.015
0.010
0.005
0

22

mCRC (nn = 938)

APC, KRAS, AMER1
q = 0.014

APC, TP53, AMER1
q = 0.00068

KRAS, TP53, PTPRT
q = 0.0004

APC, TP53, PTPRT
q = 6.0 × 10–12

0.03
0.02
0.01
12

23

0

SMAD4, TP53, PTPRT
q = 0.047

1

32

APC, TP53, KRAS, PTPRT
q = 0.0028

5.6%

0.05

0.2

0

1,206

776

102

0

b

169

0.010

0.010

0.005

0.005

0

0

Colon polyps/dysplasia

13

0.009

0.006

0

0

0.004

0.006

0.002

0.003

0

12

0

7

Early-stage CRC

Transformation(origin of malignant clone)

Driver 1

Driver 1, 2

Driver 1, 2, 3

Premalignant evolution(selection and clonal successions)

010Late-stage CRC
Early dissemination(malignant clone size < 1 cm3)

Tumor
initiation

0

Late dissemination

Driver module
1, 2, 3, X

Driver 1, 2, 3, X

Expansion ofmalignant clone

Diversification anddissemination ofmalignant clone

Early detection via cfDNA/surgical resection
Premalignant cells

Malignant cells

Adjuvant chemotherapy

Chemotherapy

Targeted and immunotherapies

Fig. 6 | Enrichment of early driver gene modules in mCRC and clinical implications of early dissemination. a, The enrichment of canonical core CRC
driver genes (APC, KRAS, TP53 or SMAD4) plus recurrent mutations in candidate drivers (AMER1, ATM, BRAF, PIK3CA, PTPRT or TCF7L2) identified in themCRC cohort was evaluated in an independent cohort of 2,751 patients with CRC. The combined bar plots (left) illustrate the overall frequency of the coremodule alone or with an additional candidate driver (X) in early-stage CRCs versus mCRCs. Individual bar plots indicate the frequency of specific modules.
q values are based on two-sided Fisher’s exact tests with Benjamini–Hochberg adjustment. b, Three stages of CRC progression are outlined: premalignancy(between initiation and transformation), early-stage CRC (between transformation and dissemination) and late-stage CRC (after dissemination). A setof potential interventions to prevent cancer mortality targets each stage could be implemented: for premalignant lesions, resection (after detection bycolonoscopy or possibly cell-free DNA (cfDNA)); for early-stage CRC, surgical resection and possibly adjuvant chemotherapy; and for late-stage CRC,
chemotherapy and/or targeted/immune therapies. Given the high rate (80% here) of early dissemination, before clinical detectability of the early-stage
CRC, detection and resection of premalignant lesions will have the greatest impact on preventing cancer mortality. For tumors that undergo dissemination
before clinical detectability, surgical resection alone, even of a small tumor, cannot prevent metastasis. Once the early-stage tumor is discovered, newly
defined metastatic modules (a) may inform patient stratification to aid the directed use of adjuvant chemotherapy.

However, the tumor phylogeny indicates early divergence of themetastatic lineage (Fig. 3c). This case illustrates that phylogenetic
divergence can occur before dissemination (Supplementary Fig. 14),emphasizing the need for a quantitative evolutionary framework totime metastasis.
The N͠ d values based on SCIMET were positively correlatedwith H (Pearson’s r = 0.63, P = 0.001; Fig. 5b), consistent with theobservation that the H metric reflects the timing of dissemination.
Additionally, both N͠ d and H were positively correlated with the time
elapsed between diagnosis of the primary CRC and distant metastasis (Fig. 5b), suggesting that metastases that are diagnosed later
likely disseminated later. Furthermore, we estimated the time span
between metastatic dissemination and surgical resection of the primary tumor using an approximate analytical function for our spatial
tumor growth model and found that dissemination often occurred
more than three years before surgery (Supplementary Fig. 24 andSupplementary Note).
1120

Metastasis-associated early driver gene modules. As noted above,
most canonical drivers were clonal and shared between paired primaries and metastases (Fig. 2), indicative of their early acquisition
before transformation. Taken together with the finding that cancer
cells seed metastases early in the majority of mCRCs in this cohort,specific combinations of early driver genes (modules) may confer
metastatic competence. In support of this view, oncogene engineering of four canonical early driver genes (APC, KRAS, TP53 andSMAD4) in wild-type primary colon organoids yielded metastases
after xenotransplantation50. Similarly, in a mouse model of CRC,
oncogenic Kras in combination with Apc and Trp53 deficiency was
sufficient to drive metastasis51.
We therefore evaluated the association between the early driver
modules defined in the mCRC cohort and metastatic proclivity
by analyzing a collection of 2,751 patients with CRC, including
938 patients with metastatic disease (stage IV) and 1,813 patients withearly-stage CRC (stages I–III) that were prospectively sequenced as
Nature Genetics | VOL 51 | JULY 2019 | 1113–1122 | www.nature.com/naturegenetics

Articles

Nature Genetics
part of the MSK-Impact32 and GENIE33 studies (Methods). Notably,
we find that numerous early driver gene modules were significantly
enriched in metastatic relative to early-stage CRCs in this independent dataset after correction for multiple-hypothesis testing (Fig. 6a,
Supplementary Fig. 25, Supplementary Table 8 and Methods). These
modules consist of a backbone of canonical core CRC drivers (combinations of APC, KRAS, TP53 or SMAD4) with one additional candidate metastasis driver (TCF7L2, AMER1 or PTPRT). Collectively,
the core modules plus an additional candidate metastasis drivershow statistically significant enrichment in metastatic versus earlystage CRCs (18% compared to 5.6%, respectively, q = 2.9 × 10−20).
Examination of the prevalence and enrichment of individual modules indicates that PTPRT mutations in combination with canonicaldrivers were almost exclusively observed in patients with metastases(Fig. 6a and Supplementary Fig. 25). Thus, PTPRT appears to be ahighly specific driver of metastasis. PTPRT mutations were previously reported in 26% of CRCs52 and loss of PTPRT in CRC and inhead and neck squamous cell cancers results in increased STAT3
activation and cellular survival53,54. It has therefore been proposedthat PTPRT mutations may be predictive biomarkers for STAT3
pathway inhibitors, highlighting new therapeutic opportunities54.
Other modules, which involved AMER1 and TCF7L2, were alsosignificantly enriched in metastatic cases, but were less specific;
perhaps because an additional driver defines the module. We therefore identify a compendium of metastasis driver modules that can
inform the stratification and therapeutic targeting of patients withaggressive disease.

Whereas drivers were not enriched in metastases when all cases wereconsidered (Fig. 2c), stratifying by the mode of tumor evolution
revealed the enrichment of private high-frequency (CCF > 20%)
driver mutations in metastases evolving under stringent selection
compared to those evolving neutrally (Supplementary Fig. 26),
suggesting that further subclonal driver mutations may occur
during the growth of some metastases. Nevertheless, a sizeable
proportion (43%) of distant metastases evolved neutrally, potentially reflecting the high fitness of the metastatic clone, consistent
with a fitness plateau56.
The finding that early dissemination—which results in successful metastatic seeding—can occur before the primary tumor isclinically detectable in the majority (80%) of patients with mCRC
in this cohort underscores the importance of detecting malignancy
at the earliest possible stage (Fig. 6b). Such small tumors fall below
the detection limits for current imaging modalities, but advancesin profiling circulating cell-free tumor DNA may ultimately enableearlier non-invasive detection57,58. Importantly, a considerable number of patients with mCRC did not exhibit early systemic spread,suggesting that colonoscopy can be beneficial in this subgroup. Ourdata also raise the possibility that patients with early-stage disease
with combinations of driver genes that confer a high risk of metastasis may particularly benefit from adjuvant chemotherapy to target micro-metastatic disease59. Although the clinical utility of thisapproach needs to be prospectively evaluated, our findings provide
a rationale for patient stratification and therapeutic targeting.

Discussion

Any methods, additional references, Nature Research reporting
summaries, source data, statements of code and data availability andassociated accession codes are available at. https://doi.org/10.1038/
s41588-019-0423-x.

We describe a theoretical and analytical framework that yields quantitative in vivo measurements of the dynamics of metastasis in apatient-specific manner, while accounting for confounding factors,
including the metastasis founder event, the mode of tumor evolution, mutation rate variation and tissue sampling bias. By analyzing genomic data from paired primary CRCs and distant metastases
to the liver and brain from five patient cohorts within this evolutionary framework, we demonstrate that metastatic seeding often
occurs early (17 out of 21 patients), when the carcinoma is clinically undetectable (~104–108 cells or 0.0001–1 cm3) and years before
diagnosis and surgery (Fig. 5 and Supplementary Figs. 21–24). The
observation that early metastatic seeding was prevalent irrespectiveof the site of distant metastasis, indicates the generalizability of theseresults. Moreover, dissemination was early even when considering
liver-exclusive and brain-exclusive metastases, which represent
extremes in terms of their prevalence and prognosis. Collectively,
these findings indicate that CRCs can be ‘born to be bad’, for whichinvasive and metastatic potential is specified early19,20,55, illuminating
the need to target the canonical drivers of tumorigenesis. However,
not all tumors will metastasize and there is an urgent need to identify biomarkers that are associated with aggressive disease.
Towards this end, we validated metastasis-associated driver modules in an independent cohort, thus defining the molecular featuresof metastasizing clones. The overlap with drivers of initiation andcombinatorial structure of these modules may explain why few drivers of metastasis have been identified to date. Although the canonical driver landscape is relatively sparse, there are nonetheless many
possible combinations of mutations that collectively disrupt key signaling pathways (WNT, TP53, TGFB, EGFR and cellular adhesion)
enabling niche independence and outgrowth at foreign sites50.
Of note, the vast majority (90%) of primary tumors in the mCRC
cohort exhibited subclonal selection consistent with the metastatic
clone having a selective growth advantage (Fig. 5a). By contrast, asmaller proportion of early stage (I–III) CRCs (33%) exhibited patterns consistent with subclonal selection23, suggesting that the modeof tumor evolution may correlate with disease stage or aggressiveness, although larger studies are needed to investigate this trend.
Nature Genetics | VOL 51 | JULY 2019 | 1113–1122 | www.nature.com/naturegenetics

Online content

Received: 21 December 2018; Accepted: 18 April 2019;
Published online: 17 June 2019

References

1. Vanharanta, S. & Massagué, J. Origins of metastatic traits. Cancer Cell 24,
410–421 (2013).
2. Turajlic, S. & Swanton, C. Metastasis as an evolutionary process. Science 352,
169–175 (2016).
3. Lambert, A. W., Pattabiraman, D. R. & Weinberg, R. A. Emerging biologicalprinciples of metastasis. Cell 168, 670–691 (2017).
4. Jones, S. et al. Comparative lesion sequencing provides insights into tumor
evolution. Proc. Natl Acad. Sci. USA 105, 4283–4288 (2008).
5. Campbell, P. J. et al. The patterns and dynamics of genomic instability inmetastatic pancreatic cancer. Nature 467, 1109–1113 (2010).
6. Yachida, S. et al. Distant metastasis occurs late during the genetic evolutionof pancreatic cancer. Nature 467, 1114–1117 (2010).
7. Yates, L. R. et al. Genomic evolution of breast cancer metastasis and relapse.
Cancer Cell 32, 169–184 (2017).
8. van de Vijver, M. J. et al. A gene-expression signature as a predictor ofsurvival in breast cancer. N. Engl. J. Med 347, 1999–2009 (2002).
9. Ramaswamy, S., Ross, K. N., Lander, E. S. & Golub, T. R. A molecularsignature of metastasis in primary solid tumors. Nat. Genet. 33, 49–54 (2003).
10. Sänger, N. et al. Disseminated tumor cells in the bone marrow of patientswith ductal carcinoma in situ. Int J. Cancer 129, 2522–2526 (2011).
11. Husemann, Y. et al. Systemic spread is an early step in breast cancer. Cancer
Cell 13, 58–68 (2008).
12. Rhim, A. D. et al. EMT and dissemination precede pancreatic tumor
formation. Cell 148, 349–361 (2012).
13. Hosseini, H. et al. Early dissemination seeds metastasis in breast cancer.
Nature 540, 552–558 (2016).
14. Siegel, R. L. et al. Colorectal cancer statistics, 2017. CA Cancer J. Clin. 67,
177–193 (2017).
15. Andres, A. et al. Surgical management of patients with colorectal cancer andsimultaneous liver and lung metastases. Br. J. Surg. 102, 691–699 (2015).
16. Vatandoust, S., Price, T. J. & Karapetis, C. S. Colorectal cancer: metastases toa single organ. World J. Gastroenterol. 21, 11767–11776 (2015).
17. Christensen, T. D., Spindler, K. L., Palshof, J. A. & Nielsen, D. L. Systematic
review: brain metastases from colorectal cancer—incidence and patientcharacteristics. BMC Cancer 16, 260 (2016).
1121

Articles
18. Fearon, E. R. & Vogelstein, B. A genetic model for colorectal tumorigenesis.
Cell 61, 759–767 (1990).
19. Sottoriva, A. et al. A Big Bang model of human colorectal tumor growth.
Nat. Genet. 47, 209–216 (2015).
20. Ryser, M. D., Min, B. H., Siegmund, K. D. & Shibata, D. Spatial mutation
patterns as markers of early colorectal tumor cell mobility. Proc. Natl Acad.
Sci. USA 115, 5774–5779 (2018).
21. Uchi, R. et al. Integrated multiregional analysis proposing a new model ofcolorectal cancer evolution. PLoS Genet. 12, e1005778 (2016).
22. Suzuki, Y. et al. Multiregion ultra-deep sequencing reveals early intermixing
and variable levels of intratumoral heterogeneity in colorectal cancer.
Mol. Oncol. 11, 124–139 (2017).
23. Sun, R. et al. Between-region genetic divergence reflects the mode and tempo
of tumor evolution. Nat. Genet. 49, 1015–1024 (2017).
24. Bozic, I., Gerold, J. M. & Nowak, M. A. Quantifying clonal and subclonalpassenger mutations in cancer evolution. PLoS Comput. Biol. 12,
e1004731 (2016).
25. Hong, W. S., Shpak, M. & Townsend, J. P. Inferring the origin of metastases
from cancer phylogenies. Cancer Res. 75, 4021–4025 (2015).
26. Naxerova, K. & Jain, R. K. Using tumour phylogenetics to identify the rootsof metastasis in humans. Nat. Rev. Clin. Oncol. 12, 258–272 (2015).
27. Zhao, Z. M. et al. Early and multiple origins of metastatic lineages within
primary tumors. Proc. Natl Acad. Sci. USA 113, 2140–2145 (2016).
28. Schwartz, R. & Schaffer, A. A. The evolution of tumour phylogenetics:
principles and practice. Nat. Rev. Genet. 18, 213–229 (2017).
29. Kim, T. M. et al. Subclonal genomic architectures of primary and metastatic
colorectal cancer based on intratumoral genetic heterogeneity. Clin. Cancer
Res. 21, 4461–4472 (2015).
30. Leung, M. L. et al. Single-cell DNA sequencing reveals a late-dissemination
model in metastatic colorectal cancer. Genome Res. 27, 1287–1299 (2017).
31. Lim, B. et al. Genome-wide mutation profiles of colorectal tumors andassociated liver metastases at the exome and transcriptome levels. Oncotarget
6, 22179–22190 (2015).
32. Yaeger, R. et al. Clinical sequencing defines the genomic landscape ofmetastatic colorectal cancer. Cancer Cell 33, 125–136 (2018).
33. The AACR Project GENIE Consortium AACR Project GENIE: powering
precision medicine through an international consortium. Cancer Discov. 7,
818–831 (2017).
34. Lee, S. Y. et al. Comparative genomic analysis of primary and synchronous
metastatic colorectal cancers. PLoS One 9, e90459 (2014).
35. Xie, T. et al. Patterns of somatic alterations between matched primary andmetastatic colorectal tumors characterized by whole-genome sequencing.
Genomics 104, 234–241 (2014).
36. Brannon, A. R. et al. Comparative sequencing analysis reveals high genomic
concordance between matched primary and metastatic colorectal cancer
lesions. Genome Biol. 15, 454 (2014).
37. Tan, I. B. et al. High-depth sequencing of over 750 genes supports linear
progression of primary tumors and metastases in most patients withliver-limited metastatic colorectal cancer. Genome Biol. 16, 32 (2015).
38. Gonzalez-Perez, A. et al. IntOGen-mutations identifies cancer drivers acrosstumor types. Nat. Methods 10, 1081–1082 (2013).
39. The Cancer Genome Atlas Network Comprehensive molecular characterizationof human colon and rectal cancer. Nature 487, 330–337 (2012).
40. Martincorena, I. et al. Universal patterns of selection in cancer and somatic
tissues. Cell 171, 1029–1041 (2017).
41. Mamlouk, S. et al. DNA copy number changes define spatial patterns ofheterogeneity in colorectal cancer. Nat. Commun. 8, 14093 (2017).
42. Yano, J. M. et al. Indigenous bacteria from the gut microbiota regulate hostserotonin biosynthesis. Cell 161, 264–276 (2015).
43. Hayakawa, Y. et al. Nerve growth factor promotes gastric tumorigenesis
through aberrant cholinergic signaling. Cancer Cell 31, 21–34 (2017).
44. Weir, B. S. & Cockerham, C. C. Estimating F-statistics for the analysis of
population-structure. Evolution 38, 1358–1370 (1984).
45. Fitch, W. M. Toward defining course of evolution: minimum change for a
specific tree topology. Syst. Zool. 20, 406 (1971).
46. Naxerova, K. et al. Origins of lymphatic and distant metastases in humancolorectal cancer. Science 357, 55–60 (2017).
47. Beaumont, M. A., Zhang, W. & Balding, D. J. Approximate Bayesiancomputation in population genetics. Genetics 162, 2025–2035 (2002).
48. Marjoram, P. & Tavaré, S. Modern computational approaches for analysing
molecular genetic variation data. Nat. Rev. Genet. 7, 759–770 (2006).
49. Sottoriva, A., Spiteri, I., Shibata, D., Curtis, C. & Tavare, S. Single-molecule
genomic data delineate patient-specific tumor profiles and cancer stem cellorganization. Cancer Res. 73, 41–49 (2013).

1122

Nature Genetics
50. Fumagalli, A. et al. Genetic dissection of colorectal cancer progression byorthotopic transplantation of engineered cancer organoids. Proc. Natl Acad.
Sci. USA 114, E2357–E2364 (2017).
51. Boutin, A. T. et al. Oncogenic Kras drives invasion and maintains metastases
in colorectal cancer. Genes Dev. 31, 370–382 (2017).
52. Wang, Z. et al. Mutational analysis of the tyrosine phosphatome in colorectal
cancers. Science 304, 1164–1166 (2004).
53. Zhang, X. et al. Identification of STAT3 as a substrate of receptor protein
tyrosine phosphatase T. Proc. Natl Acad. Sci. USA 104, 4060–4064 (2007).
54. Lui, V. W. et al. Frequent mutation of receptor protein tyrosine phosphatases
provides a mechanism for STAT3 hyperactivation in head and neck cancer.
Proc. Natl Acad. Sci. USA 111, 1114–1119 (2014).
55. Turajlic, S. et al. Tracking cancer evolution reveals constrained routes tometastases: TRACERx Renal. Cell 173, 581–594 (2018).
56. Rogers, Z. N. et al. Mapping the in vivo fitness landscape of lungadenocarcinoma tumor suppression in mice. Nat. Genet. 50, 483–486 (2018).
57. Cohen, J. D. et al. Detection and localization of surgically resectable cancers
with a multi-analyte blood test. Science 359, 926–930 (2018).
58. Tie, J. et al. Circulating tumor DNA analysis detects minimal residual disease
and predicts recurrence in patients with stage II colon cancer. Sci. Transl.
Med. 8, 346ra92 (2016).
59. Casadaban, L. et al. Adjuvant chemotherapy is associated with improved
survival in patients with stage II colon cancer. Cancer 122, 3277–3287 (2016).

Acknowledgements

C. Curtis is supported by the National Institutes of Health through the NIH Director’sPioneer Award (DP1-CA238296) and NCI Cancer Target Discovery and Development
Network (CA217851). This work was funded in part by grants from the American
Cancer Society (IRG–58-007-54), the Emerson Collective Cancer Research Fund anda gift from the Wunderglo Foundation to C. Curtis. Z.H. is supported by an Innovative
Genomics Initiative (IGI) Postdoctoral Fellowship. The project was supported in partby Cancer Center Support Grants from the National Cancer Institute to the Stanford
Cancer Institute (P30CA124435) and the University of Southern California Norris
Comprehensive Cancer Center (P30CA014089). We thank J. Caswell-Jin and A. Harpak
for critical feedback on the manuscript. This study is dedicated to the memory of
G. Borges, a tireless cancer warrior.

Author contributions

Z.H. implemented the computational and mathematical models, performed simulation
studies and statistical analyses. J.D. implemented the genomic data analysis pipeline,
analyzed and visualized genomic data and provided statistical advice. Z.M. processedclinical samples and generated the genomic data. Z.H., R.S. and J.A.S. analyzed the
genomic data. Z.H., J.D., R.S., J.A.S. and C. Curtis interpreted the data. J.S.S. contributed
to simulation studies. C.J.S., A.S.B. and P.B. performed pathology review. A.S.B. and
M.P. performed immunohistochemistry experiments. M.P., P.B., F.L., C. Cremolini, A.F.
and H.-J.L. contributed clinical samples and expertise. Z.H. and C. Curtis wrote themanuscript, which was reviewed by all authors. C. Curtis conceived and supervised
the study.

Competing interests

A.S.B. has received research support from Daiichi Sankyo and honoraria for lectures,consultation or advisory board participation from Roche Bristol-Myers Squibb, Merckand Daiichi Sankyo as well as travel support from Roche, Amgen and AbbVie. M.P. hasreceived honoraria for lectures, consultation or advisory board participation from thefollowing for-profit companies: Bayer, Bristol-Myers Squibb, Novartis, Gerson Lehrman
Group, CMC Contrast, GlaxoSmithKline, Mundipharma, Roche, Astra Zeneca, AbbVie,
Lilly, Medahead, Daiichi Sankyo and Merck Sharp & Dome. P.B. has received travel
support, honoraria for lectures, consultation or advisory board participation from thefollowing for-profit companies: Biocartis, Novartis, Pfizer, Roche and Roche Diagnostics.
C. Curtis is a scientific advisor to GRAIL and reports stock options as well as consulting
for GRAIL and Genentech.

Additional information

Supplementary information is available for this paper at https://doi.org/10.1038/
s41588-019-0423-x.
Reprints and permissions information is available at www.nature.com/reprints.
Correspondence and requests for materials should be addressed to C.C.
Publisher’s note: Springer Nature remains neutral with regard to jurisdictional claims inpublished maps and institutional affiliations.
© The Author(s), under exclusive licence to Springer Nature America, Inc. 2019

Nature Genetics | VOL 51 | JULY 2019 | 1113–1122 | www.nature.com/naturegenetics

Articles

Nature Genetics
Methods

Clinical specimens, pathology review and sequencing studies. In brief, archived
formalin-fixed paraffin-embedded (FFPE) tissue specimens from 10 patients withmetastatic CRC, including primary tumor, matched metastases and adjacent normalcolon tissue, were obtained from the Medical University of Vienna brain metastasis
biobank, which was established in accordance with ethical guidelines (approval
078/2004). Tissue specimens were collected during the course of routine clinicalcare and clinical data were retrieved by retrospective chart review. All samples
were de-identified and patients in the brain metastasis cohort were deceased prior
to initiating this study. Brain metastases were available for all patients (n = 10)
and for several patients metastases to the liver (n = 1), lung (n = 1) and regional
lymph nodes (n = 4) were also available (Supplementary Table 1). For 6 of the 10
patients, multiple specimens (n = 3–5) from both the primary tumor and metastasis
were sampled and sequenced (Supplementary Table 1). Histological sections wereindependently reviewed by expert pathologists (A.S.B., P.B. and C.J.S.). The Ki-67
proliferative index was determined using immunohistochemical staining using
the Ki-67 antibody, as previously described60. Consistent with the growth of CRC
brain metastases in an expansive rather than infiltrating fashion61, no normal brain
parenchyma was observed within the main brain metastasis lesion.
For all patients, regions of high-cellularity (>60%) were selected for DNA
isolation using the QIAamp DNA FFPE Tissue Kit (Qiagen). Libraries wereprepared using the Agilent SureSelect Human All Exon kit or Ilumina Nextera
Rapid Capture Exome kit for sequencing on the Illumina Hiseq 2000/2500 or
Nextseq 500. Paired sequencing reads were aligned to the human reference genome
build hg19 with BWA (v.0.7.10)62. Duplicate reads were flagged with Picard
Tools (v.1.111). Aligned reads were further processed with GATK 3.4.0 for local
re-alignment around insertions and deletions and base quality recalibration.
We also analyzed de-identified exome-sequencing data from patients with
mCRC in four published datasets21,29–31 using the same unified bioinformatics
framework described below. After excluding tumors with low purity (<0.4), we
retained 46 tumor specimens from 13 patients with mCRC who had paired samplesof liver metastases and refer to this as the liver metastasis cohort.
Somatic SNV detection and filtering. sSNVs were called by MuTect (v.1.1.7)63
with paired tumor and normal sequencing data. sSNVs that failed to pass theinternal filters of MuTect, had fewer than 10 total reads or 3 variant reads in thetumor sample, fewer than 10 reads in the normal sample or mapped to paralogous
genomic regions were removed. Additional Varscan (v.2.3.9)64 filters were applied
to remove sSNVs with low average variant base qualities, low average mapping
qualities among variant-supporting reads, strand bias among variant-supporting
reads and high average mismatch base quality sums among variant-supporting
reads, either within each tumor sample or across all tumor samples from thesame patient. Additional filtering removed sSNVs detected in a panel of normals
by running MuTect in single-sample mode with less stringent filtering criteria(artefact detection mode). sSNVs called in at least two normal samples were
included in the panel of normal sSNV list. For FFPE samples, sSNVs called insamples from one patient were checked against samples from all other patientsto flag those that might be artefactual. The maximum observed VAFs across all
samples from each patient were calculated based on raw output files from
MuTect. sSNVs with maximum observed VAFs between 0.01 and 0.05 in at least
two other patients were removed. Small indels were called with Strelka (v.1.0.14)
and annotated by Annovar (v.20150617)65. sSNVs and small indels in proteincoding regions were retained for downstream analyses. Additional filters wereapplied to exclude possible artefactual sSNVs due to the processing of FFPE
specimens. Specifically, artefacts among C>T/G>A sSNVs with bias in read pairorientation were filtered in each individual FFPE sample, similar to the previously
described approach66.
For patients with MRS data, we sought to exploit this information by retrieving
read counts for sSNVs across samples from the same patient. To obtain depth and
VAF information across all samples from the same patient, for each sSNV and ineach tumor sample that an sSNV was not originally called in, the total reads andvariant-supporting reads were counted using the mpileup command in SAMtools(v.1.2)67. Only reads with mapping quality ≥40 and base quality at the sSNV locus≥20 were counted and used to calculate VAF values for that sSNV.
Copy-number analysis, tumor purity and CCF estimation. Copy-number
analysis was performed using TitanCNA (v.1.5.7)68. In brief, TitanCNA uses depth
ratio and B-allele frequency information to estimate allele-specific absolute copy
numbers with a hidden Markov model, and estimates tumor purity and clonalfrequencies. Only autosomes were used in copy-number analysis. First, for eachpatient, germline heterozygous SNPs based on dbSNP build 138 were identified
using SAMtools and SnpEff (v.3.6) in the normal sample. HMMcopy (v.0.99.0)69
was used to generate read counts for 1,000-bp bins across the genome for all tumor
samples. Whole-exome sequences from multiple normal samples per patientwere pooled separately for the purpose of calculating read counts in the bins andthe pooled normal read depth data were used as controls only for the calculationof depth ratios. TitanCNA was used to calculate allelic ratios at the germline
heterozygous SNP loci in the tumor sample and depth ratios between the tumor
sample and the pooled normal data in bins that contained those SNP loci. Only
Nature Genetics | www.nature.com/naturegenetics

SNP loci within whole-exome sequence-covered regions were then used to estimate
allele-specific absolute copy-number profiles. TitanCNA was run with different
numbers of clones (n = 1–3). One run was chosen for each tumor sample based on
visual inspection of fitted results, with preference given to the results with a single
clone unless results with multiple clones had visibly better fit to the data. Resultsfrom tumor samples from the same patient were inspected together to ensure
consistency. Overall ploidy and purity for each tumor sample was calculated
from the TitanCNA results. For the public datasets including liver-exclusive
mCRCs, cases with estimated purity >0.4 in both the primary tumor and paired
metastases (Supplementary Fig. 2) were included since low purity hinders accurate
SNV/CNA calling.
Mutational CCFs were estimated with CHAT (v.1.0)70. CHAT includes afunction to estimate the CCF of each sSNV by adjusting its VAF based on local
allele-specific copy numbers at the sSNV locus. sSNV frequencies and copynumber profiles estimated from previous steps were used to calculate CCFs for allsSNVs in autosomes (using a modified function). The CCFs were also adjusted
for tumor purity. The merged CCF of each sSNV is computed by integrating CCFs
from multiple regions when MRS data are available: ∑k CCF ×
 i =1
i dik
CCF = ∑i =1 di1
CCF < 1

(1)

CCF ≥ 1

where di and CCFi are the sequencing depth and CCF estimation in region i,
respectively. Of note, the vast majority (99%) of P/M shared sSNVs have CCF (ormerged CCF) >60%, a cut-off that also optimally distinguishes the site-privateclonal and subclonal sSNV clusters (Supplementary Fig. 6). We thus use 60% as the
CCF cut-off to define clonal versus subclonal sSNVs in the PMGD analysis.
Data processing for downstream analyses. For each tumor site (primary or
metastasis) in a patient, the average CCF estimate of a sSNV is set to 0 if neitherof the following two criteria are met (1) VAF ≥ 0.03 and variant read count ≥3; (2)VAF ≥ 0.1 in any of the regions. The following additional filters were applied tosummarize the MRS P/M data in a given patient. First, filter out sSNVs without
VAF ≥ 0.05 and variant read count ≥3 or VAF ≥ 0.1 in any samples from this pairof sites. Second, filter out sSNVs with total read depth <20 from either of the twotumor sites. Third, filter out all sSNVs in chromosome regions with LOH in all
specimens from one tumor site but not in all samples from the other tumor site.
Fourth, for sSNVs not present in any specimens with LOH, filter out sSNVs thatsatisfy the following criteria in specimens from at least one of the two tumor sites:(1) absent in some samples with LOH; (2) present in any samples without LOH.
Driver enrichment analysis. Driver fold enrichment was determined based oncolorectal adenocarcinoma driver genes (defined by combining IntOGen v.2016.538
and TCGA39 including 221 genes, Supplementary Table 3) or all pan-cancer
drivers, including 369 high-confidence genes40 that had non-silent coding sSNVs/
indels out of the total number of genes with non-silent coding sSNVs/indels. The
resulting metric was normalized to the fraction of driver genes out of all genes in
the human genome. Clonal mutations (CCF > 60% in primary tumor or metastasis;merged CCF was used for MRS data) were divided into three sets that represented
shared, primary tumor-private and metastasis-private mutations, for which only
distant metastases were considered. Driver gene fold enrichment was calculated foreach set of mutations by randomly sampling 15 out of 25 P/M pairs from the wholecohort, aggregating them to calculate one driver enrichment score, and repeating
this analysis 100 times (n = 100 downsamplings) to derive a test statistic. For eachdownsampling, the driver enrichment score was calculated as:
Enrichment fold scoren(driver non − silent clonal)∕n(all non − silent clonal)=n(driver genes)∕n(total genes)(2)

where n(all non-silent clonal) and n(driver non-silent clonal) correspond to thetotal number of non-silent clonal mutations and the number of non-silent clonal
mutations in driver genes, respectively. Here n(driver genes) and n(total genes)
correspond to the total number of drivers reported for CRC (n = 221) or pancancer (n = 369) and the number of coding genes in the genome (n = 22,000),
respectively.
Orthogonal validation of early metastasis driver gene modules. Clinicalannotations and targeted sequencing data were obtained for the GENIE33 (v.3.0)
and MSK-Impact32 CRC cohorts from Synapse (http://synapse.org/genie) andcBioPortal (http://www.cbioportal.org/study?id=crc_msk_2018), respectively. The
MSK-Impact cohort includes early-stage primary CRCs, primary CRCs that are
known to have metastasized and the metastatic lesion (predominantly liver) from
1,099 patients with mCRC and a total of 1,134 samples with available sequencing
and clinical covariates, including stage, microsatellite status and time to metastasis.
As the mCRC discovery cohort did not include any cases with microsatelliteArticles
unstable tumors, these were removed, as were cases with POLE mutations.
Microsatellite stable samples were divided into early-stage non-metastatic samples(n = 57), metastatic primary tumors (n = 440) and metastatic samples (n = 498).
The GENIE cohort is composed of 39,600 samples profiled with different
targeted sequencing panels from which CRC samples were selected (oncotree
codes: COADREAD, COAD, CAIS, MACR, READ and SRCCR). In order to avoid
duplicated samples, all MSK-Impact samples from the GENIE cohort were removed,
as were duplicated samples from the same patient, resulting in 2,666 samples, 1,756
of which were from primary tumors. As the GENIE cohort does not currently
include stage or outcome information, all primaries are assumed to be nonmetastatic, although some may be stage IV or diagnosed as metastatic in the future.
All possible combinations of recurrent putative metastasis driver genes (APC,
TP53, KRAS, SMAD4, PIK3R1, BRAF, AMER1, TCF7L2, PIK3CA, PTPRT andATM) identified in the mCRC cohort were evaluated in metastatic relative toearly-stage cases using a two-sided Fisher’s exact test (with Benjamini–Hochberg
adjustment for multiple testing). The enrichment analysis was calculated forthe combined MSK-Impact and GENIE primary CRC cohort, as well as for the
MSK-Impact cohort alone (Supplementary Table 8). As the number of genes ina module increases, the specificity of the association with metastasis increases,
whereas the frequency of the module and in turn power to detect an association
decreases (Supplementary Fig. 25). Although combining datasets may potentiallyintroduce some biases, because we assume that all GENIE primary samples are
non-metastatic and microsatellite stable, this will render our analyses conservative.
Indeed, it is worth noting that although these results are already highly significant,they are likely conservative for several reasons: (1) due to the short follow-up
time, some cases with early-stage tumors may develop metastases; (2) imbalanced
sample size with nearly twice as many patients with early-stage disease versus casesof metastatic disease; (3) several putative metastasis drivers that were identified in
the mCRC cohort are not represented on the targeted sequencing panel and hencecannot be evaluated.
Phylogenetic tree reconstruction and FST computation. We ran PHYLIP71(http://www.trex.uqam.ca/index.php?action=phylip&app=dnapars) and applied
the maximum parsimony method to reconstruct the phylogeny of multiple
specimens from individual patients based on the presence or absence of SNVs
and indels. When multiple maximum parsimony trees were reported, we chosethe top ranked solution. FigTree (http://tree.bio.ed.ac.uk/software/Figtree/) wasused to visualize the reconstructed trees. We computed the FST statistic for eachprimary tumor or metastasis using the Weir and Cockerham method44 based on
the adjusted frequency of subclonal sSNVs (merged CCF < 60%) identified in MRS
data. Clonal mutations (merged CCF > 60%) did not contribute to intratumor
heterogeneity and were excluded in FST calculations.
Spatial agent-based modeling of tumor progression. We extended our previously
described three-dimensional agent-based tumor evolution framework19,23 tomodel tumor growth, mutation accumulation and metastatic dissemination
after malignant transformation under different evolutionary scenarios in P/Mpairs. Pre-malignant clonal expansions before transformation do not alter the
genetic heterogeneity within a tumor and were therefore were not modeled(Figs. 1c, 4a and Supplementary Fig. 15). We assume that dissemination occurs
after malignant transformation of the founding carcinoma cell as invasion(a cardinal feature of carcinomas) is a requirement for metastasis. We havepreviously used this framework to model primary tumor evolution23. In thismodel, spatial tumor growth is simulated by the expansion of deme subpopulations(composed of approximately 5,000 cells with diploid genome), mimicking
the glandular structures that are often found in CRCs and metastases andconsistent with the number of cells found in individual CRC glands (around2,000–10,000 cells)72. Model assumptions are detailed in Supplementary Table 5.
The deme subpopulations expand within a defined three-dimensional cubic
lattice (Moore neighborhood, 26 neighbors), through peripheral growth whilecells within each deme are well-mixed without spatial constraints and grow bya random birth-and-death process (division probability P and death probability
Q = 1 − P at each generation). The notion of peripheral growth is supported byrecent studies, which indicated that cancer cells at the periphery of the tumor
proliferate much faster than those at the center73. Moreover, peripheral growth
results in a power law model of net tumor growth (Supplementary Fig. 15b) andis supported by data on CRCs74. The first deme is generated via the same birthand-death process, beginning with a single transformed founding tumor cell. Here
we used the following parameters: P = 0.55 and Q = 0.45 for the deme expansion
in both the primary tumor and metastasis. Thus the cell birth/death probability
ratio for the founding lineage is P/Q = 0.55/0.45 ≈ 1.2. This is supported by theobservation that there is no significant difference in proliferation rates based on
Ki-67 staining of paired CRCs and brain metastases (Supplementary Fig. 10b),
as previously reported in liver metastases75. Based on these values of P and Q,
approximately 3 years are required from transformation to the diagnosis of primary
carcinoma (approximately 109 cells; Supplementary Fig. 15b). Once a deme exceeds
the maximum size (10,000 cells), it splits into two offspring demes via random
sampling of cells from a binomial distribution (Nc, P = 0.5), where Nc is thecurrent deme size.

Nature Genetics
During the growth of the primary CRC, a single cell from a random deme at
the tumor periphery is randomly chosen to seed the metastasis, which supported
by mounting pathological evidence of invasive cells in the tumor front and thefact that blood vessels are also mostly distributed in the invasive front in CRC76.
The total cell number at the time of metastatic dissemination is denoted by Nd.
The metastasis grows via the same model as the primary tumor, starting from thedisseminated tumor cell(s).
During each cell division, the number of neutral passenger mutations acquiredin the coding portion of the genome follows a Poisson distribution with mean u.
Thus, the probability that k mutations occurred in each cell division is as follows:
P (x = k ) =u ke−uk!

(3)

where an infinite sites model and constant mutation rate are assumed during tumor
progression. For simplicity, we do not simulate CNAs, LOH or aneuploidy, andall mutations are heterozygous. Under the neutral model, all somatic mutations
are assumed to be neutral passenger events and do not confer a fitness advantage,
whereas in the subclonal selection model, beneficial mutations (or advantageous
mutations) arise stochastically via a Poisson process with mean us during each cell
division. We assume us = 10−5 per cell division in the genome23,77. We investigated a
relatively strong positive selection coefficient (s = 0.1), where s specifies the increasein cell division probability per cell division when a beneficial mutation occurs
in the neutral cell lineage. The cell birth and death probabilities for a selectively
beneficial clone are Ps = P(1 + s) and Qs = 1 − Ps = 1 − P(1 + s), respectively, thus
the selective advantage is defined as s = Ps/P − 1. We selected s = 0.1, as we havepreviously shown that the resultant patterns of between-region genetic divergencecan be clearly distinguished from those arising under effectively neutral growth23.
During simulation of primary and metastatic growth, each mutation is
assigned a unique index that is recorded with respect to its genealogy and hostcells, enabling analysis of the mutational frequency in a sample of tumor cells or
the whole tumor during different stages of growth. We simulate growth until the
primary and metastasis reach a size of approximately 109 cells (or around 10 cm3)
comparable to the size of the clinical samples studied here, which ranged from4 to 15 cm in maximum diameter. To simulate each of the four scenarios of P/M
growth, namely N/N, N/S, S/N or S/S, we used a mutation rate u = 0.3 per cell
division in the exonic region (corresponding to 5 × 10−9 per site per cell division in
the 60 Mb diploid coding regions) and selection coefficients s = 0 and s = 0.1 whenmodeling neutral evolution and subclonal selection, respectively, during growth ofthe primary tumor or metastasis. Under each of the four scenarios of P/M growth,
100 time points (representing the primary tumor size at the time of dissemination,Nd) were sampled at random from a uniform distribution, log10(Nd) ~ U(2, 9), eachgiving rise to independent P/M pairs. The CCF from the whole tumor in boththe primary tumor and metastatic lesions were obtained for each sSNV (site).
CCFs > 60% in one site and CCFs < 1% in the other site were used to count the
number of primary tumor-private and metastasis-private clonal sSNVs (Lp and Lm,
respectively), consistent with the strategy used for patient samples.
SCIMET. We sought to infer two parameters that govern the dynamics ofmetastasis, namely u, the mutation rate per cell division in the exonic region, andNd, the primary tumor size at the time of dissemination based on our spatial tumor
simulation framework. The two parameters of interest (u and Nd) were randomly
sampled from a prior discrete uniform distribution, namely 10 values from 0.003
to 3 for u; and 7 values from 103 to 109 cells (on log10 scale) for Nd (Supplementary
Tables 6, 7 and Supplementary Fig. 19). Discrete prior distributions for u and Nd
were used to estimate the order of magnitude rather than the precise values ofthese two parameters. We simulated 70,000 paired primary tumors and metastases(composed of 109 cells each) under each of the four evolutionary scenarios (N/N,
N/S, S/N or S/S). After generating the virtual P/M tumors, multiple regions (n = 4),each composed of approximately 106 cells, are sampled from an octant of tumor
sphere, as was done for our clinical samples (Supplementary Fig. 19). The VAF of
all sSNVs in the sampled bulk subpopulation is considered the true VAF (denotedby fT), whereas the observed allele frequency is obtained via a statistical model thatmimics the random sampling of alleles during sequencing. Specifically, we used abinomial distribution (n, fT) to generate the observed VAF at each site given its true
frequency fT and number of covered reads n. The number of covered reads at eachsite is assumed to follow a negative-binomial distribution (negative binomial(size,
depth)) where depth is the mean sequencing depth and size corresponds to thevariation parameter78. We assume depth = 80 and size = 2 for the sequencing datain each tumor region. A mutation is called when the number of variant reads is≥3, thereby applying the same criteria as for the patient tumors. The observed
VAF for each mutation is converted to CCF and the merged CCF from four
regions were computed (equation (1)) to mimic the patient genomic data. The
nine summary statistics used to fit the CCF data are described in Supplementary
Fig. 19 and Supplementary Table 6. The median values of the posterior probability
distributions obtained from SCIMET are referred to as the inferred parametervalues (ũ and N͠ d). To be conservative, we define early dissemination as Nd (upperbound) < 108 cells (around 1 cm3 in volume) using the third quartile of the posterior
distribution as the upper bound, whereas late dissemination is defined as Nd (upper
Nature Genetics | www.nature.com/naturegenetics

Articles

Nature Genetics
bound) ≥ 108 cells (Fig. 5a). We also evaluated the robustness of SCIMET to ahigher birth/death rate ratio (Supplementary Fig. 21), collective dissemination bya cell cluster (n = 10 cells; Supplementary Fig. 22) or single-region sequencing data
(Supplementary Fig. 23). Of note, both a higher birth/death rate ratio and singleregion sequencing data would result in overestimation of the timing of metastatic
dissemination. A higher birth/death rate ratio yields a higher tumor growth rate
thus the inferred primary tumor size at the time of dissemination would be larger
than for a lower birth/death rate ratio. Single-region sampling results in a larger
number of metastasis-private clonal mutations (larger Lm and larger H) comparedwith MRS, thus the timing of dissemination would be overestimated in accordancewith the positive correlation between Lm or H and Nd. Overall, these comparisonsdemonstrate the robustness of SCIMET to different model assumptions.
We used a version of ABC based on the acceptance–rejection algorithm79 toestimate posterior probability distributions for the parameters of interest θ(u, Nd).
The ABC version of rejection sampling is as follows:
For i = 1 to K under model M (N/N, N/S, S/N or S/S):
1.
2.
3.
4.

Sample parameters θ′ from the prior distribution π(θ).
Simulate data D′ using model M with the sampled parameters θ′ and summarize D′ as summary statistics S′.
Accept θ′ if d(S′, S) < ε, for a given tolerance rate ε, where d(S′, S) is a measure of Euclidean distance between S′ and S.
Go to (1).

Using this scheme, we are able to approximate the posterior distribution by:P(θ|d(S′, S) < ε). We use a common variation of ABC47,80 in which, rather thanusing a fixed threshold ε, we sort all K distances calculated by d(S′, S) (step (3)),
and accept the θ′ that generated the smallest 100 × η percentage distances. We
use η = 0.01 so that the posterior is composed of 70,000 × 0.01 = 700 data points.
The ABC procedure is performed using the R package abc81. To determine themost probable model of tumor evolution (N/N, N/S, S/N or S/S) in P/M pairs, we
run the postpr method implemented in the R package abc, which integrates all
simulation data from the four models to run the ABC procedures (steps (1)–(4))
and outputs the probability of each model based on the posterior distribution. The
model (N/N, N/S, S/N or S/S) with the highest probability was selected.
A Monte Carlo cross-validation scheme was performed to assess theperformance of SCIMET. This procedure involves randomly sampling a
combination of parameters u′ and Nd′ (true parameters) and sampling
10 simulations of the summary statistics S′ under this parameter set toindependently run the ABC scheme. The posterior parameters u and Nd with the
maximum probability were used as parameter estimates for one simulation,
namely ũ and N͠ d. The mean value of ũ and N͠ d in 10 simulations was taken as theparameter estimate (inferred parameters) in the cross-validation. The process ofMonte Carlo sampling and SCIMET inference was repeated 200 times under eachof the four evolutionary scenarios (N/N, N/S, S/N and S/S). Comparison of theinferred versus true parameter values indicates the robustness of this approach(Supplementary Fig. 20).
Reporting Summary. Further information on research design is available in the
Nature Research Reporting Summary linked to this article.

Data availability

Data have been deposited at the European Genotype Phenotype Archive (EGA)
under accession number EGAS00001003573. Data from previously publishedstudies are available from the DDBJ (accession number JGAS00000000060)21 andthe SRA (accession numbers SRP052609, SRP074289 and SRP041725)29–31.

Nature Genetics | www.nature.com/naturegenetics

Code availability

Code used for genomic data analysis and simulation studies are available at https://github.com/cancersysbio/mCRCs and https://github.com/cancersysbio/SCIMET.

References

60. Berghoff, A. S. et al. Differential role of angiogenesis and tumour cellproliferation in brain metastases according to primary tumour type: analysis
of 639 cases. Neuropathol. Appl. Neurobiol. 41, e41–e55 (2015).
61. Berghoff, A. S. et al. Invasion patterns in brain metastases of solid cancers.
Neuro-oncol. 15, 1664–1672 (2013).
62. Li, H. & Durbin, R. Fast and accurate short read alignment with Burrows–
Wheeler transform. Bioinformatics 25, 1754–1760 (2009).
63. Cibulskis, K. et al. Sensitive detection of somatic point mutations in impure
and heterogeneous cancer samples. Nat. Biotechnol. 31, 213–219 (2013).
64. Koboldt, D. C. et al. VarScan 2: somatic mutation and copy number alteration
discovery in cancer by exome sequencing. Genome Res. 22, 568–576 (2012).
65. Wang, K., Li, M. & Hakonarson, H. ANNOVAR: functional annotation ofgenetic variants from high-throughput sequencing data. Nucleic Acids Res. 38,
e164 (2010).
66. Costello, M. et al. Discovery and characterization of artifactual mutations indeep coverage targeted capture sequencing data due to oxidative DNA
damage during sample preparation. Nucleic Acids Res. 41, e67 (2013).
67. Li, H. et al. The Sequence Alignment/Map format and SAMtools.
Bioinformatics 25, 2078–2079 (2009).
68. Ha, G. et al. TITAN: inference of copy number architectures in clonal cell
populations from tumor whole-genome sequence data. Genome Res. 24,
1881–1893 (2014).
69. Ha, G. et al. Integrative analysis of genome-wide loss of heterozygosity andmono-allelic expression at nucleotide resolution reveals disrupted pathways intriple-negative breast cancer. Genome Res. 22, 1995–2007 (2012).
70. Li, B. & Li, J. Z. A general framework for analyzing tumor subclonality using
SNP array and DNA sequencing data. Genome Biol. 15, 473 (2014).
71. Felsenstein, J. Phylogeny inference package. Cladistics 5, 164–166 (1989).
72. Siegmund, K. D., Marjoram, P., Woo, Y. J., Tavaré, S. & Shibata, D. Inferringclonal expansion and cancer stem cell dynamics from DNA methylation
patterns in colorectal cancers. Proc. Natl Acad. Sci. USA 106, 4828–4833 (2009).
73. Lloyd, M. C. et al. Darwinian dynamics of intratumoral heterogeneity: notsolely random mutations but also variable environmental selection forces.
Cancer Res. 76, 3136–3144 (2016).
74. Sarapata, E. A. & de Pillis, L. G. A comparison and catalog of intrinsic tumor
growth models. Bull. Math. Biol. 76, 2010–2024 (2014).
75. Finlay, I. G., Meek, D., Brunton, F. & McArdle, C. S. Growth rate of hepatic
metastases in colorectal carcinoma. Br. J. Surg. 75, 641–644 (1988).
76. Kather, J. N. et al. Identification of a characteristic vascular belt zone inhuman colorectal cancer. PLoS One 12, e0171378 (2017).
77. Bozic, I. et al. Accumulation of driver and passenger mutations during tumor
progression. Proc. Natl Acad. Sci. USA 107, 18545–18550 (2010).
78. Love, M. I., Huber, W. & Anders, S. Moderated estimation of fold change anddispersion for RNA-seq data with DESeq2. Genome Biol. 15, 550 (2014).
79. Tavaré, S., Balding, D. J., Griffiths, R. C. & Donnelly, P. Inferring coalescencetimes from DNA sequence data. Genetics 145, 505–518 (1997).
80. Zhao, J., Siegmund, K. D., Shibata, D. & Marjoram, P. Ancestral inference intumors: how much can we know? J. Theor. Biol. 359, 136–145 (2014).
81. Csilléry, K., François, O. & Blum, M. G. abc: an R package for approximate
Bayesian computation (ABC). Methods Ecol. Evol. 3, 475–479 (2012).

Reporting Summary
Nature Research wishes to improve the reproducibility of the work that we publish. This form provides structure for consistency and transparency
in reporting. For further information on Nature Research policies, see Authors & Referees and the Editorial Policy Checklist.

Statistical parameters

nature research | reporting summary

Corresponding author(s): Christina Curtis

When statistical analyses are reported, confirm that the following items are present in the relevant location (e.g. figure legend, table legend, main
text, or Methods section).
n/a Confirmed
The exact sample size (n) for each experimental group/condition, given as a discrete number and unit of measurement
An indication of whether measurements were taken from distinct samples or whether the same sample was measured repeatedly
The statistical test(s) used AND whether they are one- or two-sided
Only common tests should be described solely by name; describe more complex techniques in the Methods section.

A description of all covariates tested
A description of any assumptions or corrections, such as tests of normality and adjustment for multiple comparisons
A full description of the statistics including central tendency (e.g. means) or other basic estimates (e.g. regression coefficient) AND
variation (e.g. standard deviation) or associated estimates of uncertainty (e.g. confidence intervals)
For null hypothesis testing, the test statistic (e.g. F, t, r) with confidence intervals, effect sizes, degrees of freedom and P value noted
Give P values as exact values whenever suitable.

For Bayesian analysis, information on the choice of priors and Markov chain Monte Carlo settings
For hierarchical and complex designs, identification of the appropriate level for tests and full reporting of outcomes
Estimates of effect sizes (e.g. Cohen's d, Pearson's r), indicating how they were calculated
Clearly defined error bars
State explicitly what error bars represent (e.g. SD, SE, CI)
Our web collection on statistics for biologists may be useful.

Software and code
Policy information about availability of computer code
Data collection

BWA-MEM (0.7.10), Picard (1.111), GATK (3.4.0), MuTect (1.1.7), Varscan (2.3.9), Strelka (1.0.14), Annovar (20150617), DToxoG(1.14.4.0), TitanCNA (1.5.7), SAMtools (1.2), SnpEff (3.6), HMMcopy (0.99.0), CHAT (1.0), GENIE (3.0), cBioPortal (http://
www.cbioportal.org/study?id=crc_msk_2018), IntOGen (2016.5)

Data analysis

R version 3.5.0, PHYLIP (http://www.trex.uqam.ca/index.php?action=phylip&app=dnapars), FigureTree (http://tree.bio.ed.ac.uk/
software/Figuretree/), Custom software: mCRCs (https://github.com/cancersysbio/mCRCs), SCIMET (https://github.com/cancersysbio/
SCIMET)For manuscripts utilizing custom algorithms or software that are central to the research but not yet described in published literature, software must be made available to editors/reviewers
upon request. We strongly encourage code deposition in a community repository (e.g. GitHub). See the Nature Research guidelines for submitting code & software for further information.

April 2018

1

Policy information about availability of data
All manuscripts must include a data availability statement. This statement should provide the following information, where applicable:- Accession codes, unique identifiers, or web links for publicly available datasets- A list of figures that have associated raw data
- A description of any restrictions on data availability
Data have been deposited at the European Genotype Phenotype Archive (EGA) under accession number EGAS00001003573. Data from previously published studies
are available at: DDBJ: JGAS00000000060 (Uchi et al.), and the SRA: SRP052609 (Kim et al.), SRP074289 (Leung et al.), SRP041725 (Lim et al.).

Field-specific reporting
Please select the best fit for your research. If you are not sure, read the appropriate sections before making your selection.
Life sciences

Behavioural & social sciences

nature research | reporting summary

Data

Ecological, evolutionary & environmental sciences

For a reference copy of the document with all sections, see nature.com/authors/policies/ReportingSummary-flat.pdf

Life sciences study design
All studies must disclose on these points even when the disclosure is negative.
Sample size

Statistical methods were not used to predetermine sample size - rather all available samples that were suitable for inclusion were utilized.
Exome sequencing data from 118 biopsies from 23 metastatic colorectal cancer (mCRC) patients with paired liver and brain metastases wereanalyzed. Additionally, targeted sequencing data from an independent cohort of 2,751 colorectal cancer patients was utilized to corroborate
these findings.

Data exclusions

Samples with low tumor purity (<0.4) from publicly available metastatic colorectal cancer datasets were excluded since this adversely affects
mutation detection, as noted in the Methods (page 12) and Figure S2.

Replication

Replication is not applicable for patient genomic data. For the simulation studies, the number of replicates are reported in the text. In
particular, statistical inference of patient-specific parameters using SCIMET requires a large number of simulations, we performed 1000
simulations for each of the parameter combinations (n=70000 in total for each model).

Randomization

This was an observational study, no randomization was performed.

Blinding

Blinding was not considered appropriate for this study.

Reporting for specific materials, systems and methods
Materials & experimental systems

Methods

n/a Involved in the study

n/a Involved in the study

Unique biological materials

ChIP-seq

Antibodies

Flow cytometry

Eukaryotic cell lines

MRI-based neuroimaging

Palaeontology
Animals and other organisms
Human research participants

April 20182