KRAS mutant rectal cancer cells interact with surrounding fibroblasts to deplete the extracellular matrix

Somatic mutations in the KRAS oncogene are associated with poor outcomes in locally advanced rectal cancer but the underlying biologic mechanisms are not fully understood. We profiled mRNA in 76 locally advanced rectal adenocarcinomas from patients that were enrolled in a prospective clinical trial and investigated differences in gene expression between KRAS mutant (KRAS‐mt) and KRAS‐wild‐type (KRAS‐wt) patients. We found that KRAS‐mt tumors display lower expression of genes related to the tumor stroma and remodeling of the extracellular matrix. We validated our findings using samples from The Cancer Genome Atlas (TCGA) and also by performing immunohistochemistry (IHC) and immunofluorescence (IF) in orthogonal cohorts. Using in vitro and in vivo models, we show that oncogenic KRAS signaling within the epithelial cancer cells modulates the activity of the surrounding fibroblasts in the tumor microenvironment.

Somatic mutations in the KRAS oncogene are associated with poor outcomes in locally advanced rectal cancer but the underlying biologic mechanisms are not fully understood. We profiled mRNA in 76 locally advanced rectal adenocarcinomas from patients that were enrolled in a prospective clinical trial and investigated differences in gene expression between KRAS mutant (KRAS-mt) and KRAS-wild-type (KRAS-wt) patients. We found that KRAS-mt tumors display lower expression of genes related to the tumor stroma and remodeling of the extracellular matrix. We validated our findings using samples from The Cancer Genome Atlas (TCGA) and also by performing immunohistochemistry (IHC) and immunofluorescence (IF) in orthogonal cohorts. Using in vitro and in vivo models, we show that oncogenic KRAS signaling within the epithelial cancer cells modulates the activity of the surrounding fibroblasts in the tumor microenvironment.

Introduction
Colorectal cancer (CRC) is the third most common cancer diagnosed in the United States, with almost 150 000 estimated cases and more than 50 000 expected deaths in 2021 [1]. KRAS mutations are detected iñ 40% of all colorectal cancers [2][3][4] and they are associated with worse prognosis and more aggressive tumor behavior [5,6]. While CRC has been historically considered a uniform disease, recent evidence increasingly suggests that anatomical location is associated with important differences in epidemiological, molecular, and clinical features [4,7]. Within the setting of metastatic CRC, KRAS mutations have been associated with treatment resistance [8] as well as worse survival after liver metastasectomy [9]. In locally advanced rectal cancer (LARC), KRAS mutations have been identified as a biomarker of poor response to neoadjuvant chemoradiation (CRT), with KRAS and TP53 double mutants exhibiting particularly low rates of response [10][11][12]. While the resistance of KRAS-mt CRC to anti-EGFR therapy can be explained by the constitutive activation of the MAPK pathway [13], the mechanisms of resistance to nontargeted therapies (such as CRT) remain largely unexplored.
Prior transcriptomic studies have failed to identify reliable predictive features of response likely due to a multitude of factors, including small sample sizes, heterogeneous patient cohorts, and nonuniform definitions of tumor response [14,15]. In this study, we present a molecular analysis of LARC tumors from patients accrued through the Timing of Rectal Cancer Response to Chemoradiation (TIMING) clinical trial [16]. This is a phase 2, nonrandomized, multicenter clinical trial that evaluated the proportion of LARC patients achieving a pathological complete response (pCR) after being treated with neoadjuvant CRT followed by various doses of consolidative chemotherapy. This unique cohort provides an excellent opportunity to analyze a clinically homogenous population of patients that were selected and treated within the context of a clinical trial and that were evaluated based on uniform assessments of response. Since KRAS mutational status has been previously reported as a strong biomarker of response in this cohort [10], we sequenced mRNA for a subset of 76 patients and investigated differences in the transcriptomic profile of KRAS-mt and KRAS-wt tumors in order to gain insights into the mechanisms by which KRAS mutations increase tumor aggressiveness.
We identified dysregulation of the extracellular matrix (ECM) as the strongest discriminating feature between KRAS-mt and KRAS-wt LARC. The tumor stroma is composed of fibroblasts and the ECM, a three-dimensional arrangement of molecules that provide biomechanical and biochemical support to the surrounding cells [17]. Accumulating evidence suggests that the tumor stroma plays a crucial role in cancer progression and treatment resistance [18][19][20]. Using genetically engineered murine models (GEMM) and tumor-fibroblast co-cultures, we show that KRAS-mt epithelial cancer cells downregulate the expression of ECM-related genes in the neighboring fibroblasts.

