Introduction

Colorectal cancer (CRC) remains one of the major contributors to cancer incidence and mortality in Western countries. Recent global estimates indicate nearly 2 million new CRC diagnoses annually, with close to 1 million deaths worldwide. Incidence rates in high-income regions are among the highest globally, driven by lifestyle factors (e.g., diet, physical inactivity) and modulated by the introduction and expansion of screening programs [1]. Nevertheless, despite preventive strategies and improvements in early detection, CRC-related mortality remains substantial [2].

Metastatic disease accounts for a large proportion of this burden. Approximately 30% of patients present with synchronous metastases, and many others will develop metachronous dissemination during the course of their illness [1, 2]. The liver is the predominant metastatic site, affected in roughly half of all metastatic CRC cases, followed by the lungs, lymph nodes, and peritoneum. The prognosis for stage IV CRC remains poor: although advances in combination chemotherapy, targeted biologics [including anti-vascular endothelial growth factor (VEGF) and anti-epidermal growth factor receptor (EGFR) agents], and metastasis-directed surgical approaches have extended survival for selected patients, a durable cure is uncommon [3].

Historically, metastatic CRC was regarded as a single clinical entity, uniformly characterized by rapid systemic dissemination and poor outcomes. However, accumulating clinical evidence has challenged this uniform view. A subset of patients with limited metastatic burden experience unexpectedly favorable outcomes, including long-term disease control and, in some cases, prolonged survival after local treatment of metastases [4]. In contrast, others exhibit rapid, widespread progression despite intensive systemic therapy. This clinical divergence underpins the concept of oligometastatic (OM) disease, a condition in which the tumor has seeded only a limited number of metastases, often accompanied by more indolent biological behavior that may render curative-intent local interventions feasible [5]. The distinction between OM and polymetastatic (PM) disease has practical and biological implications. Clinically, identifying OM disease can inform decisions regarding surgical resection, stereotactic body radiotherapy (SBRT), or other metastasis-directed strategies. Biologically, the existence of a true OM state points to underlying molecular and microenvironmental factors that constrain metastatic dissemination, as well as to evolutionary and genetic trajectories that favor metastatic containment over widespread progression.

Cancer can be conceptualised as a dynamic evolutionary process, shaped by the interplay between genomic instability, microenvironmental constraints, and therapeutic selective pressures, which collectively drive tumor progression and clinical heterogeneity [6]. Recent advances in genomics and immunology have started to elucidate the molecular architecture of OM CRC. Prior studies suggest that both tumor-intrinsic genetic features and host immunologic characteristics may contribute to this phenotype [7, 8]. In particular, variations in the mutational landscape of primary and metastatic lesions, including regressive trajectories of key driver genes, have been implicated in metastatic containment. Additionally, emerging evidence indicates that the tumor microenvironment and host immune competence may exert selective pressures capable of eliminating or restraining aggressive subclones [8]. Therapy may also modulate the evolutionary dynamics of metastatic disease, with a complex interplay between genomic evolution, treatment-induced selective pressures, and the immune contexture collectively shaping metastatic behavior and, ultimately, the clinical phenotype [9, 10]. However, replication stress—arising from persistent proliferative signaling and defective DNA damage repair—is increasingly recognized as a fundamental driver of tumor evolution and clonal diversification [11]. Excessive stress may promote genomic plasticity and aggressive dissemination, whereas insufficient stress adaptation may activate checkpoint responses and senescence-like programs that constrain expansion. Cellular senescence, together with immune surveillance and microenvironmental pressures, may therefore contribute to maintaining a biologically contained metastatic configuration. Within this framework, OM and PM disease can be interpreted as divergent evolutionary outcomes along axes of genomic stress handling and adaptive resilience, rather than merely quantitative differences in tumor burden.

We explored genotype-phenotype relationships in a well-defined and homogeneous cohort of metastatic colon cancer (mCC) patients. To minimize technical variability, all cases were profiled using the same next-generation sequencing platform. By comparing patients with OM versus PM presentations, our objective was to identify genomic features that may serve as biomarkers of OM behavior and to provide insights with potential relevance for personalized therapeutic approaches.

