Cytogenetic signatures favoring metastatic organotropism in colorectal cancer

Abstract

Colorectal carcinoma (CRC) exhibits metastatic organotropism, primarily targeting liver, lung, and rarely the brain. Here, we study chromosomal imbalances (CIs) in cohorts of primary CRCs and metastases. Brain metastases show the highest burden of CIs, including aneuploidies and focal CIs, with enrichment of +12p encoding KRAS. Compared to liver and lung metastases, brain metastases present with increased co-occurrence of KRAS mutation and amplification. CRCs with concurrent KRAS mutation and amplification display significant metabolic reprogramming with upregulation of glycolysis, alongside upregulation of cell cycle pathways, including copy number gains of MDM2 and CDK4. Evolutionary modeling suggests early acquisition of many organotropic CIs enriched in both liver and brain metastases, while brain-enriched CIs preferentially emerge later. Collectively, this study supports a model where cytogenetic events in CRCs favor site-specific metastatic colonization. These site-enriched CI patterns may serve as biomarkers for metastatic potential in precision oncology.

Fig. 1. Pattern of CIs observed in the CRCTropism cohort.

a Chromosomal losses and gains in shades of blue and red, respectively, observed in the cohort with 314 tumors, including n = 80 primary CRC and n = 234 distant CRC metastases in the liver, lung, and brain. On the left, tumor site and sex are color-encoded according to the legend. b Number of apparent chromosomal arm aneuploidies in relation to focal CIs, visualized as a stacked bar chart (n = 234). c Distribution of chromosomal arm aneuploidies in CRC metastases (n = 191). (1) p_adjust = 2.96 × 10⁻³; (2) p_adjust = 6.94 × 10⁻⁶; (3) p_adjust = 1.31 × 10⁻². d Distribution of CIs including losses and gains for CRC metastases (n = 191). (1) p_adjust = 2.46 × 10⁻³; (2) p_adjust = 1.37 × 10⁻⁶; (3) p_adjust = 5.57 × 10⁻³; (4) p_adjust = 3.70 × 10⁻⁴; (5) p_adjust = 5.96 × 10⁻³; (6) p_adjust = 1.44 × 10⁻⁴; (7) p_adjust = 4.97 × 10⁻⁷; (8) p_adjust = 1.75 × 10⁻². c, d Data are presented as combined violin and box plots (thick line in the box corresponds to median; lower and upper edges of the box indicating the first (Q1) and third (Q3) quartiles; whiskers extend to 1.5 times the IQR). Individual data points are presented using the jitter method. Kruskal-Wallis tests with Dunn posthoc tests (two-sided) and Benjamini-Hochberg corrections for multiple comparisons were used to test for statistical significance. *, p_adjust < 0.05; **, p_adjust < 0.01, and ***, p_adjust < 0.001. Source data are provided as a Source Data file.

Fig. 2. Site-specific pattern of DNA copy number aberrations in CRC metastases based on the CRCTropism cohort.

a Organotropic map presenting CIs with site-specific enrichment (n = 191). All p-values were obtained by two-sided Fisher exact tests followed by Benjamini-Hochberg adjustment for multiple comparison (legend). The size of the bubbles encodes the frequency a given CI was observed in the cohort (legend). b Doughnut plots showing the fraction of metastases (n = 191) that harbor selected organotropic CIs (see legend for the assignment of secondary organs: li, liver; lu, lung; br, brain; chr., chromosomal arm; see also Supplementary Fig. 5 for a full representation of all organotropic CIs). The frequency of occurrence is color-encoded according to the legend shown on the right. Two-sided Fisher exact tests with Benjamini-Hochberg adjustment for multiple comparison were used (*, p_adjust < 0.05; **, p_adjust < 0.01, and ***, p_adjust < 0.001). (1) p_adjust = 3.80 × 10⁻²; (2) p_adjust = 1.14 × 10⁻²; (3) p_adjust = 2.48 × 10⁻²; (4) p_adjust = 3.80 ×  10⁻²; (5) p_adjust = 4.10 × 10⁻²; (6) p_adjust = 7.80 × 10⁻³; (7) p_adjust = 3.80 × 10⁻²; (8) p_adjust = 4.85 × 10⁻³; (9) p_adjust = 4.22 × 10⁻³; (10) p_adjust = 3.08 × 10⁻⁷. c Graphical summary of autosomal organotropic CIs with arbitrary intra-organ sites of lesions. CIs observed in more than 15% of the respective metastases are highlighted in bold. Source data are provided as a Source Data file.

Fig. 3. Organotropic genetic signatures of metastatic CRC in the MetTropism cohort.