Sample inclusion criteria
Our main cohort (LARC-TIMING cohort) consisted of 76 LARC pretreatment biopsy samples from patients accrued through the TIMING trial [16] (Clini-calTrials.gov, number NCT00335816). Local staging was performed by endorectal ultrasound or magnetic resonance imaging and patients were screened for metastatic disease with computed tomography scans. These tumors were profiled for transcriptomic analyses using gene expression microarrays, and the relevant clinicopathological information is provided as Table S1. All these patients had American Joint Committee on Cancer clinical stage II (T3-4, N0) or III (any T, N1-2) rectal adenocarcinoma with a distal tumor border within 12 cm of the anal verge. Inclusion was contingent on an adequate amount of tissue to allow for molecular analysis. Tumors with microsatellite instability (MSI) were excluded from our analysis because they comprise less than 5% of rectal cancer and present very different biologic and clinical features [2,21,22]. The study methodologies conformed to the standards set by the Declaration of Helsinki. Consent to analyze specimens from patients was obtained after approval by a central Institutional Review Board (IRB) at Memorial Sloan Kettering (MSK) and the individual IRBs at the participating institutions. Gene expression data is uploaded in GEO: (GSE170999, https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc= GSE170999).
We demonstrate the reproducibility of our findings at the transcriptional level using an external validation cohort of 71 microsatellite stable (MSS) LARC tumors from The Cancer Genome Atlas (LARC-TCGA cohort), for which RNA-Seq data are publicly available [23]. The final set of TCGA barcodes used in our analyses, together with relevant clinical features, is provided as Table S2. To further validate some of our results, we performed IHC in a set of 23 LARC pretreatment biopsy specimens (LARC-IHC cohort) from patients treated at MSK. RNA-seq was performed for validation in 15 of these 23 patients with usable biological material at the time. Due to the limited availability of pretreatment rectal cancer tissue, technical validation of some of our findings was also done on stage I-III colon adenocarcinoma tissues from patients that underwent upfront resection at MSK. This included a set of 35 colon adenocarcinoma resection specimens that were used for immunofluorescence staining (CRC-IF cohort). We also used real-time polymerase chain reaction (RT-PCR) in a set of 45 colon adenocarcinomas to validate the expression of select genes (CRC-RT-PCR cohort). Consent to analyze specimens from patients treated at MSK was obtained after approval by our IRB.

DNA and RNA extraction from tissue blocks
Pretreatment biopsies were obtained from patients in the LARC-TIMING cohort and patients in the LARC-IHC cohort at the time of diagnosis. All the specimens were stored as formalin-fixed paraffin-embedded (FFPE) blocks and ten slides of 10-micron thickness were cut per block. A pathologist reviewed the hematoxylin and eosin (H&E) slides of all the sections to confirm the boundaries of the malignant epithelia. The tumor-enriched areas were macrodissected from unstained slides guided by the H&E slides. DNA & RNA were extracted from FFPE sections using AllPrep DNA/RNA FFPE kit (Qiagen Inc., Valencia, CA, USA).

