Immunosuppressive niche engineering at the onset of human colorectal cancer

The evolutionary dynamics of tumor initiation remain undetermined, and the interplay between neoplastic cells and the immune system is hypothesized to be critical in transformation. Colorectal cancer (CRC) presents a unique opportunity to study the transition to malignancy as pre-cancers (adenomas) and early-stage cancers are frequently resected. Here, we examine tumor-immune eco-evolutionary dynamics from pre-cancer to carcinoma using a computational model, ecological analysis of digital pathology data, and neoantigen prediction in 62 patient samples. Modeling predicted recruitment of immunosuppressive cells would be the most common driver of transformation. As predicted, ecological analysis reveals that progressed adenomas co-localized with immunosuppressive cells and cytokines, while benign adenomas co-localized with a mixed immune response. Carcinomas converge to a common immune “cold” ecology, relaxing selection against immunogenicity and high neoantigen burdens, with little evidence for PD-L1 overexpression driving tumor initiation. These findings suggest re-engineering the immunosuppressive niche may prove an effective immunotherapy in CRC.

there were between and 50,000 -200,000 unique phenotypes depending on immune suppression and blockade values. For example, if the immune suppression parameter was 0, then there is only one immunosuppression phenotype, but if it is non-zero then there are 2 immunosuppression phenotypes (0 or ). We assumed that all phenotypes have the same intrinsic death rate, = 0.01.
The number of driver mutations, , determined which "species" the phenotype belonged to, which in turn determined the phenotype's division rate ( ), carrying capacity ( ), and intra-species competition (found in the competition matrix, ). If < 2, phenotype belonged to the epithelial species, which had a division rate of = 0.2 (Carulli, Samuelson, & Schnell, 2014) and carrying capacity of = 10 7 ; if 2 ≤ < 4, the phenotype was of an adenoma, phenotype had a division rate of = 0.2106 and carrying capacity of = 10 8 ; when ≥ 4 phenotype belonged to the carcinoma species and had division rate = 0.221 and carrying capacity = 10 9 . Increases in division rates were based in differences in observed average KI67 quadrat density between CRA and CRC, which was 1.05. This ratio was applied to the epithelial division rate, and then applied again to get the CRC division rate. The intra-species competition coefficients in are based on the assumption that epithelial and adenomas occupy different areas, while carcinomas are invasive.
= 1 0 1 0 1 1 1 1 1 Each time step, a phenotype could create mutants, the number of which was determined by drawing from the binomial distribution, where the number of trials was the population size and the probability of success the cellular mutation rate ( ), i.e. ( = , = ). The cellular mutation rate of was calculated as follows. Exome size was estimated based on information in (M. Lek et al., 2016), which describes the exome as having 45Mb, with 18Mb of possible synonymous variants. Therefore, we modeled the 45Mb-18Mb=27Mb base pairs that may produce non-synonymous variants. The "normal" per division base pair mutation rates were set to 2.91(10 −9 ) (Monkol Lek et al., 2016;Werner et al., 2020). Given these values, we defined the cell mutation rate to be the probability of at least one non-synonymous mutation in the exome per cell per division: 0.0756 = 1 − ( = 0, = 2.7 7, = 2.91(10 −9 )).
We assume all mutations are inherited, and so all mutants have a higher antigenicity that their parent. In addition, mutants may also acquire one of three types of beneficial mutations: driver mutations, the ability to protect from T-cell attack, and the ability to recruit immunosuppressive cells. The multinomial distribution was used to determine how many of each mutant to make. Given that the number of genes in the genome has been estimated to be = 20412 (Cunningham et al., 2019), the probability of creating a cell that mutated an immunosuppressive gene is 1 . Likewise, the probability of creating a cell that mutated an immune blockade gene is also 1 .If there are 25 possible driver genes, then the probability a mutant has an additional driver mutation is 25− . Finally, the probability that the new cell didn't mutate any of these beneficial mutations, and only acquired a neoantigen, was −27+ .
Therefore ⃗⃗ = ( = , ) is a four element vector where: ⃗⃗ 1 is the number of clones to create that do not have a beneficial mutation; ⃗⃗ 2 is the number of clones to create that acquire a driver mutation, with being the number of driver mutations has accumulated; ⃗⃗ 3 is the number of clones to create that express checkpoint inhibitors; and ⃗⃗ 4 is the number of clones to create that recruit immunosuppressive cells. The parameter is the number of genes in the genome, and driver is the number of possible driver mutations.
The vector ⃗ contains the probability for each mutation type, with Note that calculation of ⃗⃗ 2 assumes the number of driver mutations that will have an effect decreases as driver mutations are accumulated. In other words, it is assumed that if a driver is "hit" more than once, the benefit comes only with the first mutation, and subsequent hits have no effect. It is also assumed there is one gene for protecting from T-cell attack, and one gene for recruiting immunosuppressive cells. . Changes in the ecology were characterized using 17 cell markers and RNA transcripts associated with cytokine expression measured on whole slide images. Spatial analysis was conducted by registering the thin serial slices to create a spatially aligned composite image. The composite image was subsequently divided into 250 m x 250 m quadrats, and abundance estimated as being proportional to the number of pixels positive for each marker at 40x magnification. Direct spatial associations were then quantified using the quadrat counts from the registered whole slide images, allowing us to determine colocalization of cell types and/or cytokine expression. aITH was quantified using multi-region whole exome sequencing followed by TCR recognition potential predictions. The results from this analysis allowed us to test the model's predictions about the timing and type of immune escape in colorectal cancer. The tumor immune ecology was described using 17 markers to identify tumor cells (CK), cytotoxic T-cells (CD8), macrophages (CD68), M2 macrophages (CD163), M1 macrophages (iNOS), neutrophils (elastase), B-cells (CD20), PD-L1, DNA damage (γH2AX), proliferation (Ki67), inflammatory cytokines (TNF-α, IL-6), and immunosuppressive cytokines (TGF-β, IL-10). These 10x images come from a CRC, but the analysis was conducted at 40x. CK, CD8, elastase, SMA, CD68, Ki67, CD20, γH2AX, CD31, COX2, and PD-L1 were stained for in n=12 colorectal adenomas (CRA), n=26, "carcinoma-in-adenoma" (CIA), and n=15 colorectal carcinomas (CRC). IL-6, IL-10, TNF-α, TGF-β, CD163, and iNOS were stained for in n=9 CRA, n=9 CIA, and n=9 CRC.

