Wnt deprivation selects AXIN1-mutant and AXIN2-mutant organoids

On the basis of our previous study, we knew that MLH1KO human colon organoids recapitulate the MMR-deficient phenotype in vitro and accumulate substantial numbers of substitutions and insertion or deletions (indels) during in vitro culture17. Withdrawal of Wnt3a and R-spondin 1 from the culture medium allowed us previously to select CRISPR-induced APC inactivation mutant organoids13,22. Here, we attempted to similarly select Wnt signal-independent MLH1KO human colon organoids, yet without specific gene targeting (Fig. 1a). As expected, WT and MLH1KO organoids rapidly died after single-cell seeding in a medium lacking Wnt3a and R-spondin 1 (−WR). This allowed us to recover rare MLH1KO subclonal cells (MLH1KO-1) from an MLH1KO organoid line that was continuously cultured for ~350 days after single-cell cloning. No such subclones could be derived from two other MLH1 mutation clones derived from the same donor (MLH1KO-2 and MLH1KO-3) (Fig. 1b,c). No Wnt-independent subclones could be isolated from the MLH1KO-1 line at earlier time points (30, 100 and 200 days). Single surviving organoids were manually picked and expanded for downstream analysis (Fig. 1d). These organoid clones (MLH1KO-1−WR-1 and MLH1KO-1−WR-2) presented as dense and budding structures when cultured in −WR medium, whereas APCKO organoids are reported to display a much more cystic structure23,24. Similar selections were performed using media designed to inhibit or activate key pathways of the stem cell niche (EGF receptor (EGFR), BMP and p53). Despite the limited initial clonal outgrowth, we could not establish long-term expanding cultures under these conditions, suggesting that no mutations affecting these pathways were present at this stage (Extended Data Fig. 1a,b). The emergence of surviving clones in the −WR condition was also observed in the MLH1KO clonal line derived from another donor, which indicated this to be a reproducible phenomenon (Extended Data Fig. 1c,d).

Fig. 1: Niche factor deprivation strategy successfully selected AXIN1-mutant and AXIN2-mutant organoids from MLH1KO human colon organoids.

a, Strategy to select the survived organoids with deprivation of Wnt pathway factors. Organoids were dissociated into single cells and then 200,000 single cells were cultured in niche-factor-deprived medium to select mutated stem cells (red) out from WT stem cells (blue) and differentiated cells (white). W, Wnt3a; R, R-spondin 1; N, Noggin; E, EGF. Surviving clone organoids were expanded briefly for WGS and WES analysis. Accumulated mutations were determined by subtracting the mutations in the original clone (blue) from the mutations in the selected subclone (red). b, WT organoids and MLH1KO-1 organoids grew in complete medium (WRNE) but some of organoids grew out only from MLH1KO-1 organoids in −WR medium (representative pictures from n = 3 clonal organoid lines). Scale bars, 500 µm. (c) Quantification of survived organoids in −WR medium. Mean and s.d. (error bars) of n = 4 clonal organoid lines. P values (two-sided Welch’s t-test): MLH1KO-1 versus MLH1KO-2, P = 0.0001; MLH1KO-1 versus MLH1KO-3, P = 0.0002. d, Single survived organoids (MLH1KO-1−WR-1 and MLH1KO-1−WR-2) were expanded and passaged in −WR medium (representative pictures from n = 3 independent replicates). Scale bars, 500 µm. e, The Venn diagrams illustrate the number of base substitutions and indels found in WGS data of survived clones (MLH1KO-1−WR-1 and MLH1KO-1−WR-2). AXIN1 and AXIN2 candidate gene mutations are listed in the overlapping section. The numbers indicate nonsynonymous mutations and the numbers in brackets indicate all mutations found in WGS. f, RT–qPCR analysis of AXIN2 mRNA levels. AXIN2 transcripts were compared between each clone cultured in the medium with or without Wnt3a and R-spondin 1. Mean and s.d. (error bars) of n = 3 biologically independent experiments. P values (two-sided Welch’s t-test): MLH1KO-1 versus MLH1KO-1−WR-1 in −WR medium, P = 0.0048; MLH1KO-1 versus MLH1KO-1−WR-2 in −WR medium, P = 0.0037. g, Representative images of Sanger DNA sequencing of AXIN1 and AXIN2 mutations in Wnt-independent clones (MLH1KO-1−WR-1 and MLH1KO-1−WR-2).