RNA profiling
Two separate platforms (Affymetrix microarray and Hiseq) were used to analyze the LARC-TIMING cohort and the LARC-IHC cohort, respectively. For the Affymetrix platform, total RNA was amplified to generate cDNA libraries using the Ovation FFPE WTA System (NuGEN Technologies, San Carlos, CA, USA) and sent for Affymetrix U133 Plus 2.0 Array (Affymetrix, Santa Clara, CA, USA). For the Hiseq platform, total RNA was used to generate libraries for mRNA deep sequencing using an adapted version of the Illumina v1.5 protocol and optimizing for reaction volume. Sequencing was done using the Illumina HiSeq 2000 Platform (Illumina, San Diego, CA, USA).

DNA sequencing and mutation calling
Two platforms were used for DNA sequencing. For the LARC-TIMING cohort, KRAS mutations were determined by standard PCR followed by Sanger sequencing of exons 2 and 3 of KRAS. MSI status was tested by PCR for the five-marker panel (BAT-25, BAT-26, MONO-27, NR-21, and NR-24) using Promega's MSI Analysis System v1.2 (Promega, Madison, WI, USA). LARC-IHC cohort specimens and CRC-IF cohort specimens were analyzed using Memorial Sloan Kettering-Integrated Mutation Profiling of Actionable Cancer Targets (MSK-IMPACT) assay to generate somatic mutational profiles from common oncogenes as published previously [24]. In select cases where the KRAS mutation calls were discordant between Sanger sequencing and MSK-IMPACT, we elected to follow the calls made by MSK-IMPACT. MSI status for the both LARC-IHC and CRC-IF cohorts was assessed using scores computed using MSISensor [25] with a cut-off value of 10.

Computational analyses
Computational analyses were performed in the 'R' statistical environment, version 3.6.0. A background correction using the robust multiarray average method was performed on the microarray data. Low count genes with an aggregate count below 50 across all samples were filtered out. Probes without matching Entrez Gene IDs and probes that map to multiple gene symbols were removed. Data normalization was conducted using conditional quantile normalization from the R package cqn. Differential analysis of mRNA expression was conducted using the 'DESeq2' package in R [26]. Significant differences were required to exhibit | log 2 FoldChange| > 1 and FDR < 0.05. Gene set enrichment analysis (GSEA) was conducted using median expression values across probes associated with each gene with the 'fgseaMultilevel' function from the 'fgsea' package in R, which is based on the adaptive multilevel splitting Monte Carlo approach. Prerank scores were calculated as -log 10 (P-value)* sign(log 2-FoldChange). The following packages were utilized in the analysis: BiocManager 1.  Fig. S2A, fold change was calculated after normalization to human 18s rRNA using the 2 ÀΔΔC t method; one-sided Wilcoxon rank-sum test was then performed. For IHC and IF, the median value per slide was taken and a one-sided Wilcoxon rank-sum test was performed. For the RT-PCR plots in Fig. 4, each experiment was done in biological and technical triplicate and fold change of gene expression was calculated after normalization to human GAPDH using the 2 ÀΔΔC t method; one-sided t-test was performed.

Tissue staining and protein quantification
Sequential 5-micron-thick slides were prepared for H&E and IF or IHC staining. IF staining was done for DAPI, VIM, FN, and POSTN by the MSK Molecular Cytology Core. IHC staining was done for POSTN. Table S8 contains the list of antibodies used. All slides were digitally scanned with Panoramic Flash (3DHistech, Budapest, Hungary) using 20×/ 0.8NA objective and appropriate filters. Scanned data were reviewed by a pathologist who was blinded to the molecular data. Relevant regions of tumor core with adjacent stroma were annotated using Case-Viewer (3DHistech) with H&E as a guide. Annotated regions were quantified using custom macros written in FIJI/IMAGEJ (NIH, Bethesda, MD, USA). For IHC POSTN signal, DAB and hematoxylin staining were separated using the Color Deconvolution algorithm and areas of threshold positive signal were measured and normalized to the area of tissue. For IF-stained signals, the number of total cells and positive cells was tabulated. The median value per slide was taken, and a Wilcoxon rank-sum test was performed for statistical analysis.