Figure 4
Quadrat counts for all CRA samples. Dark blues indicates low abundances, while yellow reflects higher abundances.

Stain segmentation
The stains in each quadrat were separated using the following methods. Description of the tumor ecology within and across various stages was accomplished using a variety of methods. Cell abundances (assumed to be proportional to the number of pixels positive for each stain) were measured from the registered images, using the averaged quadrat counts (assumed to be proportional to the number of positive pixels) to determine if there were directional changes over time, i.e. if there were significant increases or decreases in abundance during progression from late adenoma to mature carcinoma. The significance of these trends was determined using a combination of frequentist and permutation statistical tests (Supplemental Figures 7 and 9). A similar analysis was conducted to determine if there were directional changes in spatial associations between the various cell types. Quadrat counts were used to construct a species association network (Popovic et al., 2019), whose coefficients were compared across stages. In addition to testing for these inter-stage trends, paired tests were used to determine if there were significant intra-tumor changes in cell abundance and interactions. This was accomplished by comparing the carcinoma (C-CIA) to its precursor adenoma (A-CIA) in the same CIA sample.
Several ecological tests were used to compare tumor-immune ecologies across stages. We quantified and compared the amount of ecological homogeneity within each tumor stage, using multivariate homogeneity of group dispersions (PERMDISP2)(M. J. Anderson, 2006). Permutational multivariate analysis of variance (PERMANOVA) was used to determine if there were significant differences in the structure of the ecological communities of each tumor stage, which was then visualized using constrained analysis of principal coordinates ( A distance matrix, which compares tumors based on the whole collection of markers, is required for many of the ecological tests we performed. We tested all combinations of distance metrics and normalization methods in the vegan R package (version 2.5-7) (Oksanen et al., 2018) and found that the Jaccard distance on log-scaled PPC most often had adjacent adenomas in CIA samples as being the most similar to one another, as would be expected. It should be noted that the Jaccard distance in the vegan package is not the same as a traditional Jaccard distance, but more of a variant of the Bray-Curtis distance. Thus, we performed the analysis on log-scaled data and a distance matrix constructed using the Jaccard distance. PERMANOVA, PERMDSIP2, CAP, and the Mantel tests were all conducted on this distance matrix.   Table showing significant changes in direct pairwise cell-cell/cell-cytokine spatial associations, where each row is a sample, and each column a spatial association. Shapes below each marker pair signify that there were significant changes in spatial association across the specified group. Green indicates that the association increased from the first group to the last group, while purple indicates that the spatial association decreased. Here, surprisingly, we observed that, compared to A-CIA, CRA have a higher association of CD68/CD163 (purple diamond), indicative of M2 macrophages, and a greater association between TGF-β & IL-10, suggesting CRA are more immunosuppressive than A-CIA. Interestingly, the same trend of increased CD68/CD163 and IL-10/ TGFβ are observed during the evolution of the tumor from progressed adenoma to carcinoma. A notable difference is that CD8 more strongly spatially co-localized with tumor (CK) in CRA compared to A-CIA and CRC (purple diamond and square), and the association between CD8 and tumor decreases from progressed adenoma to CRC (purple plus). This suggests that while CRA may have some immunosuppressive cells and cytokines, they remain most immunogenic. Figure S10. Significance tests for changes in spatial associations. A combination of frequentist statistics (y-axis) were used to determine if there were significant directional changes, where "Decrease" means there was a decrease in the abundance of a marker (x-axis) in each type of comparison. For example, all tests found a significant decrease in the spatial association between CD8 and CK when comparing A-CIA to C-CIA to CRC. Paired tests were used to determine how abundances changed in the same tumor, from the A-CIA to its neighboring C-CIA. Colored tiles show significance, whereas missing tiles indicate non-significance. KWrst= Kruskal-Wallis rank sum test, AKWTp=Approximative Kruskal-Wallis Test (permutation), GITp= General Independence Test (permutation), ALbLATp= Approximative Linear-by-Linear Association Test (permutation), JnTt= Jonckheere-Terpstra test.
A suite of statistical tests was used to detect significant changes in marker positivity, co-localization, and direct cell-cell interactions (Supplemental Figure 7 and 9). Several methods were used to find those results that were robust, as those that were found to be significant across multiple tests are most likely to be true. The two-sided Kruskal-Wallis rank sum test, two-sided approximate Kruskal-Wallis Test, and two-sided General Independence Test were used to determine there were differences between the tumor subtypes CRA, A-CIA, C-CIA, and CRC. Assuming that the carcinomas will take over the adenomas in the ca-in-ad samples, and thus develop into carcinomas, one can test for trends in the data by setting the timing order of A-CIA-> C-CIA -> CRC. Having ordered the subtypes, one-sided Approximative Linear-by-Linear Association Test and the Jonckheere-Terpstra test were used to determine if there a significant increase or decrease across groups (Seshan, 2018). The Conover-Iman test of multiple comparisons using rank sums was used as the post-hoc test to determine which subtypes were significantly different from one another (Dinno, 2017).
The paired-tests looking for differences between A-CIA and C-CIA included the General Symmetry test and the Wilcoxon rank sum test. As with the tests comparing all subtypes, the two-sided version of these tests was used to detect differences, and the one-sided test was used to detect the presence of trends and the direction of those trends.

Figure 11
Distribution of antigenicities (A), and percent of the tumor with immunosuppressive cell infiltrate (B), and percent of the tumor with blockade (PD-L1) (C), in benign adenomas compared to those of progressed adenomas and carcinomas using the subset of parameters that matched the data. The horizontal bin in each violin is colored by the average value of phenotypes that are within that bin. For each parameter combination that reproduced the observed patterns, the two most different malignant and benign tumors were selected for plotting.