Source data

The selected subclones (MLH1KO-1−WR-1 and MLH1KO-1−WR-2) were subjected to whole-genome sequencing (WGS) analyses. The genomic sequence of the parental MLH1KO-1 clone allowed us to subtract somatic mutations that were already present in the parental clone (Fig. 1a). As expected from our previous study, surviving clones accumulated single-base substitutions (SBSs) (8,748 ± 149 mutations) and indels (23,010 ± 110 mutations), which recapitulated the mutator phenotype of MSI-H tumors with high tumor mutational burden (TMB) (Fig. 1e). The number of actual mutations and the previously calculated mutation number per genome division allowed us to calculate an estimated culture duration (354 ± 69 days)17. This matched our actual culture duration (417 days) well. The clonality of selected organoid clones was confirmed by variant allele frequency analyses (Extended Data Fig. 2). We compared the mutations acquired by the two surviving clones to the parental MLH1KO clone. Among 6,372 base substitutions and 18,479 indels shared in two clones, we identified 74 base substitutions and 81 indels as nonsynonymous somatic coding mutations (Supplementary Table 1). Among those listed mutations, we found two compound heterozygous AXIN1 mutations (c.1051C>T, p.Gln351* and c.1523 del, p.Gly508fs*197) and a homozygous AXIN2 mutation (c.1994 del, p.Gly665fs*) carried in both subclones (Fig. 1e). These mutations are reported as hotspot mutations in the cBioPortal and in the Catalog of Somatic Mutations in Cancer (COSMIC) database as cancer-related somatic variants (COSM8405446, COSM1240884 and COSM1385326) (Extended Data Fig. 3). Biallelic loss-of-function mutations of both AXIN1 and AXIN2 lead to Wnt pathway activation and occur in CRC25. Indeed, a recent report confirmed that CRISPR KO of AXIN1 of human intestinal organoids allows growth in the absence of Wnt and R-spondin 1, similar to APCKO organoids24. In contrast, the common co-occurrence of AXIN2 mutations with AXIN1 and other Wnt pathway mutations suggests that weak pathway activation caused by AXIN2 loss of function may synergistically upregulate Wnt signal26,27,28.

We performed a medium selection experiment comparing the AXIN-mutant subclones with organoids that were CRISPR-engineered to carry mutations in other driver genes. This confirmed that organoids with AXIN1 and AXIN2 mutations upon MLH1KO acquire the same Wnt signal independency as APC-mutant organoids (Extended Data Fig. 4). Wnt pathway activation was further confirmed by reverse transcription (RT)–qPCR analysis for the Wnt target gene AXIN2, as AXIN2 was expressed at higher levels in the selected clones (MLH1KO-1−WR-1 and MLH1KO-1−WR-2) than in the WT and MLH1KO clone without Wnt3a and R-spondin 1 (Fig. 1f). The SBS in AXIN1 p.Gln351* (GCA>GTA) matched the MMR-linked mutational signature SBS 20 observed in MLH1KO organoid17,29. Both AXIN1 p.Gly508fs* and AXIN2 p.Gly655fs* result from single-nucleotide deletion in homopolymer cytosine repeats, representative of the mutation pattern known as slippage in MMR-deficient and MSI-H tumors (Fig. 1g)8,30. Importantly, these types of AXIN1 and AXIN2 mutations are reported in MSI-H CRCs and gastric cancers7,8,9,30,31. Thus, MLH1KO organoids allowed the faithful selection of some of the specific mutations seen in MMR-deficient and MSI-H tumors in in vitro culture.

Sequential mutant selections gain niche independence