RNA isolation from co-cultures and RT-PCR
Cells were released using Cell Recovery Solution and pellets harvested for RNA isolation. Total RNA was isolated using TRIzol reagent (Ambion, Austin, TX, USA, 15596018) and RNeasy Mini kit (Qiagen, 74104), and cDNA was synthesized using a TaqMan Reverse Transcription Reagents (Applied Biosystems, Foster City, CA, USA, N8080234). All targets were amplified using TaqMan 2× Universal PCR master mix (Applied Biosystems, 4304437) and gene-specific Taqman primers and probe sets (Applied Biosystems) in ViiA 7 real-time PCR system (Applied Biosystems). Primers used are listed in Table S7. Fold change of gene expression was calculated after normalization to human GAPDH, using the 2 ÀΔΔC t method.

RT-PCR analysis from frozen tissue
Total RNA was extracted from frozen tissue specimens of the CRC-RT-PCR cohort using TRIzol reagent (Ambion, 15596018) and Qiagen's RNeasy Mini kit (Qiagen Inc.). cDNA was synthesized using a TaqMan Reverse Transcription Reagents (Applied Biosystems, N8080234). All targets were amplified using TaqMan 2× Universal PCR master mix (Applied Biosystems, 4304437) and gene-specific Taqman primers and probe sets (Applied Biosystems) in ViiA 7 real-time PCR system (Applied Biosystems). Primers used are listed in Table S7.

Study Design and Cohort Description
KRAS mutational status was available for 186 LARC patient samples from the TIMING trial that had been profiled using a combination of Sanger sequencing and  (Table S5), reveal significant downregulation of ECM and matrisome-related genes in KRAS-mt tumors. (F) GSEA enrichment plot using the ESTIMATE stromal signature shows significant downregulation of stroma-related genes in KRAS-mt tumors.

KRAS mutant LARC tumors exhibit a distinct transcriptomic profile
To gain further insight into the biological role of KRAS mutations in LARC, we investigated differences in gene expression between KRAS-mt and KRAS-wt tumors. A set of 35 genes were differentially expressed, including 21 upregulated genes and 14 downregulated genes in KRAS-mt specimens ( Fig. 1D and Table S3). We demonstrated the reproducibility of our results by using RT-PCR to quantify the expression of two significantly upregulated genes (HOXB6, P = 0.02; HOXB8; P = 0.01) and two significantly downregulated genes (COL1A1, P = 0.04; SPARC, P = 0.03) in an independent set of colon adenocarcinomas (CRC-RT-PCR cohort) composed of 30 KRAS-wt and 15 KRAS-mt tumors (Fig. S2A). The differential gene expression between KRAS-mt and KRAS-wt MSS LARC was also validated using RNA-Seq from an independent set of 71 LARC specimens from TCGA (LARC-TCGA cohort) [2,23]. Both the set of upregulated and downregulated genes that we had originally identified were also differentially expressed (P = 0.001 and P < 0.001, respectively) in the LARC-TCGA cohort (Fig. S2B, Table S4).

