Clinical characteristics of the patients

The summary of the main characteristics of all patients is in Supplementary Table S1. The median age at the time of mCLM diagnosis was 65 years (range 37–78) and the group contained more men (61%) than women (39%). This unequal sex distribution is in agreement with the reported higher CRC incidence in men compared to women [1]. After the curative mCLM surgery, the median of RFS was 11 (range 1–103) months. Seven patients experienced no progression after mCLM treatment and except one had long OS > 3 years. The reported median OS of patients with CRC liver metastases after curative surgery is 24 to 37 months in the last ten years [25]. In agreement with previous studies, the median OS of our dataset was 40 (3–103) months. After mCLM resection, patients were treated predominantly with regimens FOLFOX, CAPOX, FOLFIRI, OR CAPIRI alone or with targeted therapy (cetuximab, panitumumab, or bevacizumab based on RAS/BRAF mutation screening results). Patients survival between primary tumor and mCLM resection was not affected by pTNM, stage, or grade of primary tumor or by administration of adjuvant treatment (none vs. administered). Survival between mCLM resection and second relapse or death was unaffected by the above-mentioned factors, the mCLM resection radicality, or administered chemotherapy after mCLM resection (p > 0.05). A trend towards a better outcome of the patients with single compared to those with two or more mCLM loci was observed (p = 0.08).

General description of the whole exome profiling

The average coverage was 238x for mCLM and 68x for non-malignant liver. On average, 83% of bases in metastases (41% for non-malignant tissue) were covered at least 30x, and 94% (88%, respectively) at least 10x. The duplicate rate was 53% for metastases and 37% for non-malignant liver tissue.

Somatic profile of mCLM

The total number of detected variants per mCLM sample was 1 528 ± 798 on average (ranging from 966 to 5 292, median 1 388). The amount of somatic variants fulfilling the filtering criteria (see Materials and Methods) per sample was 370 ± 464 (ranging from 82 to 2 932, median 281). From these, 7 952 silent variants were then excluded from downstream analyses. The distribution of pathogenic somatic variants in gene coding regions among 41 patients in our cohort is shown in Fig. 1. The most common class of somatic variants was the missense mutation (Fig. 1A), and the most common type was a single nucleotide variant (SNV, Fig. 1B). The most common nucleotide substitution was the C > T transition (Fig. 1C). Two patients had TMB-high status by definition and notably differed by the mutational load (total numbers of non-silent variants were 1 508 and 616, respectively) from the rest (Fig. 1D). The overall mutational summary for all samples is in Supplementary Table S2. From the list of 20 FLAGS (FrequentLy mutAted GeneS) genes that are known to be frequently mutated in cancer but are unlikely to be pathogenic [26], TTN, AHNAK2, SYNE1, MUC16, and OBSCN were found among genes altered at ≥ 20% in mCLM. Due to their FLAGS status, these genes will not be discussed further.

Fig. 1

The summary of the distribution of the overall variants in mCLM. Only protein-changing variants were considered, 7 952 silent variants were excluded from the analysis. (A) The classification of variants according to their functional effect (missense mutation, frameshift deletion/insertion, nonsense mutation, splice site, inframe deletion/insertion, translation start site, or nonstop mutation). On the x-axis, the counts are in the log scale. The most prevalent variants were missense. (B) The types of variants (TNP stands for trinucleotide variant; SNP, single nucleotide polymorphism; INS, insertion; DNP, dinucleotide variant; DEL, deletion). On the x-axis, counts are in the log scale. The most common type of variant was the SNP. (C) The type of nucleotide substitution. The most frequent substitution was the C > T transition. (D) The counts and distribution of the variants for the indicated samples; the dashed line represents a median (124 variants per sample excluding silent variants). Two TMB-H patients are separated to deflate the y-axis

From the rest of the genes, the most frequently mutated gene was TP53, which was detected in 76% of patients (31/41; 34 variants in total) (Fig. 2). Missense was the most common variant type in TP53, followed by nonsense, frameshift indels, and splice site variants (lollipop in Supplementary Fig. S1B). The second most frequently mutated gene was APC, with variants detected in 66% of patients (27/41; 42 variants in total) in our cohort (Fig. 2). The most frequent APC alterations were nonsense variants followed by frameshift deletions, missense variants, and frameshift insertions (lollipop in Supplementary Fig. S1A). In other genes, e.g., KRAS the missense variants prevailed (lollipop in Supplementary Fig. S1C).

Fig. 2

Oncoplot of most mutated 20 genes (Top 20)

Analysis of mutation co-occurrence (Supplementary Fig. S2) revealed that ABCA13 were altered together with MYT1L (p < 0.01), FREM2, UNC80, and PIK3CA (all p < 0.05). Variants in UNC80 also co-occurred with alterations in FBXW7 (p < 0.01) and PCLO (p < 0.05). Furthermore, the alterations in FREM2 co-occurred with variants in PIK3CA, COL6A5, and HRNR (all p < 0.05) and those in HRNR with alterations in EYS and FAT4 (both p < 0.05). Moreover, alterations in TP53 and FREM2 were mutually exclusive (p < 0.01). Interestingly, alterations in TP53 and PIK3CA and in KRAS and RYR2 were mutually exclusive too (p < 0.05).

