Tumor-intrinsic expression of the autophagy gene Atg16l1 suppresses anti-tumor immunity in colorectal cancer

Abstract

Microsatellite-stable colorectal cancer (MSS-CRC) is highly refractory to immunotherapy. Understanding tumor-intrinsic determinants of immunotherapy resistance is critical to improve MSS-CRC patient outcomes. Here, we demonstrate that high tumor expression of the core autophagy gene ATG16L1 is associated with poor clinical response to anti-PD-L1 therapy in KRAS-mutant tumors from IMblaze370 (NCT02788279), a large phase III clinical trial of atezolizumab (anti-PD-L1) in advanced metastatic MSS-CRC. Deletion of Atg16l1 in engineered murine colon cancer organoids inhibits tumor growth in primary (colon) and metastatic (liver and lung) niches in syngeneic female hosts, primarily due to increased sensitivity to IFN-γ-mediated immune pressure. ATG16L1 deficiency enhances programmed cell death of colon cancer organoids induced by IFN-γ and TNF, thus increasing their sensitivity to host immunity. In parallel, ATG16L1 deficiency reduces tumor stem-like populations in vivo independently of adaptive immune pressure. This work reveals autophagy as a clinically relevant mechanism of immune evasion and tumor fitness in MSS-CRC and provides a rationale for autophagy inhibition to boost immunotherapy responses in the clinic.

Fig. 1. Elevated ATG16L1 expression associates with poor outcome in immunotherapy of non-MSI-high CRC harboring oncogenic KRAS mutations.

a Kaplan–Meier curves indicating the association between ATG16L1 transcript levels and overall survival in atezolizumab and regorafenib treated patients (IMblaze370). Median cutoff was used to determine high and low levels separately within KRAS mutant and KRAS wildtype tumors. P-values were obtained from log-rank tests. Log-rank hazard ratios (HR) are provided with 95% confidence intervals in parentheses. b Scatter plots for KRAS mutant tumors showing the correlation between ATG16L1 transcript levels and tumor microenvironment signatures for immune, stroma and epithelial cells. Pearson correlation coefficients and two-sided t-test P-values are shown (n = 181 samples). c Immunohistochemical (IHC) staining for ATG16L1 protein in tumor biopsies obtained from IMblaze370. Dotted lines indicate tumor epithelial margins, asterisk depicts stromal component. 22 tumor biopsies were analyzed; representative micrographs are shown. Scale bar = 100μm. d, Expression of depicted lineage markers (x-axis) in major cellular compartments (y-axis) of CRC tumor tissue, analyzed by single-cell RNA sequencing (GSE146771). e Comparison of ATG16L1 transcript levels in each cellular compartment analyzed in (d). f Immune cell subsets enriched in ATG16L1-low tumors within IMblaze370, as determined by CIBERSORT gene signatures. Arrows indicate T and NK cell subsets. Unadjusted P-values from two-sided t-tests are shown (KRAS mutant n = 181, KRAS wildtype n = 113 samples). Dashed lines denote significance threshold at P < 0.05. All analysis restricted to non-MSI-high tumors. *P  <  0.05, **P  <  0.01, and ***P  <  0.001. IMblaze370 RNAseq data have been deposited to the EGA under accession number EGAS00001005952, and GSE146771 is publicly available from GEO (see Data Availability section of Methods). Source data (including exact P-values) for panel (f) are provided as a Source Data file.

Fig. 2. ATG16L1 drives CRC growth in the liver by promoting resistance to cellular immunity.