KRAS mutations in LARC are associated with changes in the stroma
Genes upregulated in KRAS-mt tumors included HOXB6 and HOXB8, two homeobox genes that have been previously reported to be dysregulated in colorectal cancer [30]. We observed increased levels of S100P, which encodes a member of the S100 calcium-binding proteins, and whose upregulation has been linked to increased metastatic potential and decreased chemosensitivity in colorectal cancer [31]. KRAS-mt tumors also exhibited higher expression levels of interleukin 33 (IL-33), a regulator of stromal cell activation that has been reported to contribute to the formation of a protumorigenic microenvironment [32]. We performed GSEA [33] using the Canonical Pathways collection from the Molecular Signatures Database (MSigDB) from the Broad Institute [34] and found that ECM remodeling and matrisome-related biological processes [35] clearly stood out as the top pathways being differentially regulated in KRAS-mt vs KRASwt tumors (Fig. 1E, Table S5). In particular, the set of genes downregulated in KRAS-mt tumors was strongly enriched in genes encoding ECM remodeling proteins such as periostin (POSTN), fibronectin I (FN1), thrombospondin 2 (THBS2), osteopontin (SPP1), osteonectin (SPARC), adlican (MXRA5), and several collagen alpha chains (COL1A1, COL1A2, COL12A1) ( Fig. 1D and Tables S3,S5). As the ECM is mainly regulated by fibroblasts, we postulated that changes in the stromal cells of the tumor microenvironment may be responsible for this pronounced dysregulation of ECM-related genes. We evaluated this hypothesis using the previously published ESTIMATE gene signature for stromal infiltration [19], and we found significantly reduced stromal gene expression in KRAS-mt compared to KRAS-wt tumors (Fig. 1F). Furthermore, our GSEA results showing significant differences in the ECM-related pathways and stromal signature were validated using the LARC-TCGA cohort (Fig.  S2C,D, Table S6).
We then employed IHC for orthogonal validation of our results by quantifying differences in the protein expression of POSTN, a key regulator of collagen cross-linking of the ECM, which serves as a ligand for alpha-V/beta-3 and alpha-V/beta-5 integrins to support adhesion and migration of epithelial cells [36]. In normal colorectal tissue, POSTN is expressed primarily in the lamina propria surrounding the crypts, but not in the epithelial cells or in the submucosa [36] (Fig. S3). We stained MSS LARC biopsies (n = 12 KRAS-wt; n = 11 KRAS-mt), which were obtained from an independent patient cohort. We found that POSTN staining intensity was significantly lower in KRAS-mt compared to KRAS-wt tumors (P = 0.03, Fig. 2A,B). We performed RNA-Seq in a subset of 15 cases (n = 9 KRAS-wt; n = 6 KRAS-mt), and we found that the IHC staining intensities were significantly correlated with their associated levels of mRNA expression (P = 0.02, Fig. 2C).
To further analyze changes in the stromal composition of KRAS-mt vs KRAS-wt tumors, we used transcriptomic signatures for a set of 14 distinct stromal cell populations that had been identified by performing single-cell RNA sequencing of 5993 cells from 23 colorectal cancer patients in a recently published study [37]. Consistent with the results observed in our LARC-TIMING cohort, we found that the gene expression of a majority of these stromal cell types was depleted in KRAS-mt tumors (Fig. S4A,B). In particular, differences were driven by reduced transcriptional activity of fibroblasts, including both myofibroblast-specific signatures and signatures for other stromal fibroblasts. By contrast, transcriptional activity for diverse subtypes of endothelial cells showed no significant differences or was actually increased in the KRAS-mt samples. Among the myofibroblast-specific signatures, the strongest decrease was observed in type 1 and 2 myofibroblasts, which were characterized by genes involved in ECM synthesis and degradation, respectively. Of note, no significant differences were observed in the levels of type 3 myofibroblasts, which were characterized by ACTA2 and RGS5 expression and therefore thought to have a pericyte origin.

