Geographic and age-related variations in mutational processes in colorectal cancer

Study design

981 colorectal cancers (45.7% female) were collected from intermediate-incidence countries with ASRs of 13-20/100,000 people (Iran, Thailand, Colombia, Brazil) and high-incidence countries with ASRs >24 (Argentina, Canada, Russia, Serbia, Czech Republic, Poland, Japan), including the highest ASR of 37 in Japan¹ (Fig. 1a; Supplementary Table 1). Of the 981 cases, 320 were from the proximal colon, 333 from the distal colon, 326 from the rectum, and 2 from unspecified subsites (Fig. 1b). There were 132 early-onset cases, which were 1.88-fold enriched in the distal colon and rectum compared to the proximal colon (p=0.006). All cancers and their matched normal samples underwent whole genome sequencing, achieving a median coverage of 53-fold and 27-fold, respectively.

• Download figure
• Open in new tab

Fig. 1. Geographic, clinical, and molecular characterization of the Mutographs colorectal cancer cohort.

a, Geographic distribution of the 981 patients across four continents and 11 countries, with an indication of the total number of cases as well as the percentage of early-onset (eo) cases below 50 years of age. Countries were colored according to their age-standardized incidence rates (ASR) per 100,000 individuals. The designations employed and the presentation of the material in this publication do not imply the expression of any opinion whatsoever on the part of the authors or their institutions concerning the legal status of any country, territory, city or area or of its authorities, or concerning the delimitation of its frontiers or boundaries. b, Tumor subsite distribution of the cohort across the colorectum, with an indication of the total number of cases and the percentage of early-onset cases. Different subsites were colored according to the percentage of early-onset cases. Additional 2 cases had unspecified subsites. c, Scatter plots indicating the distribution of molecular subgroups across the sequenced tumors according to the total number of single base substitutions (SBS) and small insertions and deletions (indels; ID), as well as the percentage of genome aberrated (PGA). Cases for which tumor purity was insufficient to determine an accurate copy number profile or without large copy number alterations (65/981) were excluded from the SBS - PGA panel. d, Box plots indicating the distribution of SBS and ID across early-onset (under 50 years; purple) and late-onset (50 years or older; green) microsatellite stable (MSS) colorectal tumors. Statistically significant differences were evaluated using multivariable linear regression models adjusted by sex, country, tumor subsite, and tumor purity. The line within the box is plotted at the median, while the upper and lower ends indicate the 25^th and 75^th percentiles. Whiskers show 1.5 × interquartile range, and values outside it are shown as individual data points. e-g, Average mutational profiles of early-onset and late-onset MSS colorectal tumors for SBS (SBS-288 mutational context; e), ID (ID-83 mutational context; f), and copy number alterations (CN-68 mutational context; g).

Mutation burden and molecular classification

The 981 colorectal cancers were divided into known molecular subtypes based on their somatic mutation burdens and profiles. Consistent with prior studies^20,28, two main subtypes were identified: DNA mismatch repair proficient cancers, also known as microsatellite stable (MSS), and DNA mismatch repair deficient cancers, often referred to as tumors showing microsatellite instability (MSI). MSS samples (n=802, 81.8%; Fig. 1c) were characterized by a lower burden of single base substitutions (SBS; median: 12,054) and small insertions and deletions (ID; median: 1,451), and a higher burden of large-scale genomic aberrations (median: 53.5% of genome altered). In contrast, MSI samples (n=153, 15.6%) exhibited higher SBS and ID burdens (median: 95,426 and 125,100, respectively) with limited genomic aberrations (median: 7.0%). As expected, the average mutational profiles of MSS and MSI colorectal tumors were different (Extended Data Fig. 1a-b).

MSI samples were predominantly found in the proximal colon (OR=12.2, p=3.8×10^-27) and were more common in early-onset cases (OR=2.6, p=0.001). Notably, 31/153 MSI cases (20.3%), including 13/28 MSI early-onset cases (46.4%), carried germline pathogenic variants in DNA mismatch repair genes consistent with Lynch syndrome (Supplementary Table 2). After excluding all cases attributed to Lynch syndrome, there was no enrichment of MSI cancers in early-onset cases (p>0.05). Deficiencies of other DNA repair mechanisms were observed in 24/981 cancers (2.4%), including ultra-hypermutated cases with mutations in POLE (n=10, 1.0%) and POLD1 polymerases (n=3, 0.3%), homologous recombination deficient (HRD) cases (n=7, 0.7%), and cases with mutations in the base excision repair genes MUTYH (n=1, 0.1%), NTHL1 (n=2, 0.2%), and OGG1 (n=1, 0.1%) (Methods; Supplementary Table 3-4; Supplementary Fig. 1-3).

The mutational catalogues of DNA repair deficient cancers are dominated by somatic mutations resulting from the failed repair process, rendering it difficult to characterize mutational processes unrelated to this failure²⁹. To enable investigation of the latter, we therefore focused the main analyses on DNA repair proficient colorectal cancers, while reporting DNA repair deficient cases in the Supplementary Note. Two cases treated with chemotherapy for prior cancers were also excluded as their mutation profiles were dominated by the mutational signatures of chemotherapy agents^19,30 (Supplementary Fig. 4). The remaining cohort consisted of 802 treatment-naïve DNA repair proficient colorectal cancers, including 97 early-onset cases.

After adjustment for sex, country, tumor subsite, and tumor purity (Methods), early-onset cancers showed reduced burdens of SBS (fold-change [FC]=0.92, p=0.045) and ID (FC=0.90, p=0.018; Fig. 1d) but not of doublet base substitutions (DBS), copy number alterations (CN), or structural variants (SV) when compared to late-onset cases (p>0.05). Nevertheless, the average mutation spectra of early-onset and late-onset cancers were remarkably similar for all types of somatic mutations (cosine similarity>0.97; Fig. 1e-g; Extended Data Fig. 1c-d). Mutation burden also varied substantially for specific countries when compared to all others, including Canada (lower SBS and ID burdens), Poland (higher SBS and DBS), Japan (lower SBS, ID, DBS), Iran (lower ID), and Brazil (higher ID and CN; Extended Data Fig. 2). However, mutation profiles were generally consistent across all countries (Extended Data Fig. 3).