Materials and methods

Results

Clinical outcome validation and molecular stratification of OM versus PM disease

Consistent with the expected clinical behavior of OM, this cohort exhibited a markedly superior OS compared with the PM cohort, thereby supporting the appropriateness and robustness of the patients’ selection (Figure 1).

Display full size

Figure 1. Kaplan-Meier survival curves for oligometastatic (OM) and polymetastatic (PM) colorectal cancer patients. Survival probability (Y-axis) is plotted over time (X-axis) with 95% confidence intervals (CIs). OM patients are shown in dark blue, and PM patients in green. Hazard ratio (HR) and 95% CI are reported in the results section. The survival curves are significantly different, with a log-rank P < 0.0001.

After a median follow-up of 52 months, a total of 55 death events were recorded (OM versus PM: 16/50 versus 39/77). Median OS was not reached in the OM group, whereas it was 29 months in PM patients. Survival analysis confirmed a highly significant difference between the two cohorts, with OM patients experiencing a substantially reduced risk of death (HR = 0.1974; 95% CI: 0.1095–0.3559; P < 0.0001 by log-rank test). This pronounced survival advantage reinforces the notion that long-term survival represents a defining clinical hallmark of the OM phenotype and confirms that the two populations are well differentiated from a prognostic standpoint. A multivariate analysis confirmed that the prognostic impact of the OM phenotype is strong and remains independent of other clinical and/or molecular variables (Table S2).

Building upon this validated clinical stratification, we next investigated whether such phenotypic divergence was mirrored by differences in the distribution of key molecular biomarkers implicated as oncogenic drivers in this disease (Table 2).

Table 2. Distribution of key molecular alterations stratified by metastatic status.

Molecular variable	No. (%)	OM	PM	P
APC status
Wt	40 (31.5)	10	30	0.0250
Mut	87 (68.5)	40	47
BRAF p.V600E status
Absent	118 (92.9)	49	69	0.0730
Present	9 (7.1)	1	8
BRAF all mutations
Wt	115 (90.6)	47	68	0.2862
Mut	12 (9.4)	3	9
ERBB2 amplification
Yes	10 (7.9)	2	8	0.1932
No	117 (92.1)	48	69
PIK3CA status
Wt	92 (72.4)	35	57	0.6212
Mut	35 (27.6)	15	20
RAS status
Wt	69 (54.3)	33	36	0.0341
Mut	58 (45.7)	17	41
TP53 status
Wt	43 (33.9)	14	29	0.2628
Mut	84 (66.1)	36	48
MSS
Stable	117 (92.1)	42	75	0.0064
Unstable	10 (7.9)	8	2
TMB (mut/Mb)
≥ 10	35 (27.6)	19	16	0.0345
< 10	92 (72.4)	31	61

Display full size

APC: adenomatous polyposis coli; BRAF: v-Raf murine sarcoma viral oncogene homolog B1; MSS: microsatellite stability; PIK3CA: phosphatidylinositol-4,5-bisphosphate 3-kinase catalytic subunit alpha; RAS: Rat sarcoma viral oncogene family (including KRAS and NRAS); TMB: tumor mutational burden; TP53: tumor protein p53; OM: oligometastatic; PM: polymetastatic; ERBB2: v-Erb-B2 Avian Erythroblastic Leukemia Viral Oncogene Homolog 2; Wt: wild-type.