Decreased expression of ECM-related genes in KRAS-mt tumors is not explained by differences in tumor cellularity
To resolve whether or not the changes in the stromal gene signature expression were simply due to differences in the fraction of stromal versus cancer cells in KRAS-mt vs KRAS-wt tumors, we performed multiplex IF staining on KRAS-wt (n = 23) and KRAS-mt (n = 12) MSS colon cancer specimens (Fig. 2D,E). Levels of vimentin (VIM), an intermediate filament protein expressed by fibroblasts, were not different in KRAS-mt and KRAS-wt tumors, demonstrating that the number of fibroblasts was comparable. This is consistent with the lack of significant differences in VIM mRNA expression between KRAS-mt and KRAS-wt tumors in our LARC-TIMING cohort (q = 0.2). By contrast, POSTN and FN1, which are two ECM-associated genes that were differentially expressed in our KRAS signature, had diminished protein expression in KRAS-mt compared to KRAS-wt specimens (Fig. 2D,  E). We also analyzed computational estimates of tumor purity (defined as the ratio of cancer vs noncancer cells in the sequenced samples) for the LARC-TCGA cohort that had been obtained using the ABSOLUTE algorithm based on SNP6 array data from TCGA [38], and no significant differences were observed between KRAS-mt and KRAS-wt specimens (Fig. S2E). These results are all consistent with the observed transcriptional signature being predominantly related to a fibroblast cell state change rather than cellularity effects.

Mutational activation of KRAS induces the observed stromal signature
To demonstrate that the introduction of a novel KRAS mutation can lead to the type of downregulation of stromal and ECM genes that we have observed in our LARC patients, we used a GEMM of CRC whereby Apc can be conditionally suppressed using a doxycycline-regulated shRNA and mutant Kras (G12D) is activated under a Cre recombinase [28] (Fig. 3A). Consistent with our observations in human tissues, we observed lower POSTN and FN1 protein expression, but no changes in VIM in KRAS-mt (n = 3) vs KRAS-wt (n = 3) mouse tumors (Fig. 3B,C). ECM- related genes that were downregulated in KRAS-mt human tumors also exhibited lower levels of expression in tumors from KRAS-mt mice (Fig. 3D).
Furthermore, tumors from KRAS-mt mice were associated with a diminished stromal signature based on GSEA using the ESTIMATE gene set (Fig. 3E).

Contribution of fibroblasts to the KRASassociated gene signature
The KRAS-associated gene signature described above could be driven by intrinsic changes in the KRAS-mutated epithelial cancer cells or it could be the result of changes in the levels of expression of surrounding stromal cells. To investigate this question, we examined mRNA expression levels of KRAS-wt and KRAS-mt Caco2 cells cultured in isolation. As the genomic background of Caco2 is KRAS-wt [29], we created a KRAS-mt line by inducing a KRAS G12V mutation on one allele. The differences in expression of matrisome and ECM-related genes that we had observed in solid tumors were not replicated in KRAS-mt vs KRAS-wt Caco2 cells (Fig. 3D). Similarly, no significant differences were observed in terms of the stromal signature (Fig. 3F). These results suggest that the KRAS-associated gene signature that we had identified in LARC reflects changes in the tumor stroma, rather than in the epithelial cancer cells.
To better understand the contribution of the stromal component to the KRAS-associated gene signature, we co-cultured mouse fibroblasts (NIH-3T3) with human Caco2 KRAS-mt or Caco2 KRAS-wt tumor cells. Since the co-culture is a mixture of mouse (fibroblast) and human (cancer cell line) cells, we used mouse-specific RT-PCR primers to quantify the expression of stromal genes from the fibroblasts. Expression of both mouse Fn1 and Postn, but not Vim, was diminished in the Caco2 KRAS-mt co-culture compared to the Caco2 KRAS-wt co-culture (Fig. 4A). The reduction in mouse Fn1 and Postn expression was not observed when the mouse fibroblasts and the Caco2 cells were co-cultured physically separated by a microporous transwell membrane. These experiments suggest that the KRAS-associated changes in gene expression that we have described above require the presence of fibroblasts and that these fibroblasts need to be in close proximity with epithelial cancer cells.
To further ensure the reproducibility of our results, we replicated this experiment in a 3-dimensional culture setting using colonic organoids isolated from genetically engineered mice. KRAS-wt (shApc/Kras wt ) and KRAS-mt (shApc/Kras G12D ) mouse organoids were cultured with the human fibroblast cell line CCD-18co in direct co-culture or indirect transwell setting (Fig. 4B). We also employed a different set of KRAS-wt and KRAS-mt mouse organoids (Apc −/ − Trp53 −/− Kras wt and Apc −/− Trp53 −/− Kras G12D ) as these co-occurring mutations are predominant in colorectal cancer and thought to be more biologically relevant.
Quantification of the stromal genes expressed from the human fibroblasts was done with human-specific RT-PCR primers. We found reduced expression of the stromal genes FN1 and POSTN in the human fibroblasts co-cultured with KRAS-mt compared to KRASwt organoids (Fig. 4B). The difference in FN1 and POSTN gene expression was not observed when the organoids and fibroblasts were co-cultured in a transwell setting. Therefore, our data are consistent with a model where KRAS-mt rectal cancer modulates the gene expression of surrounding fibroblasts in the microenvironment.