Repertoire of mutational signatures

A total of 16 SBS, 10 ID, 4 DBS, 6 CN, and 6 SV de novo mutational signatures were extracted from the 802 MSS colorectal cancers and subsequently decomposed into a combination of previously reported reference signatures and potential novel signatures (Supplementary Note; Supplementary Tables 5-15). The 16 de novo SBS signatures encompassed 15 COSMICv3.4 signatures (Extended Data Fig. 4a; Supplementary Table 10), including those previously associated with clock-like mutational processes (SBS1, SBS5)³¹, APOBEC deamination (SBS2, SBS13)³¹, deficient homologous recombination (SBS3)³¹, reactive oxygen species (SBS18)³², exposure to the mutagenic agent colibactin synthesized by Escherichia coli and other microbes carrying a ∼40kb polyketide synthase (pks) pathogenicity island (SBS88)^33,34, and mutational processes of unknown causes (SBS8, SBS17a/b, SBS34, SBS40a, SBS89, SBS93, SBS94)^{5,19,32,34,35}. Three previously described signatures of unknown origin²¹ (SBS_F, SBS_H, SBS_M; Extended Data Fig. 4b) and a novel signature (SBS_O; Extended Data Fig. 4c) were also detected. SBS_O corresponds to a refined version of a previously reported signature of unknown etiology (SBS41; Methods)¹⁹. With respect to ID, DBS, CN, and SV, most de novo extracted mutational signatures were highly similar to, or directly reconstructed by, COSMICv3.4 reference signatures (Extended Data Fig. 4d-f and 5; Supplementary Table 10) with the exception of an ID signature (ID_J), characterized by deletions of isolated Ts and insertions of Ts in long repetitive regions resembling a previously reported signature³⁴ (Extended Data Fig. 4e), and three novel signatures from large mutational events (CN_F, SV_B, SV_D; Extended Data Fig. 5b&d), which were extracted due to the extended contexts used in our signature analysis (Methods).

Geographic variation in mutational signatures

Despite the similar mutation profiles across countries (Extended Data Fig. 3), several signatures exhibited varying prevalence when comparing one country to all others (Fig. 2a; Supplementary Fig. 5; Supplementary Table 16). Notably, SBS89 (OR=28.0, q=0.001), DBS8 (OR=8.9, q=3.2×10^-4), and the novel ID_J (OR=9.6, q=6.2×10^-5) were at elevated frequencies in Argentina when compared to all other countries (Fig. 2b). Signatures SBS89, DBS8, and ID_J also showed a strong tendency to co-occur (p<1.7×10^-11) suggesting they may arise from the same underlying mutational process. In Colombia (Fig. 2c), higher frequencies were observed for SBS94 (OR=19.7, q=3.2×10^-5), the novel SBS_F (OR=10.7, q=2.0×10^-4), and DBS6 (OR=12.5, q=0.028) when compared to all other countries, with evidence of co-occurrence of SBS94 with SBS_F (p=0.017) and DBS6 (p=1.9×10^-4). Enrichments were also found for SBS2 (OR=2.0, q=0.041) and SBS_H (OR=2.3, q=0.001) in Russia and CN_F (OR=3.5, q=3.9×10^-4) in Brazil, whereas depletions were identified for DBS2 in Thailand (OR=0.38, q=0.008) and for DBS4 in Colombia (OR=0.06, q=0.034; Fig. 2a). Overall, the results indicate international differences in the prevalence of certain mutational processes involved in colorectal cancer development.

• Download figure
• Open in new tab

Fig. 2. Geographic variation of mutational signatures in microsatellite stable colorectal cancers.

a, Dot plot indicating the variation of mutational signature prevalence in specific countries compared to all others. Statistically significant enrichments were evaluated using multivariable logistic regression models adjusted by age of diagnosis, sex, tumor subsite, and tumor purity. Firth’s bias-reduced logistic regressions were used for regression presenting complete or quasi-complete separation. Data points were colored according to the odds ratio (OR) of the enrichment, with their size representing statistical significance. P-values were adjusted for multiple comparisons based on the total number of mutational signatures considered per variant type and the total number of countries assessed, and reported as q-values. Q-values<0.05 were considered statistically significant and marked in red. b-c, Geographic distribution of the ID_J (b) and SBS_F (c) mutational signatures. Countries were colored based on the signature prevalence. d, Volcano plots indicating the association of mutational signature activities with the age-standardized incidence rates. Statistically significant associations were evaluated using multivariable linear regression models adjusted by age of diagnosis, sex, tumor subsite, and tumor purity. P-values were adjusted for multiple comparisons based on the total number of mutational signatures considered per variant type and reported as q-values. Horizontal lines marking statistically significant thresholds were included at 0.05 (dashed orange line) and 0.01 q-values (dashed red line). e-f, Scatter plots indicating the association of the mutations attributed to the SBS88 and ID18 mutational signatures with the age-standardized incidence rates (ASR) across countries for colorectal cancer (e), and independently for colon and rectal cancers (f). Data points were colored based on signature prevalence, with their size indicating the total number of cases per country. Statistically significant associations were evaluated using the sample-level multivariable linear regression models used in d (e), and similar multivariable linear regression models adjusted by age of diagnosis, sex, and tumor purity (f).