Because these observations indicated that functional selection of spontaneous oncogenic mutations in MLH1KO organoids allowed us to model MMR-deficient tumorigenesis in vitro, we performed sequential selection of these mutant clones (MLH1KO-1−WR-1 and MLH1KO-1−WR-2) in −Noggin, −EGF +gefitinib (−EGF+Gef) or +Nut3 medium. We, thus, hoped to find BMP pathway mutations, Ras and Raf pathway mutations and TP53 mutations, respectively. Surviving organoids readily grew out from the MLH1KO-1−WR-1 clone in the −Noggin condition. In contrast, only one organoid survived and grew into a large budding organoid from the MLH1KO-1−WR-1 clone in +Nut3 medium (Fig. 2a,b). These surviving organoids were picked and clonally expanded for further genomic analyses (Extended Data Fig. 5). As expected, a missense mutation in TP53 (c.473G>T, p.Arg158Leu (Hom)) was found in the Nut3-selected clone (MLH1KO-1−WR-1+Nu) by whole-exome sequencing (WES), which was confirmed by Sanger sequencing. Aberrant p53 protein expression and increased p21 expression were confirmed by western blot (Fig. 2c). This TP53 mutation is mostly reported in lung tumors but also seen in persons with CRC (Extended Data Fig. 3)32,33.

Fig. 2: Further combinational mutant selection recapitulates sequential independence from niche factors.

a, MLH1KO-1−WR-1 organoid grew in −WR medium but some of organoids grew out in −WR−Noggin medium (representative pictures from n = 3 clonal organoid lines) and a single organoid grew out in −WR+Nut3 medium. Scale bars, 500 µm. All cells died out in −WR−EGF+Gef medium. b, Quantification of survived organoids in each selection medium. Mean and s.d. (error bars) of n = 4 clonal organoid lines. P values (two-sided Welch’s t-test): −WR−Noggin versus −WR−EGF+Gef, P < 0.0001; −WR−Noggin versus −WR+Nut3, P < 0.0001; −WR−EGF+Gef versus −WR+Nut3, P = 0.0002. c, Representative western blot analysis to detect p53, p21 and GAPDH in WT and MLH1KO-1 organoid lines with WT p53 (MLH1KO-1 and MLH1KO-1−WR-1) and with mutant p53 (MLH1KO-1−WR-1+Nut3) treated with DMSO or 10 μM Nut3 for 24 h. The analysis was performed with n = 2 biologically independent experiments for one organoid clone. d, MLH1KO-1−WR-1+Nut3 organoids were selected in −WR+Nut3−Noggin or −WR+Nut3−Noggin−EGF+Gef medium (representative pictures from n = 3 clonal organoid lines). Scale bars, 500 µm. e, Quantification of survived organoids in −Noggin and −EGF+Gef medium. Mean and s.d. (error bars) of n = 4 clonal organoid lines. P values (two-sided Welch’s t-test). f, Mutational spectra of all base substitutions observed for each selected clone. g,h, Heat map showing the cosine similarity scores for each indicated clone and COSMIC SBS signature (g) and indel signature (h). Arrows indicate signatures associated with deficiency in DNA MMR in previous analyses29,38.

Source data

Interestingly, the Noggin-independent clones arising from the MLH1KO-1−WR-1 clone (MLH1KO-1−WR-1−N) carried ACVR2A mutations (c.285 del, p.Thr96* (Hom)) or BMPR2 mutations (c.255G>A, p.Trp85* (Het) and c.1748 del, p.Asn583fs* (Het)). SMAD4 mutation, the most common BMP pathway driver mutation in sporadic CRC, was not encountered (Extended Data Fig. 3). Another Noggin-independent clone selected from the MLH1KO-1−WR-2 clone (MLH1KO-1−WR-2−N) also carried partially the same ACVR2A (c.285 del, p.Thr96* (Het) and c.1310 del, p.Lys437fs(Het)) and BMPR2 (c.1748 del, p.Asn583fs* (Hom)) mutations. These mutations occurred through one nucleotide deletion in the polyadenine tract (ACVR2 A8, BMPR2 A11), as reported in MSI-H CRCs and recently defined as a representative indel mutational signature for MSI (ID2 signature)3,10,12,29,34,35. In particular, ACVR2A mutations were recently found and functionally validated as driver mutations for CRC36. The BMPR2 mutation is reported to be mutually exclusive with SMAD4 in sporadic CRC34.