a Schematic of CRISPR engineering and generation of CRC organoids. b Immunoblot analysis of the indicated proteins in WT or Atg16l1 KO CRC organoids Three independent clones of each genotype are shown. c Hydrodynamic tail vein (HTV) injection of CRC organoids for liver growth model. d–h Bioluminescence imaging (BLI) signal quantification from CRC organoids implanted in the livers of immunodeficient hosts. d C57BL/6 J (BL6) immunocompetent hosts. ***P = 0.0001, ****P < 0.0001. e NOD-SCID/Gamma (NSG) immunodeficient hosts; WT control BL6 hosts shown in panel (d). **P = 0.0027, ***P = 0.0001, ****P < 0.0001. f CD8⁺ T cell or NK cell depletion in BL6 immunocompetent hosts. Week 5 BLI data are shown. Within tumor genotypes, treatment groups are compared using Kruskal-Wallis test with Dunn’s multiple comparisons tests. Atg16l1 KO groups are also compared to their corresponding WT groups using two-sided Mann-Whitney tests; **P = 0.0021 and ****P < 0.0001. For all Atg16l1 KO groups, n = 10 per condition. For WT groups, n = 9 for isotype control-treated mice, and n = 10 each for anti-CD8 and anti-NK1.1 treated mice. g IFN-γ KO hosts; WT control BL6 hosts shown in panel (d). ns, not significant. h Direct comparison across studies of BLI signal (from panels d–g) re-plotted as values normalized to the medians of WT organoid groups (after 4 weeks of growth). For WT/untreated hosts, n = 14 (WT organoids) and n = 11 (KO organoids), ****P < 0.0001; for CD8 T cell depleted hosts, n = 10 per group, ****P < 0.0001; for NK-depleted hosts, n = 10 per group, **P = 0.0015; for IFN-γ KO hosts, n = 11 per group, *P = 0.0233; for NSG hosts, n = 14 (WT organoids) and n = 13 (KO organoids), **P = 0.0027. WT and Atg16l1 KO organoids were compared in each condition using two-sided Mann–Whitney tests. Atg16l1 KO groups from each immunodeficient condition were also compared to Atg16l1 KO tumor growth in WT/untreated hosts using Mann-Whitney tests; ^†P = 0.0048, ^††P = 0.0032, ^†††P < 0.0001. i, j Representative BLI images of mice (top) and tumor burden in livers (bottom) from (i) BL6 (n = 14 WT; n = 11 Atg16l1 KO) hosts or (j) NSG (n = 14 WT; n = 13 Atg16l1 KO) hosts. In panels (d–h), lower and upper hinges in box plots correspond to first and third quartiles, while whiskers extend to minima and maxima. Individual data points indicate separate mice (biological replicates). In panels (d, e, and g), groups were compared using two-sided Mann–Whitney tests, with P-values adjusted for multiple comparisons using the Holm-Sidak method. All data are representative of 2-3 independent experiments. Source data are provided as a Source Data file.

Fig. 3. Tumor-intrinsic loss of Atg16l1 alters the phenotype and composition of CRC cells in the liver.

a Visualization of eight scRNA-seq clusters in UMAP dimensions: Prolif (n = 5154), Sensory/secretory (SS) (n = 1956), Mature enterocyte (ME) (n = 2378), Neuroendocrine (NE) (n = 1046), Enterocyte progenitor (EP) (n = 2467), IFN resp (n = 369), Stem (n = 1320), and Goblet/Paneth (GP) (n = 573) clusters. b, Heatmap showing the average expression of top 10 markers for each cluster in panel (a). c Density plot for the WT (left) and Atg16l1 KO (right) conditions. Cells were harvested from n = 2 mice for each condition. Low and high density are denoted by blue and red, respectively. d Proportion of cells in each cluster among all cells in the respective sample (Atg16l1 KO or WT). e Barplot indicating the proportion change of each cluster in the Atg16l KO sample with respect to the WT condition. For each cluster, the y-axis shows the log-transformed value for (proportion in KO)/(proportion in WT) ratio. Adjusted P-values are derived from two-sided Pearson’s chi-squared test for two proportions as implemented in the prop.test R function, adjusted after false discovery rate correction. f GSEA for Atg16l1 KO vs WT contrasts in each cluster using MSigDB Hallmark gene sets. Only gene sets with significant results in at least two contrasts are shown (clusterProfiler, FDR-adjusted P-values for enrichment scores derived via permutation test). g GSEA of Atg16l1 KO vs WT organoids in vitro (bulk RNA-seq) for the top 100 genes marking each scRNA-seq cluster. The x-axis denotes clusters while the y-axis shows normalized enrichment scores (normalized to mean enrichment of random samples of the same size). h GSEA enrichment plot for GP cluster from panel (g). For panels (g and h), normalized enrichment scores were computed with fgsea using voom+limma derived fold changes; P-values for enrichment scores were derived via permutation test; unadjusted P-values shown. *P  <  0.05, **P  <  0.01, and ***P  <  0.001. (IFN interferon). Source data and exact P-values for panels (e and g) are provided as a Source Data file.