To explore the broader epidemiological implications of international variation in mutational processes, as previously done for kidney cancer⁵, we evaluated the relationships between ASR and mutational signatures (Fig. 2d; Supplementary Table 17). Independent of covariates, colibactin-induced mutational signatures, SBS88 and ID18, as well as clock-like signature SBS1 and novel signature SBS_H, associated with an increasing rate of ASR for colorectal cancer, whereas novel signature CN_F associated with a reduced ASR rate (q<0.05; Fig. 2d-e; Extended Data Fig. 6a). For SBS88 and ID18, the association was linked with the ASR for rectal cancer (q=0.088 and q=0.008; Fig. 2f; Supplementary Table 18). In contrast, for SBS1, SBS_H, and CN_F the association was particularly strong for the ASR of colon cancer (q=0.009, q=0.015, and q=0.057; Extended Data Fig. 6b). Colibactin-associated signatures were also found elevated in patients from countries with high ASR rates for early-onset colorectal cancer (Extended Data Fig. 6c).

Colibactin induced mutational signatures are enriched in early-onset colorectal cancer

In addition to examining the global distribution of mutational signatures, the substantial number of early-onset colorectal cancer cases enabled evaluating the association between mutational signatures and age at diagnosis. Although the average mutation profiles of early-onset and late-onset colorectal cancer cases were similar (Fig. 1e-g), the prevalence of some mutational signatures was associated with the age of diagnosis, independently of country of origin (Fig. 3a; Supplementary Table 19), genetic ancestry or ethnicity (Supplementary Fig. 6-8). As expected, late-onset cases showed enrichment in signatures known to accumulate linearly with age in normal colorectal crypts³⁶, including SBS1, SBS5, ID1, and ID2 (Fig. 3a-b). Unknown etiology indel signatures ID4, ID9, and ID10 also showed associations with late-onset cases (Fig. 3a-b).

• Download figure
• Open in new tab

Fig. 3. Variation of mutational signatures with age of onset in microsatellite stable colorectal cancers.

a, Volcano plots indicating the enrichment of mutational signature prevalence in early-onset and late-onset cases. Statistically significant enrichments were evaluated using multivariable logistic regression models for age of onset categorized in two subgroups (early-onset, <50 years of age; and late-onset, ≥50) and adjusted by sex, country, tumor subsite, and tumor purity. Firth’s bias-reduced logistic regressions were used for regression presenting complete or quasi-complete separation. P-values were adjusted for multiple comparisons based on the total number of mutational signatures considered per variant type and reported as q-values. Horizontal lines marking statistically significant thresholds were included at 0.05 (dashed orange line) and 0.01 q-values (dashed red line). b, Line plots indicating mutational signature prevalence trend across ages of onset, using five different age groups. Signatures significantly enriched in early-onset or late-onset cases (as shown in a) were colored in purple and green, respective, whereas signatures not varying significantly with age were colored in grey. c, Bar plots indicating mutational signature prevalence across age groups, with indication of the total number of cases where signatures were detected. Statistically significant trends were evaluated using multivariable logistic regression models for age categorized in five subgroups (0-39, 40-49, 50-59, 60-69, ≥70) and adjusted by sex, country, tumor subsite, and tumor purity. Firth’s bias-reduced logistic regressions were used for regressions presenting complete or quasi-complete separation. d-e, Box plots indicating the variation in age of onset according to the presence of colibactin mutational signatures (either SBS88, ID18, or both) in all microsatellite stable cases (d) and across tumor subsites (e). Statistically significant differences were evaluated using multivariable linear regression models adjusted by sex, country, tumor purity, and tumor subsite (only for the analysis of all cases, d). The line within the box is plotted at the median, while the upper and lower ends indicate the 25^th and 75^th percentiles. Whiskers show 1.5 × interquartile range, and values outside it are shown as individual data points.

By contrast, enrichment in early-onset cancers was observed for colibactin-induced signatures. Signatures SBS88 and ID18 were 2.5 and 4 times more common, respectively, in colorectal cancers diagnosed below than above the age of 50 (q=0.006 and q=3.7×10^-7, respectively; Fig. 3a-b). The primary associations of early-onset cases with SBS88 and ID18 were further supported by the successive decline in the prevalence of these signatures with increasing age of diagnosis (p-trend=1.3×10^-4 and p-trend=2.0×10^-7, respectively; Fig. 3c; Supplementary Table 20). A similar effect was observed using a complementary motif enrichment analysis for detecting SBS88, similarly to a recent study²⁶ (p-trend=1.0×10^-7; Extended Data Fig. 7a-b). On the basis of the strong co-occurrence of SBS88 and ID18 (p=7.4×10^-63), as well as previous functional³³ and population studies^22,26,34, we defined exposure to colibactin by the presence of either SBS88 or ID18. Colibactin exposure was found in 21.1% of all colorectal cancers (169/802) and was associated with earlier age of onset (median age: 62 vs. 67, p=1.6×10^-8; Fig. 3d), an effect more evident in the distal colon (median age: 57 vs. 66, q=5.2×10^-7) and rectum (median age: 63 vs. 66, q=0.025; Fig. 3e). Overall, colibactin exposure had a strong inverse correlation with age, being 3.3 times more common in colorectal cancers diagnosed in individuals younger than 40 compared to those over 70 (p-trend=2.7×10^-7, Extended Data Figure 7c).

Signatures of unknown etiology SBS_M and ID14 (Fig. 3a-c) were also enriched in early-onset cases, and SBS89 similarly exhibited a higher prevalence in younger individuals (5.8 times more prevalent in early-onset compared to late-onset patients with p-trend=0.047), albeit based on a very small number of cancers harboring the signature (9/802, 1.1%; Fig. 3c). Interestingly, SBS_M showed an elevation in distal colon and rectum tumors compared to proximal colon similar to the one observed in colibactin-associated signatures SBS88 and ID18, previously reported²² (Supplementary Fig. 9).

Colibactin mutagenesis is an early event in colorectal carcinogenesis