While most molecular alterations were similarly distributed between OM and PM patients, several noteworthy patterns emerged. v-Raf murine sarcoma viral oncogene homolog B1 (BRAF) p.V600E mutations were numerically more frequent in the PM cohort (8 versus 1 cases), although this difference did not reach statistical significance (P = 0.0730). Likewise, v-Erb-B2 Avian Erythroblastic Leukemia Viral Oncogene Homolog 2 (ERBB2) amplification, phosphatidylinositol-4,5-bisphosphate 3-kinase catalytic subunit alpha (PIK3CA) mutations, and TP53 alterations showed no significant differences between groups. In contrast, APC mutations were significantly more prevalent in OM disease, whereas APC wild-type tumors were relatively enriched in the PM group (P = 0.0250). RAS mutations, conversely, were significantly more frequent in PM patients compared with OM patients (41 versus 17 cases; P = 0.0341). Microsatellite stability (MSS) status also differed markedly between cohorts, with MSI-high tumors being significantly enriched in OM patients (P = 0.0064). Consistently, a high TMB (TMB ≥ 10 mut/Mb) was more commonly observed in the OM group, whereas low TMB (< 10 mut/Mb) predominated in PM patients (P = 0.0345).

Genotype-phenotype correlation with Phenolyzer and GO analyses

To further characterize the mutational landscape associated with metastatic burden, we analyzed the 20 genes harboring the most frequent Tiers 1–3 mutations in the two study groups. Within this framework, we specifically focused on identifying group-restricted alterations, herein referred to as private mutational events, defined as mutations exclusively detected in one group but absent in the other. Using this approach, exclusive mutations in ERCC1, MST1, myeloblastosis oncogene (MYB), and nuclear receptor corepressor 1 (NCOR1) were identified in the PM cohort, whereas FAT1, isocitrate dehydrogenase 2 (IDH2), low-density lipoprotein receptor-related protein 1B (LRP1B), and paired box 8 (PAX8) emerged as exclusively mutated genes in the OM group. The gene-specific mutational burden for these private events across the analyzed cohort is reported in Table S3 and Table S4.

To explore potential functional or biological relationships among the genes most frequently altered in the two groups, a Phenolyzer-based network analysis was performed (Figure 2).

Display full size

Figure 2. Phenolyzer-based molecular characterization of polymetastatic (PM) and oligometastatic (OM) colorectal cancer. Panels A and C show the interaction networks generated by Phenolyzer for PM (A) and OM (C) patients, highlighting the functional interconnections among prioritized genes in each metastatic phenotype. Nodes represent genes, and edges indicate predicted functional or regulatory interactions. Panels B and D display the corresponding gene prioritization as bar plots for PM (B) and OM (D) patients. In both panels, the Y-axis lists the genes, and the X-axis shows the corresponding Phenolyzer scores. The bars represent the Phenolyzer score for each gene, reflecting both its predicted relevance to metastatic colorectal cancer and its connectivity within the interaction network (see methods).

Despite the qualitative differences in the genes identified as group-specific, this analysis did not reveal distinct relational or network patterns between the two metastatic settings. In both cohorts, APC, TP53, and EGFR consistently emerged as the most highly interconnected nodes within the gene interaction network, indicating a shared core of molecular interdependencies across different metastatic phenotypes.

GO analysis revealed a selective enrichment of molecular characteristics related to cellular energy handling and signal transduction in the PM versus OM patients (Figure 3). Specifically, the most significantly overrepresented GO terms belonged to the MF category and included ATP binding, adenyl ribonucleotide binding, adenyl nucleotide binding, and protein kinase activity.

Display full size

Figure 3. Gene Ontology enrichment analysis comparing poly-metastatic and oligometastatic colorectal cancer. The figure shows statistically significant enriched pathways (Y-axis) plotted against the enrichment score (X-axis), which quantifies the degree of pathway overrepresentation in poly-metastatic compared with oligo-metastatic disease.

Discussion

The concept of OM disease has progressively evolved from a purely descriptive clinical entity into a biologically grounded state, characterized by limited metastatic dissemination and a more indolent natural history compared with widespread PM disease [10]. In this study, we sought to interrogate the clinical, molecular, and functional correlates of OM versus PM mCC in order to better understand whether these two phenotypes represent quantitatively different stages along a metastatic continuum or qualitatively distinct biological conditions. Our results provide several important suggestions that collectively support the latter hypothesis, while simultaneously highlighting the complexity underlying metastatic behavior. However, while our data are consistent with a biologically distinct OM configuration, they should not be interpreted as demonstrating biological inevitability; rather, they suggest probabilistic tendencies within a dynamic evolutionary landscape.