These double-mutant lines were selected again in additional conditions. Some surviving clones arose in −Noggin or −EGF+Gef medium from the TP53 mutated clone (MLH1KO-1−WR-1+Nu-1). Again, surviving organoids readily grew out in −Noggin medium (Fig. 2d,e). We selected two Noggin-independent clones in −Noggin medium and sent these for WES analyses (Extended Data Fig. 6). Interestingly, we confirmed that one (MLH1KO-1−WR-1+Nu−N-1) acquired duplication of a somatic mutation in ACVR2A (c.285 del, p.Thr96*(Hom)) and the other one (MLH1KO-1−WR-1+Nu−N-2) acquired compound heterozygous mutations in ACVR2A (c.285 del, p.Thr96*(Het) and c.1310 del, p.Lys437fs). As the earlier clone (MLH1KO-1−WR-1+Nu-1) harbored a single ACVR2A (c.285 del, p.Thr96*) mutation, this means that the clone obtained another ACVR2A mutation during the selection step, resulting in loss of heterozygosity (LOH). Of note, this clone did not carry the BMPR2 homozygous mutations observed in the previously selected Noggin-independent clones (MLH1KO-1−WR-1−N, MLH1KO-1−WR-2−N), suggesting that the clone acquired Noggin independency solely through LOH of the ACVR2A gene. This was consistent with the previous findings identifying ACVR2A as a driver gene in the transforming growth factor-β (TGFβ) pathway36 and reporting the BMPR2 as contributing cooperatively to the development of sporadic CRC but not as single driver mutation10,34,37. We again extracted mutational signatures from each selected clone. All clones showed the expected mutational spectra (Fig. 2f)17. We next extracted mutational signatures for comparison to the COSMIC database, using cosine similarity as a measure of closeness29,38. All selected organoids showed mutational signatures that resembled signatures SBS6, SBS15, SBS20 and SBS44, which are associated with defective DNA MMR (d-MMR) and occur in microsatellite unstable tumors38 (Fig. 2g). Again, signature ID2 (representing the slippage of DNA MMR deficiency) was found in all mutant clones (Fig. 2h). These data clearly suggest that the mutational signature of d-MMR induced by MLH1 mutation remained the dominant driver of mutagenesis regardless of the accumulated genetic mutations.

Long-term cultured d-MMR organoids complete tumor progression

As early attempts to select out mutant clones carrying EGF, Ras and Raf pathway mutations failed, we extended culturing times of the triple-pathway mutant clone (MLH1KO-1−WR-1+Nu−N-1) up to 300 days to allow for the accumulation of further critical driver mutations. Recently, differential drug sensitivities of EGFR inhibition were reported among various different RAS-mutant patient-derived CRC organoids (PDOs). The pan-human EGFR inhibitor afatinib (Afa) can eliminate ERK activity oscillations in KRAS-mutant PDOs39,40. To select Ras and Raf pathway mutant clones more reliably, we used Afa as an alternative EGFR inhibitor in addition to Gef. Some organoids grew out in −EGF+Gef or −EGF+Afa conditions (Fig. 3a). These surviving organoids could be continuously expanded under both −EGF+Gef and −EGF+Afa conditions (Fig. 3b).

Fig. 3: MMR-deficient organoids acquired quadruple colorectal driver mutations in long-term culture.

a, MLH1KO-1−WR+Nu−N organoids were selected in −WRN+Nu−EGF+Afa/Gef medium (MLH1KO-1−WR+Nu−N−E+Afa/Gef). Outgrowing clones are shown in the inset. Representative pictures from n = 3 independent replicates. Scale bars, 500 µm. b, Survived MLH1KO-1−WRN+Nu−EGF+Afa/Gef organoids were picked up and continuously expanded in each selection medium. Representative pictures from n = 3 clonal organoid lines. Scale bars, 500 µm. c, Sanger DNA sequencing result of NRASQ61K (Het) mutation in the selected clone (MLH1KO-1−WR+Nu−N−E+Afa). d, Representative western blot analysis to detect p44/42 MAPK and β-actin in NRASWT and NRASQ61K organoid lines cultured in WRNE medium or −WRNE medium + 10 μM Nut3 and 1 μM Afa. The analysis was performed with n = 2 biologically independent experiments for one organoid clone. e, Representative H&E staining and IHC analysis of E-cadherin, Ki67, β-catenin and EGFR in the parental clone (NRASWT) and the mutant clone (NRASQ61K). The analysis was performed with n = 2 clonal organoid lines. Scale bar, 100 µm. f, Dose–response curves showing the sensitivity of NRASWT and NRASQ61K clones to inhibitors of the EGFR (Afa), MEK (selumetinib) and SHP2 (SHP099). The y axis represents the viability of the samples while the x axis represents the drug concentration (power of 10 µM). Mean and s.d. (error bars) of n = 3 technical replicates for n = 2 independent organoid clones for each genetic background. g, Heat maps reporting the viability of the drug screening performed on NRASWT and NRASQ61K organoids as normalized to vehicle-treated organoids. Color mapping ranges from low viability (red) to high viability (green). The analysis was performed with n = 3 technical replicates for n = 2 independent organoid clones for each genetic background. h, Timescale phylogenic tree of the selected mutant clones that arose from the MLH1KO-1 single organoid clone. Organoids were selected under each condition (colored rectangles). Circles represent survived clones and dashed-line circles represent nonrecovered clones. Horizontal bars indicate the culture duration of cloning and expanding each clone. All mutations in each clone are listed in Supplementary Table 3.