Discussion
In this study, we showed that KRAS-mt LARC tumors are associated with worse response rates across all treatment groups and CMS clusters, highlighting the profound biological effects driven by KRAS. We found that KRAS-mt tumors exhibited diminished expression of ECM-related genes compared to KRASwt tumors. We used GEMMs to show that the introduction of a KRAS mutation in the CRC cells altered the expression of ECM genes from fibroblasts in a manner that was consistent with the patterns that we had observed in human tumors. Our results using IF, as well as DNA-based estimates of tumor purity, were consistent with a change in the transcriptional program of fibroblasts, rather than a cellularity effect. We also analyzed the subtypes of stromal cells and discovered that the gene expression levels of stromal fibroblasts and myofibroblasts were diminished in KRAS-mt tumors, suggesting that oncogenic KRAS affects a diverse group of fibroblasts and not just a specific subtype.
RAS-mediated modulation of the ECM in CRC has clinical and translational relevance because remodeling of the ECM is thought to play important roles in tumor growth, invasion, and treatment resistance [39,40]. The CMS3 subgroup-which is enriched in KRAS-mt tumors-has been shown to exhibit diminished transcriptomic signatures of stromal infiltration [21]. In this respect, our finding that a KRAS mutation is associated with reduced expression of ECM proteins is consistent with the CMS classification. However, not all tumors classified as CMS3 are KRAS-mt, which suggests that other mechanisms may also contribute to the modulation of the ECM in CRC.
The effects of KRAS signaling in the tumor microenvironment of solid tumors have been extensively studied [41], but the role of KRAS in modulating the ECM remains poorly understood. The conditioned medium from a KRAS G12V transformed colon epithelial cell line has been shown to play a role in modulating the migration but not the expressive state of myofibroblasts [42]. Similarly, in our indirect co-culture experiments, we did not observe changes in the gene expression of the fibroblasts. However, our direct co-culture experiments showed that KRAS-mt epithelial cancer cells downregulated the expression of ECM genes in the surrounding fibroblasts, suggesting that physical proximity of the two cells may be necessary for this mechanism. Contrary to the work in colorectal cancer, studies using pancreatic adenocarcinoma models have found that KRAS G12D mutations activate stromal cells and upregulate multiple ECM proteins [43]. These findings suggest that the mechanisms behind KRAS modulation of the ECM may be cancer-type specific. The unique set of human rectal cancer specimens that we have analyzed provides an anatomically homogenous population in which to study both cancer cell-intrinsic and stromal changes driven by KRAS mutations. Gene expression microarrays were used for transcriptional profiling of patients in our main cohort because this was the most prevalent sequencing platform at the time when specimens from the TIMING trial were originally collected and processed; however, we have shown that our results can be reliably reproduced using alternative approaches such as RT-PCR and RNA-Seq. Genomic profiling of KRAS mutational status was done using Sanger sequencing in a large subset of our cases, and therefore, we cannot exclude the possibility that other important members of the RAS signaling pathway such as NRAS or BRAF were mutated in a small fraction of our KRAS-wt samples. While our patient cohort presented a diverse catalogue of KRAS mutations, we selected KRAS G12D and KRAS G12V to study KRAS-mt tumors in our experimental models. This might limit the generalizability of these results as different KRAS mutations have been hypothesized to have diverse tumorigenic effects [12,44]. TP53 and KRAS double mutant tumors have been reported to exhibit particularly low rates of response to CRT in LARC [12], but we have not analyzed them as a separate entity in this study due to the small sample size. Finally, a number of previous studies have shown that KRAS mutant cancer cells have immunoregulatory effects that extend beyond the ECM changes reported in our study to affect other elements of the tumor microenvironment [41]. Future analyses combining single-cell sequencing technologies and multiplex immunofluorescence will make it possible to better investigate this clinically relevant question.