From a clinical standpoint, the marked survival advantage observed in OM patients represents a critical internal validation of cohort selection. The median OS was not reached in the OM group, in contrast to a median survival of 29 months in PM patients, with a nearly 80% relative reduction in the risk of death. This pronounced difference is consistent with the growing body of literature indicating that long-term survival is a defining hallmark of OM disease and reinforces the notion that OM is not merely an early temporal snapshot of PM progression [5]. Rather, the observed outcomes strongly suggest that OM patients harbor intrinsic disease features that limit metastatic competence and/or enable prolonged disease control through local and systemic interventions. At the same time, it remains plausible that OM disease represents a metastable equilibrium—an intermediate evolutionary state maintained by a balance between tumor expansion and host-imposed constraints—rather than a permanently fixed biological endpoint.

Clinico-pathological characteristics further support this interpretation. OM tumors were significantly enriched for low-to-intermediate histological grade, whereas high-grade tumors were overrepresented in the PM cohort, consistent with a more aggressive biological phenotype. Importantly, no significant differences were observed in primary tumor T stage or nodal involvement, underscoring that the divergent metastatic behavior cannot be simply attributed to differences in local tumor extent at diagnosis. Instead, metastatic burden appears to reflect distinct systemic disease properties that manifest downstream of primary tumor growth. The higher utilization of locoregional approaches, such as metastasectomy, stereotactic radiotherapy, and ablative techniques, in OM patients is likely a consequence, rather than a cause, of this more favorable biology, enabling durable disease control when metastatic dissemination remains limited.

Against this clinically validated backdrop, we examined whether molecular differences might underlie the divergent metastatic phenotypes. At the level of canonical driver alterations, several noteworthy patterns emerged. APC mutations were significantly enriched in OM disease, whereas RAS mutations were significantly more frequent in PM disease; BRAF p.V600E mutations also showed a numerical enrichment in PM patients, in line with their well-established association with aggressive tumor biology and poor prognosis [22]. However, aside from these alterations, most established driver events—including ERBB2 amplification, PIK3CA mutations, and TP53 alterations—were comparably distributed between OM and PM cohorts. These findings suggest that while certain oncogenic drivers may contribute to metastatic aggressiveness, driver dominance alone is insufficient to fully explain the stark clinical divergence between the two groups.

More revealing was the differential distribution of MSS and TMB. MSI-high status and elevated TMB were significantly enriched in the OM cohort, whereas PM disease was predominantly characterized by MSS and low mutational burden. This observation is particularly intriguing, as it shifts the interpretative focus away from individual oncogenic lesions toward broader genomic and immunological landscapes. MSI-high tumors are known to generate a higher neoantigen load, leading to enhanced immune recognition and, in many cases, more effective immune-mediated tumor control [23, 24]. In this context, the association between OM disease, MSI, and high TMB raises the possibility that immune surveillance may play a central role in constraining metastatic dissemination, thereby maintaining an OM state. However, this immune-surveillance hypothesis has not been directly tested in tumor tissues from this cohort and should therefore be considered inferential. Dedicated studies and prospective validation are required to substantiate this interpretation. Nevertheless, immune editing is inherently dynamic, and immune-mediated containment may fluctuate over time, potentially allowing eventual escape and transition toward more aggressive phenotypes. Within this evolving framework, additional components of the tumor microenvironment may further modulate the balance between containment and dissemination. Tumor-associated macrophages, for instance, represent key regulators of immune contexture, and their polarization toward an M2-like phenotype has been implicated in immune suppression and tumor progression, including in CRC [25]. In parallel, metastatic dissemination is increasingly recognized as a process sustained by complex intercellular communication networks, including extracellular vesicle-mediated signaling between tumor cells and the surrounding microenvironment. These vesicles can influence pre-metastatic niche formation, immune modulation, and organotropism, thereby contributing to metastatic competence [26]. Integrating these perspectives may help refine the biological interpretation of OM versus PM states, extending beyond intrinsic tumor features to encompass dynamic microenvironmental and systemic interactions.