To time the imprinting of SBS88 and ID18, mutations were categorized as early clonal, late clonal, or subclonal during the development of each cancer and the contribution of each mutational signature to each category was determined (Methods). SBS88 and ID18 were both enriched in early clonal compared to late clonal mutations (q=4.2×10^-4 and q=6.1×10^-5; Fig. 4a), as well as a similar trend in clonal compared to subclonal mutations (q=0.138 and q=0.058; Extended Data Fig. 8a), consistent with the presence of these mutational signatures in normal colorectal epithelium³⁴. This enrichment in earlier evolutionary stages was similar to the one observed for other well-known clock-like signatures like SBS1, SBS5, or ID1 (Fig. 4a-b), as previously shown in tumors^37,38 and normal tissues³⁴, and in contrast to signatures known to preferentially generate late clonal and subclonal mutations, such as SBS17a/b³⁸. Interestingly, the enrichment of colibactin signatures in early clonal mutations was observed for both early-onset (q=0.004 for SBS88 and q=2.0×10^-4 for ID18) and late-onset colorectal cancer cases (q=0.020 and q=0.024; Extended Data Fig. 8b).

• Download figure
• Open in new tab

Fig. 4. Colibactin mutagenesis as an early event in microsatellite stable colorectal cancer evolution.

a, Box plots indicating the fold-change of the relative contribution per sample of each signature between early clonal and late clonal single base substitutions (SBS, left) and small insertions and deletions (ID, right). SBS signatures that generated early clonal SBSs in fewer than 50 samples and also generated late clonal SBSs in fewer than 50 samples were excluded from the analysis. Similarly, ID signatures that generated early clonal IDs in fewer than 20 samples and late clonal IDs in fewer than 20 samples were also excluded. Signatures were sorted by their median fold-change. The line within the box is plotted at the median, while the upper and lower ends indicate the 25^th and 75^th percentiles. Whiskers show 1.5 × interquartile range, and values outside it are shown as individual data points. b, Bar plot indicating the lack of concordance between colibactin exposure status determined by the presence of colibactin-induced mutational signatures SBS88 or ID18, and the microbiome pks status. Statistical significance was evaluated using a multivariable Firth’s bias-reduced logistic regression model (due to quasi-complete separation) adjusted by age of diagnosis, sex, country, tumor subsite, and tumor purity. c-d, Distribution of the age of onset (c) and cases across age groups (d) based on the detection of colibactin-positive samples using genomic and microbiome status. The genomic status is defined by the presence of SBS88 or ID18 signatures, while the microbiome status (pks) is determined by coverage of at least half of the pks island, and suggests ongoing or active pks+ bacterial infection. Statistical significance was evaluated using a multivariable linear regression model adjusted by sex, country, tumor subsite, and tumor purity.

Since colibactin is produced by bacteria carrying the pks pathogenicity island, we investigated whether colorectal cancer cases with SBS88 or ID18 harbored pks+ bacteria based on sequencing reads from the cancer sample that did not map to the human genome but mapped to the pks locus (Methods). Consistent with a prior observation³⁹, there was no association between the presence of SBS88 or ID18 and that of pks+ bacteria (Fig. 4b; Extended Data Fig. 9). Similarly, no microbiome association was observed for the other signatures enriched in early-onset colorectal cancers (Supplementary Note). Moreover, we observed a younger age of diagnosis for cases with SBS88 or ID18 but without an identified pks+ bacteria (p=1.3×10^-7; Fig. 4c-d). While the reasons are unclear, one likely explanation is the imprinting of SBS88 and ID18 on the colorectal epithelium during an early period of life when pks+ bacteria were present, followed by the natural plasticity of the microbiome over subsequent decades, leading to the loss or gain of pks+ bacteria.

Colibactin exposure and driver mutations

Using the IntOGen framework⁴⁰, 46 genes under positive selection were identified, with eight mutated in more than 10% of cancers: APC, TP53, KRAS, FBXW7, SMAD4, PIK3CA, TCFL2, and SOX9 (Fig. 5a; Supplementary Table 21). Forty-three of the 46 genes have been previously reported as colorectal cancer driver genes^22,40, two in other cancer types (MED12, NCOR1)⁴⁰, and a putative novel colorectal cancer driver gene, CCR4, was identified with mutations indicating inactivation of the encoded protein. Mutations affecting these 46 cancer driver genes were annotated as driver mutations using a multi-step process based on the mutation type and the mode of action of the gene (Methods). An elevation in the total number of driver mutations was observed in late-onset compared to early-onset cases (FC=1.21, p=5.4×10^-5; Fig. 5b). In addition, an enrichment in APC driver mutation carriers was also found for late-onset cases (OR=2.7, q=0.027; Fig. 5c-d; Supplementary Table 22), as previously reported⁴¹, whereas no hotspot driver mutation (defined as those affecting the same genomic position in at least 10 cases) was associated with age of onset (q>0.05; Supplementary Table 23). No statistically significant differences across countries were found for driver mutations within cancer driver genes or for hotspot driver mutations (q>0.05; Supplementary Tables 24-25).

• Download figure
• Open in new tab

Fig. 5. Variation of driver mutations with age of onset and association with colibactin mutagenesis in microsatellite stable colorectal cancers.

a, Bar plot indicating the prevalence of driver mutations affecting the 48 bioinformatically detected driver genes in microsatellite stable colorectal cancers. Genes were colored according to their status as known cancer driver genes for colorectal cancer, known cancer driver genes for other cancer types, or newly detected cancer driver genes. b, Box plots indicating the distribution of total driver mutations across early-onset (under 50 years of age; purple) and late-onset (50 or over; green) tumors. Statistical significance was evaluated using a multivariable linear regression model adjusted by sex, country, tumor subsite, and tumor purity. The line within the box is plotted at the median, while the upper and lower ends indicate the 25^th and 75^th percentiles. Whiskers show 1.5 × interquartile range, and values outside it are shown as individual data points. c, Volcano plot indicating the enrichment of driver mutations in cancer driver genes in early-onset and late-onset cases. Statistically significant enrichments were evaluated using multivariable logistic regression models adjusted by sex, country, tumor subsite, and tumor purity. Firth’s bias-reduced logistic regressions were used for regressions presenting complete or quasi-complete separation. P-values were adjusted for multiple comparisons based on the total number of cancer driver genes considered and reported as q-values. Horizontal lines marking statistically significant thresholds were included at 0.05 (dashed orange line) and 0.01 q-values (dashed red line). d, Line plot indicating the prevalence of driver mutations in cancer driver genes across ages of onset, using five different age groups. Cancer driver genes significantly enriched in late-onset cases (as shown in c) were colored in green, whereas genes not varying significantly with age of onset were colored in grey. e, Bar plots indicating the proportion and number of driver mutations probabilistically assigned to colibactin-induced and other mutational signatures, including single base substitutions (e) and small insertions and deletions (indels; f). Driver mutations were divided into different groups, including the APC c.835-8A>G splicing-associated driver mutation, other APC driver mutations, TP53 driver mutations, and driver mutations affecting other cancer driver genes.