We then performed an analysis of CNV in mCLM samples. On average, the tumors bore 55.9 ± 30.4 CNVs (ranging from 16 to 151, median 53). The average size of the deletions/amplifications was 11.7 ± 5.0 Mb. Most common CNVs were single-copy amplifications (23.3% of all CNVs), followed by single-copy deletions (21.8%). Less common were CNVs with > 3 copies (9.5%) and homozygous losses (1.2%) (Supplementary Table S4).

Next, relative contributions of reference SBS mutational signatures (COSMIC database v3.3) were determined for each mCLM sample. All 79 available reference SBSs were assessed in the first round (Supplementary Table S5). The top 10 signatures by overall contribution (in the whole cohort) were refitted to get their final relative contributions, except for SBS54, which was excluded due to being a suspected sequencing artefact by COSMIC. Reducing the mutational catalog to 9 signatures (Fig. 3, Supplementary Table S6) led to a negligible decrease in fitting accuracy in the vast majority of samples, as demonstrated by cosine similarity of the original mutational catalogs and the reconstructed catalogs from the top 9 signatures and from all 79 reference signatures (Fig. 3). A high correlation (r2 > 0.5, p < 0.001) between signature pairs SBS1-SBS15 (negative correlation) and SBS24-SBS31 (positive) in mCLM was observed.

Fig. 3

Relative contribution of top mutational signatures in mCLM. Top – relative exposures of top 9 SBS signatures in each sample. Bottom: top_SBSs – cosine similarity of the mutational catalog reconstructed from 9 most dominant signatures (above) to the original data. full_COSMIC_v3.3 – cosine similarity when fitting the full set of 79 COSMIC reference signatures (version 3.3). Few samples showed a substantial decrease in fitting accuracy by reduction of the number of signatures from 79 to 9

Germline profile of non-malignant liver samples

We performed the analysis of rare (allele frequency in gnomAD < 0.05) and deleterious (stop-gain, stop-loss, frameshift insertion or deletion, and changing the splice site or transcription start site predicted by VEP or listed as pathogenic in ClinVar or CADD score > 25) germline variants similarly to mCLM. In total, we found 5 573 variants. The median count was 184 (140–248) per patient. The median count of non-silent variants per patient was 163 and the most frequently mutated genes were CTU2, GGT3P, and AGAP6 (Fig. 4A, B, Supplementary Table S7). Interestingly, these and some other genes, e.g., CCDC7, PRAMEF10, ZNF101, SPTBN5, or DHRS4L2 had a high rate of variants with predicted HIGH functional effect (lollipops in Supplementary Fig.S3). In the case of CTU2, all samples had the same variant (rs11278302) with a HIGH functional effect (splice donor variant) predicted. Due to its high frequency (95%), it was further not considered clinically relevant.

Fig. 4

Distribution of rare and deleterious germline variants in non-malignant liver of mCLM patients. (A) plot depicting TOP20 most mutated genes fulfilling filtering conditions, (B) load of germline variants (median 163)

We also provide co-occurrence analysis of altered genes, together with somatic mutations. Except for already depicted interactions between somatic variants (Supplementary Fig. S2), germline variants in AGAP6 or ZNF101 co-occurred with somatic mutations in KRAS (p = 0.048 and p = 0.045, respectively) while CATSPER2 was mutually exclusive (p = 0.002). MOL7 was mutually exclusive to APC (p = 0.042). Germline variants in several genes co-occurred mutually too, e.g., CDH23 and ZNF101 (p = 0.007) (Supplementary Fig.S4).

Clinical relevance of germline and somatic profiles of the patients

Although there is no clear consensus about the definition of early recurrence [27], we divided patients by the 6 months RFS cut-off [28] to provide surrogate information about poor response to mCLM therapy. At the time of analysis, for 12 patients the relapse appeared within six months after curative surgery. Patients in this subgroup, had notably more variants in ABCA13 and UNC80, and fewer in APC and RYR1 (Supplementary Fig.S5A) but none of these differences were statistically significant. Moreover, patients with short RFS seemed to have exclusively nonsense variants in the TP53 tetramer domain (n = 2 vs. none in long RFS, Supplementary Fig.S5B), although this could be by chance and should be subject to validation using larger datasets. Stratification of patients according to subsequent systemic therapy after mCLM resection was not possible due to high heterogeneity.

Of the whole sample set, only two patients, 8164T and 7960T, were classified as TMB-high and had considerably high MSI status; however, only 7960T was classified MSI-high based on the strict 20% cut-off. Both had the MMR-D status using their somatic mCLM profile. One patient suffered from early recurrence and died six months after surgery. The other had OS of 82 months without recurrence signs by the censoring time point. Thus, these characteristics do not seem to provide prognostic information in our patient set.