a Pattern of chromosomal arm aneuploidies in liver, lung, and brain metastases determined using ASCETS, with n = 792 metastases subjected to analysis. For each metastatic site, the frequency of the chromosomal arm aneuploidy (bottom, deletion; top, amplification) is color-coded as shown in the color scale at the right. b Profile of tumor mutational burden (TMB) and fraction genome altered (FGA) for each liver, lung, and brain metastasis (FGA: p = 8.23 × 10⁻⁹ [Kruskal-Wallis test], p_adjust = 4.66 × 10⁻⁹ for brain metastases vs. liver and p_adjust = 5.50 × 10⁻⁹ for brain metastases vs. lung metastases [Dunn posthoc test with Benjamini-Hochberg correction]; TMB: p = 0.113 [Kruskal-Wallis test]; n = 792). c Organotropic map of gene-level gains and losses derived from GISTIC 2.0 (n = 792). Significant CIs are indicated in blue, with the event frequency represented by the bubble size (legend). The doughnut plot in the top right shows the fraction of metastases with +12p as determined using ASCETS (n = 758), where the occurrence frequency is color-encoded according to the legend on the right. (1) p_adjust = 5.19 × 10⁻³. d Co-occurrence probability mapping of genetic changes observed in CRC metastases (n = 792), with red shades indicating significant positive associations and blue shades depicting significant negative associations, based on co-occurrence modeling (two-sided), followed by Benjamini-Hochberg adjustment. The inset presents close-ups of selected genetic changes (amp, amplification; del, deletion; mut, mutation). e Fraction of metastases (n = 792) harboring KRAS amplifications determined using GISTIC 2.0 (color-coded according to the legend). (1) p_adjust = 3.49 × 10⁻³. f Bar plot depicting the frequency of KRAS amplifications among metastases (n = 352) with oncogenic KRAS mutations. (1) p_adjust = 0.019; (2) p_adjust = 0.019. c, e, f Adjusted p-values were obtained by two-sided Fisher exact test with Benjamini-Hochberg adjustment for multiple comparisons. *, p_adjust < 0.05; **, p_adjust < 0.01. Source data are provided as a Source Data file.

Fig. 4. Rewiring of KRAS mutated and amplified CRCs as indicated by the TCGA cohort.

a Venn diagram showing the relationship between CRCs with KRAS mutations (blue circle) and amplifications (red circles). b Box plot of normalized KRAS expression levels, with p-values from two-sided Wilcoxon rank-sum test and Benjamini-Hochberg adjustment for multiple comparison (n = 64 and n = 133 for CRC with KRAS mutation and amplification and CRC with KRAS mutation, respectively). (1) p_adjust = 1.32 × 10⁻³. c Multiple linear regression indicating KRAS_CNV as a significant predictor of phosphorylated MAP2K1/2, along with pBRAF_S445, and pMAP2K1 with similar variable importance scores (Score_VI), with n = 160. d. Two-sided GSEA using the M2900 gene set, derived from cells overexpressing oncogenic KRAS. The normalized enrichment score (NES), nominal p-value and false discovery rate (FDR)-adjusted p-value are shown (n = 197). e Volcano plot of differential gene expression, colored by chromosomal arm (see legend). The inset highlights genes on chromosome 12. Statistical significance was assessed using two-sided Wilcoxon rank-sum tests with Benjamini-Hochberg correction for multiple comparisons (n = 197). The horizontal dashed line indicates p_adjust = 0.05, while vertical dashed lines correspond to a 1.25-fold change. f Two-sided GSEA using C1 positional gene sets on chromosomal bands enriched in KRAS mutated and amplified CRCs vs. CRC with KRAS mutation (n = 197); circle size indicates NES and color the FDR-adjusted p-value. An ideogram of chromosome 12 is provided for orientation. g Enrichr analysis of upregulated genes in CRCs with KRAS mutation and amplification vs. CRC with KRAS mutation (n = 197), using enriched genes on chromosome 12 and other chromosomes. Circle size reflects the Enrichr score and color indicates the Benjamini-Hochberg-adjusted p-value. The Gene Ontology and MSigDB Hallmark genes sets were tested. h. Winter hypoxia score comparison using two-sided Wilcoxon rank-sum test (n = 45 and n = 106 for CRC with KRAS mutation and amplification and CRC with KRAS mutation, respectively; (1) p_adjust = 8.75 × 10⁻⁴). i Heatmap of median ssGSEA scores using KEGG_MEDICUS gene sets, with Benjamini-Hochberg adjusted p-values from Kruskal-Wallis tests encoded in grayscale (n = 412). m&a, KRAS mutated and amplified CRCs; mut, KRAS mutated CRCs; wt, KRAS wildtype CRCs; del, KRAS deletion CRCs (without KRAS mutation). j Box plots of ssGSEA scores for selected gene sets, analyzed by Kruskal-Wallis tests with Dunn posthoc test (two-sided) and Benjamini-Hochberg adjustment (n = 64 for KRAS mutated and amplified CRCs; n = 133 for KRAS mutated CRCs; n = 181 for KRAS wildtype CRCs; n = 34 for KRAS deletion CRCs (without KRAS mutation); (1) p_adjust = 6.78 × 10⁻³; (2) p_adjust = 1.69 × 10⁻³; (3) p_adjust = 1.59 × 10⁻⁷; (4) p_adjust = 6.19 × 10⁻⁴; (5) p_adjust = 5.27 × 10⁻⁴; (6) p_adjust = 3.20 × 10⁻⁶; (7) p_adjust = 3.12 × 10⁻²; (8) p_adjust = 7.63 × 10⁻³; (9) p_adjust = 1.17 × 10⁻³). Analyses compare CRCs with KRAS mutation and amplification to those with KRAS mutation alone (b, d–h) as well as to wildtype CRC and CRC with KRAS deletions without KRAS mutation (i, j). (b, h, j) Data are presented with individual data points and a box-and-whisker-plot (thick line, median; box edges, Q1/Q3; whiskers, 1.5 ×  IQR). *, p_adjust < 0.05; **, p_adjust < 0.01, and ***, p_adjust < 0.001. Source data are provided as a Source Data file.