Despite identifying group-specific “private” mutational events, such as ERCC1, MST1, MYB, and NCOR1 in PM disease and FAT1, IDH2, LRP1B, and PAX8 in OM disease, network-based analysis using Phenolyzer revealed no meaningful differences in the relational architecture of altered genes between the two cohorts. However, the Phenolyzer-based network analysis should be interpreted as exploratory and hypothesis-generating rather than mechanistically definitive. While it enables a systems-level prioritization of genes and functional interconnections associated with distinct metastatic phenotypes, it does not establish causal biological relationships. Dedicated in vitro and in vivo investigations will be required to validate the functional relevance of the identified candidate genes and enriched pathways, and to clarify their potential contribution to metastatic competence and disease trajectory. Among the private genes identified, ERCC1 emerges as the only alteration directly involved in the maintenance of genomic integrity, given its central role in nucleotide excision repair and interstrand crosslink repair. Dysfunction of ERCC1 is mechanistically linked to defective DNA damage resolution and the accumulation of genomic lesions, thereby promoting genomic instability. However, it is important to note that ERCC1 loss may be at least partially functionally compensated by alternative DNA repair pathways or redundant repair complexes, potentially attenuating its phenotypic impact on tumor behavior in certain molecular contexts [24, 27]. In contrast, LRP1B does not act as a primary driver of genomic instability but rather represents a marker of an intrinsically unstable genomic background. LRP1B-mutant tumors are noteworthy for exhibiting a higher TMB and neoantigen load, along with upregulation of antigen-presentation and inflammatory gene signatures and increased CD4⁺/CD8⁺ T-cell infiltration—features that plausibly enhance tumor immunogenicity and have been associated with improved clinical courses [28–30]. Despite these differences, in both OM and PM disease, APC, TP53, and EGFR emerged as the most highly interconnected nodes, reflecting a shared core of molecular dependencies typical of CRC biology. This convergence suggests that, at a systems level, both phenotypes rely on similar oncogenic frameworks, despite diverging clinically. Importantly, in the absence of mechanistic validation, the identification of these biological drivers should be interpreted as hypothesis-generating rather than definitive.

Furthermore, GO enrichment analysis revealed a selective functional skewing in PM disease toward molecular programs related to energy utilization and signal transduction. The overrepresentation of GO terms such as ATP binding, adenyl ribonucleotide/nucleotide binding, and protein kinase activity points to an increased reliance on energy-dependent processes and kinase-driven signaling cascades in PM tumors. Collectively, these features are hallmarks of heightened proliferative capacity, metabolic adaptability, and dynamic signal integration, all of which are critical for sustaining widespread metastatic growth across multiple anatomical sites. Interestingly, PM tumors display a kinase-enriched, RAS-skewed, and bioenergetically demanding molecular profile that is permissive for sustained Nuclear Factor kappa-light-chain-enhancer of activated B cells (NF-κB) activation, potentially reinforcing survival signaling, stress tolerance, metastatic plasticity, and stemness [31]. However, it should be emphasized that this inference is based on functional enrichment rather than direct measurement of pathway activation, and alternative explanations—such as microenvironmental cytokine gradients, stromal interactions, or hypoxia-induced signaling—may also contribute to this phenotype. In contrast, the absence of comparable GO enrichment signals in OM disease aligns with a less metabolically and signaling-intensive phenotype. This relative functional restraint is coherent with the clinical and genomic features observed in OM patients, including higher MSI/TMB and prolonged survival, and may reflect a tumor state that is more effectively constrained by immune surveillance and less dependent on aggressive bioenergetic and kinase-driven programs. Yet, this restrained profile may not represent a stable biological endpoint; rather, it could correspond to a transient balance between proliferative drive and external constraint, susceptible to disruption by additional genomic events or environmental pressures. Importantly, these GO-level differences emerge despite a shared core oncogenic network, suggesting that metastatic behavior may be modulated not by distinct pathway architectures, but by differential functional weighting of common molecular components. Such modulation likely reflects an interplay between intrinsic tumor cell programs and extrinsic forces, including immune editing, stromal conditioning, vascular dynamics, and therapy-induced stress. Overall, the GO analysis provides a functional bridge between genomic observations and clinical behavior, indicating that PM disease could be characterized by a molecular bias toward energy-intensive and kinase-mediated processes that plausibly support extensive metastatic competence. Conversely, OM disease appears to lack this functional enrichment, reinforcing the notion that OM biology could represent a qualitatively distinct, system-level state rather than a mere quantitative reduction of metastatic spread.