The contributions of SBS88 and ID18 to driver mutations were assessed using probabilistic assignment of signatures to individual mutations⁴². SBS88 accounted for 64.3% of the colibactin-induced⁴³ APC splicing variant c.835-8A>G in colibactin-exposed samples, compared to only 3.9% and 3.8% of driver substitutions in APC or other cancer genes (Fig. 5e). Similarly, ID18 accounted for 25.3% of APC driver indels and 16.9% of other driver indels in colibactin-exposed cases (Fig. 5f). Overall, SBS88 and ID18 accounted for 8.3% of all SBS and ID driver mutations, and 15.5% of all APC driver mutations in colibactin positive cancers. Nevertheless, no differences were observed between early-onset and late-onset colibactin positive colorectal cancer in the proportion of driver mutations assigned to specific mutational signatures (Extended Data Fig. 10a-b). In addition, a prior study observed that SBS88 is also responsible for mutations in chromatin modifier genes³⁹, and we were able to validate this as well as show a similar effect for the colibactin-associated indel signature, ID18 (Extended Data Fig. 10c-d). Interestingly, using a similar methodology, we observed an elevated number of driver mutations assigned to SBS94 and SBS_F in Colombia, as well as SBS89 and ID_J in Argentina, compared to other countries (Extended Data Fig. 10e-g).

Recruitment of patients and informed consent

The International Agency for Research on Cancer (IARC/WHO) coordinated case recruitment through an international network of 17 collaborators from 11 participating countries in North America, South America, Asia, and Europe (Supplementary Table 26). The inclusion criteria for patients were ≥18 years of age (ranging from 18 to 95, with a mean of 64 and a standard deviation of 12), confirmed diagnosis of primary colorectal cancer, and no prior treatment for colorectal cancer. Informed consent was obtained for all participants. Patients were excluded if they had any condition that could interfere with their ability to provide informed consent or if there were no means of obtaining adequate tissues or associated data as per the protocol requirements. Ethical approvals were first obtained from each Local Research Ethics Committee and Federal Ethics Committee when applicable, as well as from the IARC/WHO Ethics Committee.

Bio-samples and data collection

Dedicated standard operating procedures, following guidelines from the International Cancer Genome Consortium (ICGC), were designed by IARC/WHO to select appropriate case series with complete biological samples and exposure information⁴⁶, as described previously^5,6,8 (Supplementary Table 26). In brief, for all case series included, anthropometric measures were taken, together with relevant information regarding medical and familial history. All biological samples from retrospective cohorts were collected using rigorous, standardized protocols and fulfilled the required standards of sample collection defined by the IARC/WHO for sequencing and analysis. Potential limitations of using retrospective clinical data collected using different protocols from different populations were addressed by a central data harmonization to ensure a comparable group of exposure variables (Supplementary Table 26). All patient-related data were pseudonymized locally through the use of a dedicated alpha-numerical identifier system before being transferred to the IARC/WHO central database.

Expert pathology review

Original diagnostic pathology departments provided diagnostic histological details of contributing cases through standard abstract forms, together with a representative hematoxylin-eosin-stained slide of formalin-fixed paraffin-embedded tumor tissues whenever possible. IARC/WHO centralized the entire pathology workflow and coordinated a centralized digital pathology examination of the frozen tumor tissues collected for the study as well as formalin-fixed paraffin-embedded sections when available, via a web-based approach and dedicated expert panel following standardized procedures as described previously^5,6. A minimum of 50% viable tumor cells was required for eligibility for whole-genome sequencing. In summary, frozen tumor tissues were first examined to confirm the morphological type and the percentage of viable tumor cells. A random selection of tumor tissues was independently evaluated by a second pathologist. Enrichment of tumor component was performed by dissection of the non-tumoral part, if necessary.

DNA extraction

A total of 1,977 primary colorectal cancer patients were enrolled into the study, including biological samples for 1,946 cases and sequencing data (FASTQ) for 31 cases from Japan. Of these, 906 samples (45.8%) were excluded due to insufficient viable tumor cells (pathology level) or inadequate DNA (tumor or germline). Extraction of DNA from fresh frozen primary tumor and matched blood/normal tissue samples was centrally conducted at IARC/WHO (except for samples from Japan) following a similarly standardized DNA extraction procedure. Germline DNA was extracted from whole blood (n=1,015), except for a small subset of Canadian cases (n=25) where only adjacent normal tissue was available, following previously described protocols and methods^5,6. As a result, DNA from 1,040 cases was sent to the Wellcome Sanger Institute for whole-genome sequencing.

Whole-genome sequencing