Fig. 4. Tumor-intrinsic loss of Atg16l1 reprograms myeloid cell composition and phenotypes in the tumor microenvironment.

a Visualization of ten scRNA-seq clusters in UMAP dimensions: TREM2 MF (n = 7116), VSIG4 MF (n = 890), MARCO MF (n = 434), IFN responsive MONO/MF (n = 560), MONO (n = 2767), cDC1 (n = 969), cDC2 (n = 2297), CCR7⁺ DC (n = 331), pDC (n = 1115), NEUT (n = 4770) clusters. b Heatmap showing the average expression of top 10 markers for each cluster in panel (a). c Density plot for the WT (left) and Atg16l1 KO (right) groups. Cells were harvested from n = 2 mice for each condition. Low and high density are denoted by blue and red, respectively. d Proportion of cells in each cluster among all cells in the respective sample (Atg16l1 KO or WT). e Barplot indicating the proportion change of each cluster in the Atg16l1 KO sample compared to the WT group. For each cluster, the y-axis shows the log-transformed value for (proportion in KO)/ (proportion in WT) ratio. P-values were calculated from two-sided Pearson’s chi-squared test for two proportions as implemented in the prop.test R function, and adjusted by false discovery rate correction. f GSEA for Atg16l1 KO vs WT contrasts in each cluster using MSigDB Hallmark gene sets. Only gene sets with significant results in at least four contrasts are shown (clusterProfiler, P-values for enrichment scores derived via permutation test, FDR-adjusted P-values < 0.05). *P  <  0.05, **P  <  0.01, and ***P  <  0.001. (MF macrophages, MONO monocytes, DC dendritic cells, pDC plasmacytoid dendritic cells, NEUT neutrophils). Source data and exact P-values for panel e are provided as a Source Data file.

Fig. 5. Atg16l1 loss sensitizes CRC organoids to TNF + IFNγ-induced cell death.

a GSEA using MSigDB IFNγ response hallmark gene set for differential expression between Atg16l1 KO (n = 2) and WT (n = 2) CRC organoids treated (in vitro bulk RNA-seq data) with p-value = 4.78 × 10⁻¹⁵ (normalized enrichment score was computed with fgsea using voom+limma derived fold changes; P-value for enrichment score was derived via permutation test). b Differential expression of representative IFNγ-induced genes between Atg16l1 KO (n = 2) and WT (n = 2) CRC organoids in either untreated or IFNγ-treated conditions. c WT or Atg16l1 KO CRC organoids treated as indicated and stained with propidium iodide (PI). Scale bar = 1000 μm. Images are representative of 4–6 technical replicates per condition (see Source Data for Fig. 5d). d, f Cell death assayed by live imaging of WT versus Atg16l1 KO organoids (d), and Atg16l1 KO versus Atg16l1 + RIPK3 double KO organoids (f) treated with combinations of TNF + IFNγ + zVAD for 48 h. PI staining is measured by fluorescence intensity/μm². Groups compared using two-way ANOVA. e Immunoblot analysis of the indicated phosphorylated and total proteins in WT or Atg16l1 KO CRC organoids stimulated with combinations of TNF + IFNγ for 4 or 18 h. * indicates non-specific bands. Data for panels c–f are representative of three independent experiments. Summary data are shown as mean ± s.e.m. Source data for panels (d–f) are provided as a Source Data file.