Collectively, these observations may argue against an overly simplistic model in which OM disease is driven by a unique set of dominant oncogenic pathways. Rather, they could support a more nuanced interpretation in which OM and PM disease diverge through higher-order biological mechanisms integrating genomic complexity, tumor-host interactions, and immune dynamics. Within this framework, driver mutations may contribute to establishing the oncogenic potential of the tumor, while the ultimate metastatic trajectory might be influenced by additional layers of regulation that extend beyond the cancer cell itself.

Several limitations of this study should be acknowledged. The retrospective nature of the analysis and the sample size may limit the statistical power to detect subtle molecular differences, particularly for low-frequency events. Moreover, our analyses were based primarily on genomic data, without direct assessment of immune infiltration, spatial heterogeneity, or functional immune activity. An additional limitation of the genotype-phenotype analysis is the absence of validation in an independent or publicly available dataset; therefore, external validation in future studies is warranted. Another point deserving consideration concerns the interpretation of bioinformatic enrichment analyses derived from large transcriptomic repositories and curated genomic datasets. Differential expression and pathway enrichment results may be influenced by technical and compositional biases inherent to bulk data, including sample heterogeneity and variability in preprocessing pipelines. These factors can introduce distortions that affect biological inference and limit the robustness of conclusions. Accordingly, our findings should be interpreted with caution, acknowledging the intrinsic limitations of bulk transcriptomic analyses as highlighted in prior methodological studies [32, 33]. Future prospective studies integrating multi-omics approaches—including transcriptomic, proteomic, immunophenotypic, and spatial profiling analyses—will be essential to more comprehensively define the biological basis of the OM phenotype. Applied to well-annotated and adequately powered cohorts, such integrative strategies would enable a more in-depth functional characterization of the candidate genes and molecular networks emerging from our findings, thereby clarifying their biological significance and potential clinical implications. Importantly, although our analysis captures baseline genomic features of primary tumors, metastatic behavior likely emerges from a dynamic interplay between intrinsic genomic architecture and subsequent stress-adaptive responses induced by treatment and other external exposures. A more comprehensive characterization of these external interferences, including systemic therapies, locoregional interventions, inflammatory cues, metabolic stress, and even non-oncologic perioperative factors, will, therefore, be necessary. Such influences may reshape tumor cell states, modulate stem-like compartments, and alter stress tolerance and survival programs, ultimately impacting metastatic trajectories [34]. Disentangling baseline genomic determinants from environmentally induced adaptive reprogramming will be critical to achieve a biologically integrated understanding of OM versus PM evolution.

In conclusion, our findings reinforce the concept that OM and PM CRC represent clinically and biologically distinct entities. While differences in specific driver alterations are detectable, they do not fully account for the profound divergence in clinical outcomes. Instead, broader biological and immunological features, closely linked to MSI status and TMB, appear to play a pivotal role in shaping metastatic behavior. Advancing our understanding of these mechanisms will be critical not only for refining prognostic stratification but also for identifying patients in whom limited metastatic disease truly reflects an indolent, controllable state rather than an early manifestation of aggressive systemic spread.