Further, we analyzed prognostic associations of somatic variants in frequently altered genes and enrichment of any of the thirteen oncogenic pathways previously identified in CRC [29]. For this comparison, we used separately patients with variants classified as having a HIGH/MODERATE, or exclusively HIGH predicted functional effect against those without such variants. One patient was excluded from both RFS and OS analyses as lost to follow-up and two other patients were excluded from RFS analyses due to lung metastasis and absence of relapse-free period. Thus, 38 patients entered RFS and 40 OS analyses.

Most importantly, patients with HIGH functional effect only somatic variants in homologous recombination repair (HRR) genes [30] (n = 5; one in ATR, BRCA2, NBN, and two in BRCA1, Supplementary TableS8) had significantly prolonged RFS compared to those without such variants (p = 0.021, Fig. 5A). Including patients with pathogenic germline variants (two in ATM and one in RAD51C) showed the same trend for RFS (n = 8, p = 0.040, Fig. 5B). OS analysis was not significant (p = 0.903 and 0.400, respectively). Concerning MMR genes (MLH1, MLH3, MSH2, MSH3, MSH6, and PMS2 [31]), five patients had either somatic profile/HIGH functional effect variant (n = 3, two with TMB-H status and one with MSH3 variant) or pathogenic germline variants (n = 2, both in MLH1). However, no significant prognostic association was observed.

Fig. 5

Kaplan-Meier plots of patient survival stratified by the carriage of variants in HRR, oncodriver pathways, and individual genes. (A) RFS analysis of somatic variants with HIGH functional effect in the HRR gene panel [30]. (B) RFS analysis of (A) complemented with rare pathogenic HRR variants. RFS analyses of somatic variants with HIGH or MODERATE functional effect in the MYC (C), Notch (D), and Hedgehog (E) pathways. OS analysis of somatic variants with HIGH or MODERATE functional effect in the JAK-STAT pathway (F), VIPR2(G), and MUC16 (H). Red line represents patients carrying the variant, and the blue line those without. HR = hazard risk

Patients with HIGH or MODERATE functional effect somatic variants in the MYC pathway (n = 4, MYC, MGA, MLXIP, and MLX) had significantly worse RFS (p = 0.004, Fig. 5C), and those with somatically altered genes in the Notch (n = 16) or Hedgehog (n = 13) pathways had prolonged RFS (p = 0.008 and 0.012, respectively, Fig. 5D, E). Again, these associations were non-significant in OS analyses (p > 0.05). Importantly, patients with somatic variants in the JAK-STAT pathway (n = 14) had significantly prolonged OS (p = 0.009, Fig. 5F), but no association with RFS was found. The list of followed genes in pathways is in Supplementary TableS9.

Regarding individual genes, patients carrying the somatic KRAS G12D variant (n = 6) had significantly poorer RFS compared to the rest of the patients (p = 0.044, Supplementary Fig.S6A). This association was supported by the OS analysis of an external dataset (MSK panel data, n = 97) resulting in the same trend, i.e. shortened OS for patients with KRAS G12D variant compared to the rest (p = 0.020, Supplementary Fig.S6B). No association was observed in patients stratified to wild type only vs. G12D variant (p > 0.05). APC, TP53, FAT4, PIK3CA, FBXW7, or other frequently mutated genes listed in the Cancer Gene Census (https://cancer.sanger.ac.uk/census) were not prognostic individually. However, patients bearing somatic variants in VIPR2 (n = 6, missense) had significantly poorer OS than patients without such variants (p < 0.001, Fig. 5G). On the other hand, patients with somatic variants in MUC16 (n = 8, one nonsense, and seven missense variants), encoding the CA125 tumor antigen, had significantly prolonged OS (p = 0.038, Fig. 5H). No patient had pathogenic germline or HIGH/MODERATE functional effect somatic variants in oncodrivers BRAF or NRAS that would be informative about the therapy. One patient had a somatic variant in POLE, together with STK11, one in PTEN, one in MUTYH, and two in SMAD4.

We then performed survival analyses of patients stratified using the TMB divided by the median and CNV number, size, and types. None of them was significant.

Next, survival analyses were done with top 9 SBS signatures divided by the median relative contribution but again no significant association was found.

Finally, we analyzed the prognostic significance of rare and deleterious germline variants in genes listed in Fig. 4A. Out of 19 individual genes and 13 oncodriver pathways fulfilling the above conditions (Supplementary TableS9), only two were prognostic. Patients with alterations in CCDC7 (n = 11) had significantly poorer RFS compared to wild-type carriers (p = 0.040, Supplementary Fig.S7A) and similarly, patients with wild-type KMT2E (n = 31) had significantly worse OS than variant carriers (p = 0.017, Supplementary Fig.S7B). Six patients had germline frameshift deletion p.Lys1250Aspfs*8 (rs371466318) in CCDC7, while the rest of the patients with alterations had frameshift insertions (n = 2, p.Leu77Phefs*6, rs146679927), or other deletions (n = 2, p.Asn919Leufs*11, rs202220321 and one novel p.Asn861Ilefs*21). The most frequent alteration did not influence RFS or OS. None of the CCDC7 variants were reported in ClinVar. For KMT2E, all patients had the missense variant pGly999Cys (rs117986340), deleterious according to the SIFT predictor but benign according to ClinVar.