Conclusions
In conclusion, we present the first transcriptomic analysis of samples accrued through the TIMING trial for patients with LARC. Our results show pronounced remodeling of the ECM in KRAS-mt tumors and altered transcriptional programs in their surrounding fibroblasts mediated by KRAS. While more work is clearly needed to better understand how this contributes to the malignant potential of KRAS-mt rectal cancer, our analyses constitute a step forward in the characterization of the dysregulated biological processes orchestrated by KRAS.
writing-original draft, visualization, and project administration; and JGA involved in conceptualization, resources, writing-review/editing, and funding acquisition.

Supporting information
Additional supporting information may be found online in the Supporting Information section at the end of the article. Fig. S1. CMS clustering of samples in the LARC-TIMING cohort. Heatmap shows varying levels of gene expression across samples stratified by consensus molecular subtype (CMS). Within each CMS subtype, samples were further grouped according to KRAS mutational status. Genes were grouped vertically into relevant sets using a previously published set of gene expression signatures of biological relevance in colorectal cancer. Fig. S2. Validation of KRAS-mt transcriptional signature using tumors from the LARC-TCGA and the CRC-RT-PCR cohorts. A: RT-PCR validation of selected differentially expressed genes in the CRC-RT-PCR cohort (n = 30 KRAS-wt, n = 15 KRAS-mt). Experiments were done in technical triplicates; three independent experiments were performed. Fold change was calculated after normalization to human 18s rRNA using the 2 ÀΔΔCt method. Error bars show standard error of the mean. a one-sided Wilcoxon ranksum test was performed. P < 0.05 is denoted by asterisk. B: External validation of upregulated and downregulated gene sets using the LARC-TCGA cohort (n = 45 KRAS-wt, n = 26 KRAS-mt). The volcano plot is used to visualize concordance between genes that were upregulated and downregulated in KRASmt vs KRAS-wt across the two cohorts. C: Validation of the top 15 pathways using the LARC-TCGA cohort. D: Validation of the stromal signature using the LARC-TCGA cohort. P-values were computed using gene set enrichment analysis (GSEA). E: Tumor purity estimates of KRAS-wt (n = 45) and KRAS-mt (n = 26) samples in the LARC-TCGA cohort, computed using the ABSOLUTE algorithm (P = 0.38). Box whisker plot is shown; whiskers represent 1.5× the interquartile range. P-value is computed by one-sided Wilcoxon rank-sum test. Fig. S3. Analysis of VIM, POSTN, and FN1 protein expression in normal colon. A: Representative images of IF staining in normal human colon. Scale bars correspond to 2mm on the left panels and 50 μm on the right panels. B: Representative images of IF staining in normal mouse colon. Scale bars correspond to 200 μm on the left panels and 50 μm on the right panels.  Table S1. Clinicopathological Features of LARC-TIMING Patients. Table S2. List of LARC-TCGA Barcodes. Table S3. Differentially Expressed Genes in KRAS-mt vs KRAS-wt Patients from the LARC-TIMING Cohort. Table S4. Differentially Expressed Genes in KRAS-mt vs KRAS-wt Patients from the LARC-TCGA cohort.