Source data

Through WES analyses, we found that surviving clonal organoids under both inhibitor conditions (MLH1KO-1−WR-1+Nu−N−E+Afa/Gef) acquired a heterozygous somatic mutation in NRAS (c.181C>A, p.Gln61Lys(Het) COSM580) (Fig. 3c and Extended Data Fig. 6), a known oncogenic mutation in the Ras and Raf pathway3,41. While KRAS mutations are reportedly more common in low-MSI (MSI-L) relative to MSI-H CRC, the NRAS mutation is equally frequent between MSI-L and MSI-H CRC tumors42. To confirm that the selected NRASQ61K mutant organoids (MLH1KO-1−WR-1+Nu−N−E+Afa/Gef) indeed exert constitutive ERK signaling activity, transcription of ERK target genes (ETV4, ETV5, DUSP5, DUSP6 and CCND1) was compared between NRASWT and NRASQ61K mutant organoids by RT–qPCR (Extended Data Fig. 7a). Furthermore, western blot analysis confirmed that NRASQ61K mutant organoids increased phosphorylation of ERK1 and ERK2 (Fig. 3d). Taken together, we obtained quadruple-pathway mutant organoids, which sequentially acquired MSI-H mutations in AXIN1 and AXIN2 (Wnt pathway), TP53, ACVR2A and BMPR2 (BMP pathway) and NRAS (Ras and Raf pathway) and, thus, became completely independent from the pertinent niche factors. Both triple-pathway and quadruple-pathway mutant organoids were highly proliferative and, as indicative of an active Wnt signaling pathway, accumulated expression of nuclear β-catenin (Fig. 3e). The quadruple-pathway mutant organoids appeared as more compact/less cystic structures as seen in our previous CRISPR-engineered quadruple-pathway mutant organoids13. To investigate the acquisition of drug resistance in the quadruple-pathway mutant organoids, we performed drug screening with various pathway inhibitors. We, thus, observed the expected sensitivity of NRASWT organoids toward Afa and Gef and the resistance of NRAS-mutant organoids. In contrast, the mutant organoids remained sensitive to the MEK inhibitors selumetinib and trametinib, the mTOR inhibitor everolimus, the PI3Kα inhibitor alpelisib and the multiple-CDK (cyclin-dependent kinase) inhibitor dinaciclib (Fig. 3f,g and Extended Data Fig. 7b–d). Recently, Ras status was reported to correlate with SH2 domain-containing phosphatase 2 (SHP2) inhibitor sensitivity in several cancers; in particular, the NRASQ61K mutation in neuroblastoma confers resistance to SHP2 inhibitors43,44. Consistent with this, our NRASQ61K mutant organoids exhibited resistance to SHP2 inhibitor SHP099 treatment (Fig. 3f,g). We constructed a mutant phylogenic tree for each sequential mutant clone (Fig. 3h and Supplementary Table 3). The acquisition of identical mutations was independently observed in different clones. The available data indicated that these recurrent mutations arose independently yet were because of a shared mutational mechanism. Taken together, simple selection in an MLH1KO mutant organoid background allowed the sequential isolation of spontaneous quadruple-pathway mutant organoids, starting from a single MLH1KO stem cell, over a period of no more than 30 months.

MSI status and mutational tendency in d-MMR organoids