Fig. 5. Transcriptional profile of CRCs with deletions of MLH1 and BRCA1 in the TCGA cohort.

a Volcano plot showing differential gene expression between MLH1 deletion CRC and MLH1 control CRCs, with genes colored by chromosomal arm (inset highlights genes on chromosome 3; horizontal dashed line, p_adjust = 0.05; vertical dashed lines, 1.25-fold change). Significance was determined by two-sided Wilcoxon rank-sum tests with Benjamini-Hochberg adjustment (n = 276). b Two-sided GSEA using C1 positional gene sets indicating that genes encoded on chromosomal arm 3p (ideogram on the left) are significantly downregulated in CRC with MLH1 deletion compared to MLH1 control CRCs (n = 276). Circle size represents NES, color indicates FDR-adjusted p-value. c G:profiler analysis of genes upregulated in MLH1 deletion CRCs, categorized by Gene Ontology (GO), including molecular function (MF), biological process (BP), and cellular component (CC) (n = 276). P-values are adjusted using the g:SCS method. d Differential gene expression of the HOXA and HOXD loci and flanking genes, with circle size based on log₂(fold change) and color reflecting Benjamini-Hochberg-adjusted p-values from two-sided Wilcoxon rank-sum tests (open circles indicate p_adjust > 0.05; n = 276). e CNV analysis of the HOXA and HOXD loci (two-sided Fisher exact test with Benjamini-Hochberg correction; n.s. not significant, (1) p_adjust = 0.757; (2) p_adjust = 0.941; n = 276)). f Volcano plot visualizing differential gene expression for BRCA1 deletion CRC vs. BRCA1 control CRC (n = 213). Inset shows genes on chromosome 17. Significance determined by two-sided Wilcoxon rank-sum tests with Benjamini-Hochberg adjustment. g Two-sided GSEA using C1 positional gene sets for chromosome 17 (ideogram) reveals that genes encoded on 17q are significantly downregulated in CRCs with BRCA1 deletion compared to BRCA1 control CRC (n = 213). Circle size represents NES, color indicates FDR-adjusted p-value. h G:profiler analysis of genes upregulated in BRCA1 deletion CRCs compared to BRCA1 control CRC (n = 213) with GO, KEGG (Kyoto Encyclopedia of Genes and Genomes), and Reactome (REAC) pathways (p-values adjusted by g:SCS method). i Differential gene expression of the HOXA locus and surrounding genes, with the top showing upregulated genes in BRCA1 deletion CRCs and the bottom showing downregulated ones; circle size reflects log₂(fold change), while color indicates Benjamini-Hochberg-adjusted p-values from two-sided Wilcoxon rank-sum test (open circles denote p_adjust > 0.05; n = 213). Panels a–e compare CRCs with MLH1 deletion to control CRCs, while (f–i) compare CRCs with BRCA1 deletion to control CRCs. Source data are provided as a Source Data file.

Fig. 6. Oncogenetic tree models of metastatic CRC based on maximum likelihood estimation.

Analyses were done using the CRCTropism cohort for the following metastasis sites: a. liver (n = 98), b. lung (n = 69), and c. brain (n = 33). CIs were classified into nine categories from early to late according to proximity to the root (gray boxes). Organotropic events are marked by asterisks. Source data are provided as a Source Data file.

Fig. 7. Reconstruction of phylogenetic trees using multi-sample cytogenetic data from individual CRC patients of the CRCTropism cohort.

a–h The clinical course (top) and the reconstructed phylogenetic tree (bottom) are shown, with the tree scaled according to the number of CIs. Maximum parsimony was used for tree inference. CI events are colored based on the classification introduced in Fig. 6. If a particular event could not be uniquely linked (indicated by brackets), it is displayed at the earliest suggested occurrence. High-level amplifications are indicated in bold. The site of the primary tumor (p) and metastases (m, metachronous metastasis; s, synchronous metastasis) is also listed. Dx, diagnosis of primary cancer; Met, metastasis. Source data are provided as a Source Data file.