Fluidigm SNP genotyping with a custom panel was performed to ensure that each pair of tumor and matched normal samples originated from the same individual. Whole-genome sequencing (150bp paired-end) was performed on the Illumina NovaSeq 6000 platform with a target coverage of 40x for tumors and 20x for matched normal tissues. All sequencing reads were aligned to the GRCh38 human reference genome using the Burrows-Wheeler Aligner MEM (BWA-MEM; v0.7.16a and v0.7.17)⁴⁷. Post-sequencing quality control metrics were applied for total coverage, evenness of coverage, contamination, and total number of somatic single base substitutions (SBSs). Cases were excluded if coverage was below 30x for tumor or 15x for normal tissue. For evenness of coverage, the median over mean coverage (MoM) score was calculated. Tumors with MoM scores outside the range of values determined by previous studies⁴⁸ to be appropriate for whole-genome sequencing (0.92-1.09) were excluded. Conpair⁴⁹ was used to detect contamination, cases were excluded if the result was greater than 3%⁴⁸. Lastly, samples with <1,000 total somatic SBSs were also excluded. A total of 981 pairs of colorectal cancer and matched-normal tissue passed all criteria. Comparing the clinicopathological characteristics between the included and excluded patients revealed very similar traits (Supplementary Table 27), and comparable to those expected for each country according to GLOBOCAN metrics (obtained from https://gco.iarc.who.int/today/en/dataviz/; Supplementary Fig. 11).

Germline variant calling

Germline SNVs and indels were derived from whole-genome sequencing from the normal paired material for each individual using Strelka2 with appropriate quality-control criteria⁵⁰. Variant calls were then derived into genotypes for each individual and annotated using ANNOVAR⁵¹.

Somatic variant calling

Variant calling was performed using the standard Sanger bioinformatics analysis pipeline (https://github.com/cancerit). Copy number profiles were determined using ASCAT⁵² and BATTENBERG⁵³ when tumor purity allowed. SNVs were called with cgpCaVEMan⁵⁴, indels were called with cgpPINDEL⁵⁵, and structural rearrangements were called using BRASS (https://github.com/cancerit/BRASS). CaVEMan and BRASS were run using the copy number profile and purity values determined from ASCAT when possible (complete pipeline, n=916). When tumor purity was insufficient to determine an accurate copy number profile (partial pipeline, n=31) CaVEMan and BRASS were run using copy number defaults and an estimate of purity obtained from ASCAT. Finally, for a subset of cases which had no large copy number alterations (copy number normal pipeline, n=34), CaVEMan and BRASS were run using copy number defaults and an estimate of purity calculated by the median variant allele frequency (VAF) of indels multiplied by two. For SNVs, additional filters (ASRD≥140 and CLPM=0) were applied in addition to the standard PASS filter to remove potential false positive calls. To further exclude the possibility of caller-specific artifacts being included in the analysis, a second variant caller, Strelka2⁵⁰, was run for SNVs and indels. Only variants called by both the Sanger variant calling pipeline and Strelka2 were included in subsequent analysis.

Generation of mutational matrices

Mutational matrices for single base substitutions (SBS), indels (ID), doublet base substitutions (DBS), copy number alterations (CN), and structural variants (SV) were generated using SigProfilerMatrixGenerator with default options (v1.2.0)^56,57.

Microsatellite instability validation

The presence of microsatellite instability (MSI) in colorectal cancers was validated using the QX200 Droplet Digital PCR System (Bio-Rad, Hercules, CA, USA) for the detection of five microsatellite markers (BAT25, BAT26, NR21, NR24, and Mono27) commercially pooled in three primer–probe mix assays, as previously described⁵⁸. Briefly, samples were tested in duplicate, and each reaction comprised 1× ddPCR Multiplex Supermix for probes (Bio-Rad), 1X primer–probe mix, and 10 ng of extracted tumor DNA, in a total volume of 22 μl. MSI-positive, negative, and no-template (nuclease-free water) controls were included in each experiment. Droplet generation and plate preparation for thermal cycling amplification were performed using the QX200 AutoDG Droplet Digital PCR System (Bio-Rad). The following PCR protocol was applied on a C1000 Touch Thermal Cycler (Bio-Rad): 37 °C for 30 min, 95 °C for 10 min, followed by 40 cycles of denaturation at 94 °C for 30 s, annealing at 55 °C for 1 min, with a final extension at 98 °C for 10 min. Following PCR amplification, fluorescence signals were quantified using the QX200 Droplet Reader (Bio-Rad), and data were analyzed with QuantaSoft Analysis Pro v.1.0.596.0525 (Bio-Rad) software. Positive and negative controls served as guides to call markers and delineate clusters. For each assay, the cluster at the bottom left of the x–y plot was designated as the negative population. Clusters located vertically and horizontally from the negative cluster were identified as the mutant population, while clusters located diagonally from the negative cluster represented the wild-type population. Tumors were characterized for the MSI phenotype by analyzing the results for all five markers using the following criteria: MSI positive if two or more mutant microsatellite markers were observed, and microsatellite stable (MSS; i.e., MSI negative) when none or only one of the microsatellite markers was altered (Supplementary Note).

Extraction and decomposition of mutational signatures

Mutational signatures were primarily extracted using SigProfilerExtractor³⁵, based on nonnegative matrix factorization, and validated by mSigHdp⁵⁹, based on hierarchical Dirichlet process mixture models. For SigProfilerExtractor (v1.1.21), de novo mutational signatures were extracted from SBS, DBS, and ID mutational matrices using 500 NMF replicates (nmf_replicates=500), nndsvd_min initialization (nmf_init=”nndvsd_min”), and default parameters. Extractions were performed separately on the subsets of 802 MSS and 153 MSI cases (Supplementary Note). De novo SBS mutational signatures were extracted for both SBS-288 and SBS-1536 contexts, which, beyond the common SBS-96 trinucleotide context using the mutated base and the 5’ and 3’ adjacent nucleotides^57,60, also consider the transcriptional strand bias and the pentanucleotide context (two 5’ and 3’ adjacent nucleotides), respectively⁵⁷. The results were largely concordant, with the SBS-288 de novo signatures allowing additional separation of mutational processes. The SBS-1536 results can be found in the Supplementary Note. Therefore, the SBS-288 de novo signatures were taken forward for further analysis (Supplementary Table 5). Previously established mutational contexts DBS-78 and ID-83^19,57 were used for the extraction of DBS and ID signatures (Supplementary Tables 6-7).

Copy number signatures were extracted de novo using SigProfilerExtractor with default parameters and following an updated context definition benefitting from WGS data (CN-68) (Supplementary Table 8), which allowed to further characterize CN segments below 100kbp in length (in contrast to current COSMICv3.4 reference signatures using the CN-48 context, which were based on SNP6 microarray data and therefore without the resolution to characterize short CN segments)⁶¹. SV signatures were extracted using a similarly refined context, with an in-depth characterization of short SV alterations below 1kbp (SV-38 context, in contrast to current COSMICv3.4 signatures based on the SV-32 context⁶²; Supplementary Table 9).

After de novo extraction was completed, SigProfilerAssignment⁴² v0.0.29 was used to decompose the de novo extracted SBS, ID, DBS, CN, and SV mutational signatures into COSMICv3.4 reference signatures based on the GRCh38 reference genome⁶³ (Supplementary Table 10). When possible, SigProfilerAssignment matched each de novo extracted mutational signature to a set of previously identified COSMICv3.4 signatures (Supplementary Note). For the SBS-288, CN-68, and SV-38 signatures, this required collapsing the high-definition classifications into the standard SBS-96, CN-48, and SV-32 mutational classifications, respectively. Four of the de novo extracted MSS SBS signatures did not match any previous COSMICv3.4 signatures, with three of them (SBS_F, SBS_H, and SBS_M) showing a strong similarity with previously reported signatures in the UK population²¹ (cosine similarity > 0.93), and one (SBS_O) reflecting a cleaner version of a previously reported COSMICv3.4 signature (SBS41). To validate the latter, we performed a decomposition of the current mutational profile of signature SBS41 using the decomposed signatures from our analysis, obtaining a confirmation that SBS41 can be reconstructed by a linear combination of SBS_O (contributing 19.00% of the mutational profile), SBS93 (62.54%), and SBS34 (12.60%), and SBS5 (5.86%) with a cosine similarity of 0.91. Notably, SBS93, first identified in gastric tumors³⁵, was unknown at the time SBS41 was first reported¹⁹. For the MSS cohort, one ID (ID_J), one CN (CN_F), and two SV signatures (SV_B and SV_D) were additionally not decomposed into previously known signatures, and therefore considered as novel (Supplementary Table 10). The novel SV signature SV_D, identified in the MSS cohort, was also considered for the decomposition of de novo SV signatures extracted in the MSI cohort. In the MSI cohort, four of the de novo extracted SBS signatures (SBS_I_MSI, SBS_M_MSI, M, SBS_N_MSI, and SBS_O_MSI) as well as one de novo DBS signature (DBS_B_MSI) did not match COSMICv3.4 signatures, with SBS_M_MSI showing a strong similarity with a previously reported signature in the UK population²¹ (cosine similarity=0.89), and the other four signatures considered as novel (Supplementary Note).

mSigHdp⁵⁹ extraction of SBS-96 and ID-83 signatures was performed on the 802 MSS subset using the suggested parameters and using the country of origin to construct the hierarchy. SigProfilerAssignment was subsequently used to match mSigHdp de novo signatures to previously identified COSMIC signatures. A comparison of the signatures extracted from mSigHdp and SigProfilerExtractor can be found in the Supplementary Note.

Attribution of mutational signatures to individual samples

Known COSMIC signatures and de novo signatures that were not decomposed into COSMIC signatures (Supplementary Table 10; Supplementary Note) were attributed for each sample using MSA⁶⁴ (v2.0) for SBS, ID, and DBS, whereas SigProfilerAssignment⁴² was used for CN and SV. A conservative approach was used for MSA attributions utilizing the (params.no_CI_for_penalties=False) option for the calculation of optimum penalties. Pruned attributions were used for the final analysis, where confidence intervals were applied to each attributed mutational signature and any signature activity with a lower confidence limit equal to 0 was removed.

Attribution of mutational signatures to individual somatic mutations

SBS and ID mutational signatures were probabilistically attributed to individual somatic mutations using the MSA activities per sample, based on Bayes’ rule and the specific mutational context for the mutation, as previously described⁴². Briefly, to calculate the probability of a specific mutational signature being responsible for a mutation in a given mutational context and in a particular sample, we multiplied the general probability of the signature causing mutations in a specific mutational context (obtained from the mutational signature profile) by the activity of the signature in the sample (obtained from the signature activities), and then normalized this value dividing by the total number of mutations corresponding to the specific mutational context (obtained from the reconstructed mutational profile of the sample). The signature with the maximum likelihood estimation was assigned to each individual somatic mutation.

Driver gene analysis

Consensus de novo driver gene identification was performed by IntOGen⁴⁰, which combines seven state-of-the-art computational methods to detect signals of positive selection across the cohort. The genes identified as drivers with a combination q-value<0.10 were classified according to their mode of action in tumorigenesis (i.e., tumor suppressor genes or oncogenes) based on the relationship between the excess of observed nonsynonymous and truncating mutations computed by dNdScv⁶⁵ and their annotations in the Cancer Gene Census⁶⁶.

To identify potential driver mutations, we selected SBS or ID mutations that fulfilled any of the following criteria: mutations classified as “Oncogenic” or “Likely Oncogenic” by OncoKB⁶⁷; mutations classified as drivers in the TCGA MC3 drivers study⁶⁸; truncating mutations in driver genes annotated as tumor suppressors; recurrent missense mutations (seen in at least three cases); mutations classified as “Likely Drivers” by boostDM (score >0.50)⁶⁹; or missense mutations classified as “Likely Pathogenic” by AlphaMissense⁷⁰ in driver genes annotated as tumor suppressors. Six of the IntOGen-identified driver genes did not carry any potential driver mutations according to our strict criteria and were therefore excluded from subsequent analysis. In summary, 60 driver genes were identified (46 and 31 for MSS and MSI cases, respectively; Supplementary Table 21; Supplementary Note).

Evolutionary analysis

DPClust⁵³ was run on all complete pipeline MSS samples with Battenberg data (n=774) to identify clonal structure in each sample. The DPClust output was used in running MutationTimeR³⁸ to annotate somatic mutations as early clonal, late clonal, subclonal, or NA clonal. Samples with at least 256 early clonal and late clonal SBSs or 100 early clonal and late clonal IDs were retained and split into separate VCF files (n=574 for SBS; n=430 for ID). MSA⁶⁴ was run on the resulting VCF files to identify the active mutational signatures in the early clonal and late clonal mutations. SBS signatures that were found to generate early clonal SBSs in fewer than 50 samples and also generated late clonal SBSs in fewer than 50 samples were excluded from the analysis. Similarly, ID signatures generating early clonal IDs in fewer than 20 samples and late clonal IDs in fewer than 20 samples were also excluded. Wilcoxon signed-rank tests were used to assess the differences in the relative activity of each signature between the early clonal and late clonal mutations. P-values were adjusted across signatures using the Benjamini-Hochberg method⁷¹, and adjusted p-values were reported as q-values. This process was repeated with the same thresholds for SBSs and IDs to also assess the difference in the relative activity of each signature between clonal and subclonal mutations (n=133 for SBS; n=64 for ID). Due to the lower numbers, signatures that were found to generate clonal somatic mutations in fewer than 10 samples and also generated subclonal somatic mutations in fewer than 10 samples were excluded from the analysis.

Motif analysis

MutaGene⁷² was used to find the number of mutations with the WAWW[T>N]W motif, previously associated with colibactin mutagenesis³³, in each sample, regardless of the DNA strand. This value was then divided by the total number of W[T>N]W mutations per sample to identify the percentage of W[T>N]W mutations with the colibactin mutational motif.

Microbiome analysis

To identify microbial reads that map to the pks island (pks), non-human reads were aligned to the IHE3034 genome (RefSeq assembly: GCF_000025745.1) using Bowtie2⁷³. IHE3034 is a pks E. coli strain that contains the pks island with all 19 clb genes in the clbA-clbS gene cluster. Prior to alignment, poor quality reads were filtered using fastp⁷⁴, and the remaining human reads were removed by excluding those that mapped to GRCh38, T2T-CHM13v2.0, and the 47 pangenomes⁷⁵. A sample was considered pks+ if it had at least one read across at least 8 out of the 19 genes in the clbA-clbS gene cluster. Genome coverage circos plots were generated using Reads Per Kilobase per Million (RPKM) values and visualized with the circlize R package⁷⁶.

Regressions

To compare the mutation burden of different variant types, a linear regression of the mutation burden logarithm (base 10) was considered, using age, sex, tumor subsite, country, and tumor purity as independent variables. For mutational signature-based analyses, signature attributions were dichotomized into presence and absence using confidence intervals, with presence defined as both lower and upper limits being positive and absence as the lower limit being zero. If a signature was present in at least 70% of cases (SBS1, SBS5, SBS18, ID1, ID2, ID14 and CN2 for MSS cases; ID1, ID2, DBS_B_MSI, CN1, and SV_D for MSI cases), it was dichotomized into above and below the median of attributed mutation counts. The binary attributions served as dependent variables in logistic regressions. Regressions with variables presenting complete or quasi-complete separation⁷⁷ were performed using Firth’s bias-reduced logistic regressions based on the logistf R package. To adjust for confounding factors, sex, age of diagnosis, tumor subsite, country, and tumor purity were added as covariates in all regressions, serving as independent variables for the regressions. The tumor subsite variable was categorized as proximal colon (ICD-10-CM codes C18.0, C18.2, C18.3, and C18.4), distal colon (C18.5, C18.6, and C18.7), or rectum (C19 and C20), unless otherwise specified. One MSI tumor from an unspecified subsite was removed for the multivariable regression models in MSI cases. The age of diagnosis variable was generally considered as a numerical variable, or categorized into two (early-onset, <50 years old; and late-onset, ≥50) or five subgroups (0-39, 40-49, 50-59, 60-69, ≥70), depending on the analysis performed, with specific indications in the corresponding figure legends. Similarly, regressions for driver mutations in cancer driver genes and hotspot driver mutations (present in at least 10 cases) were done using the same logistic regression models but replacing signature by driver mutation prevalence across samples.

Regressions with colorectal cancer incidence were performed as linear regressions with signature attributions with confidence intervals not consistent with zero as dependent variables, and age-standardized rates (ASR) of colorectal cancer (and independent ASR of colon and rectal cancer) obtained from the Global Cancer Observatory (GLOBOCAN)¹, sex, age of diagnosis, tumor subsite, and tumor purity as independent variables. Regressions were performed on a sample basis.

Regressions with colibactin presence (based on genomic and/or microbiome-derived detection) were performed as linear regressions with age of diagnosis as the dependent variable, and sex, tumor subsite, country, and tumor purity as independent variables.

Additional statistical analyses

For regressions of signatures, driver mutations in cancer driver genes, and hotspot driver mutations, p-values were adjusted for multiple comparisons based on the total number of decomposed reference mutational signatures considered per variant type (i.e., 19 SBS, 7 DBS, 11 ID, 9 CN, and 11 SV signatures for MSS cases; 18 SBS, 10 DBS, 2 ID, 4 CN, and 4 SV for MSI cases), cancer genes (46 for MSS; 31 for MSI), or hotspot driver mutations (38 for MSS; 14 for MSI) using the Benjamini-Hochberg method⁷¹. For country enrichment analyses, the mutation burdens and binary attributions of mutational signatures were compared for each country against all others. Therefore, p-values were also adjusted for multiple comparisons based on the total number of countries assessed (a total of 11 countries). Adjusted p-values were reported as q-values, with q-values<0.05 considered statistically significant. For age of diagnosis-based regressions of colibactin presence across tumor subsites, p-values were adjusted and reported as q-values based on the total number of tumor subsites assessed (a total of 3 tumor subsites). For the age of diagnosis trend enrichment analysis of signatures, p-trends were reported, with p-trends<0.05 considered statistically significant. For evidence of co-occurrence or mutual exclusivity of two signatures, two-sided Fisher’s exact tests were used, and p-values were reported, with p-values<0.05 considered statistically significant.