To assess the MSI profile of our MMR-deficient organoids, we compared the somatic mutations acquired in our quadruple-pathway mutant organoids to those of CRCs from The Cancer Genome Atlas (TCGA) database (TCGA-COAD). Leveraging two bioinformatic MSI predicting tools, MSMuTect and MANTIS45,46,47, we derived MSI scores from WES data. Notably, MANTIS enabled the extraction of TCGA-COAD samples annotated as MSI-H and microsatellite-stable (MSS), facilitating group-specific comparisons. In the MANTIS analysis, the MSI scores of the quadruple-pathway mutant organoids mirrored those of the MSI-H group, significantly diverging from the MSS counterparts (Fig. 4a). This observation clearly indicated that MLH1KO mutant organoids recapitulated the phenotype of MMR-deficient colorectal tumors. Furthermore, to evaluate the fidelity of the MSI-H tumor landscape in the MMR-deficient organoids, we checked whether the frequent microsatellite loci48 of MSI-H cancers were targeted in our organoids. As expected, culturing of the mutant clones over time led to increased numbers of somatic mutations in the microsatellite loci found in MSI cancers (Fig. 4b). We then calculated the expected probability for acquiring specific mutations in in vitro cultured MLH1KO organoids, taking into account the specific mutational signatures and the kinetics of mutation accumulation. First, by analyzing all mutational signatures in collected mutant clones, we visualized the ‘average’ mutational signature as the mutational frequency in 96 channels (Fig. 4c). The mutational signature of the mutant clones was compared to the mutational signatures of in vitro cultured organoids and in vivo tissue derived from small intestine and colon, as previously reported (Fig. 4d)49. Next, candidate mutations known as driver mutations in the Ras, Raf, MAPK, PI3K and ErbB pathways were listed from the previous report in MSI-H cancers50. The relative probability of occurrence in each specific mutation was calculated from the average mutational signature derived from the previous approach51 and compared to the mutational rate in normal intestinal tissue and organoids (Fig. 5). As expected, KRASG12D and KRASG12V, known as the most frequent mutations in the Ras pathway, were unlikely to occur in an MMR-deficient background. BRAF mutations, which represent the most common mutations of the Ras and Raf pathway in sporadic MSI-H, were also calculated as being of low probability, which is consistent with a previous report that BRAFV600E mutation frequency is low in MMR mutation carriers3,5,52. Selection in −EGF+Afa medium resulted in NRASQ61K mutation despite being previously reported as of lower probability albeit detected in MMR-deficient CRCs53. In contrast, the high relative probability of AXIN1Q351fs corresponded to our selection result (described above). The TP53R158L mutation found in our clone showed a relatively low probability, corresponding to the occurrence of a single surviving clone in the Nut3 selection step (Fig. 2a,b). The report of TP53R158L mutations specifically in carcinogenesis related to inflammatory bowel disease suggests that the TP53 mutation is influenced by various backgrounds and mutational signatures32. Despite the low probability of TP53 mutations in our MLH1KO mutant organoids, we were able to select a TP53 missense mutation. Together, these data suggest that driver mutations inferred from the mutational signature are acquired with high probability while other oncogenic mutations with relatively low probability can also be acquired (albeit at lower frequency), shaping the genetic diversity of MSI-H CRCs.

Fig. 4: Signatures of MSI in mutant organoids.

a, MANTIS scores were calculated in the mutant organoid clones (n = 8 samples) and MSI-H (n = 38 samples) and MSS (n = 38 samples) of TCGA-COAD. We used the stepwise difference threshold set by MANTIS, which defines a score above 0.4 as MSI-H (unstable) and below this threshold as MSS (low) using TCGA-COAD. MSMuteTect2 scores were also calculated in the mutant organoid clones (n = 8 samples) and all samples used in TCGA as MSMuTect2 (n = 262 samples). The minimum values are the smallest number per MANTIS and MSMuTect2 scores. Sample numbers are written below each group. The first quartile above the whiskers represents the data point that separates the lowest 25% of the data from the rest. The center line per box plot represents the median value among the data points. The third quartile just on top of the box plot separates the lowest 75% of the data points from the highest 25%. The maximum value represents the largest score of each tool. b, MSI mutations (pink, mutated) harbored in each mutant organoid clone (columns). The top 50 most differentially unstable microsatellite loci (rows) were selected from the previous study