Introduction

Colorectal cancer (CRC), the third cause of death among cancer patients, is a highly heterogeneous disease, with a slow initial progression that favors the accumulation of mutations leading to a complex phenotype. Differences that exist both between and within tumors of the same cancer type are a major hurdle towards proper treatment selection and for developing more targeted therapies. Depending on the perspective under which these differences are investigated, various categorization paradigms have emerged. The systematization of clinical and histopathological parameters led to the definition of current TNM staging system (1), which presently constitutes the gold standard for diagnosis and prognosis. The development of high throughput molecular technologies brought a novel perspective and set the stage for the appearance of molecular taxonomies categorizing the tumors into subgroups sharing common molecular traits (2–7) with consensus molecular subtypes (CMS)(8) representing their common denominator. While these studies were based on whole-tumor (bulk) gene expression data, the developments in single-cell sequencing further refined the CMS classes adding two intrinsic epithelial subtypes (iCMS2/3) to the picture (9). Other studies combined genomics and transcriptomics data and an alternative classification emerged (10).

Whole transcriptome expression profiling of tissue sections is generally performed on RNA extracted from regions of interest covering diverse cell collections. By consequence, the expression levels associated with various transcripts represent, in the end, a weighted mean of contributions of each cell type, being driven by the most abundant ones. The signals from less abundant cell types are reduced or even silenced and are, therefore, overlooked. In the case of solid tumors, this approach requires a representative region, enriched in tumoral cells, to be selected in the tissue section(s) and used for RNA extraction. This is the predominant approach to tissue expression profiling that fueled the myriad of studies over the last two decades and led to significant progress in understanding the various cancers. Newer technologies such as single cell sequencing and spatial transcriptomics allow for a much finer selection of cells to be interrogated (11,12). However, while powerful, these techniques rely on fresh tissue and have still to find their place in routine clinical practice.

The importance of a morphological perspective on the molecular classification has been acknowledged from the beginning, Jass (13) already identifying several morphological features associated with the five groups proposed (e.g. serration, mucinous and poor differentiation were highly present in two of the five groups), but also noted that these features were not sufficient for predicting the groups. Later, Budinská et al. (3) proposed six morphological patterns (morphotypes) as major histological descriptors and showed that a two-tier histological score is strongly associated with the five molecular subtypes identified. Interestingly, a pure data-driven image-based classifier for the same molecular subtypes resulted in selecting remarkably similar morphological motifs (14,15). Müller et al. (16) reviewed the TCGA and CMS subtypes and their links with some morphological aspects, most notably the serrated phenotype. It is worth mentioning that in all these cases the evaluation of the morphological features referred to the whole tumor section; for example, a tumor was considered of mucinous morphology if the mucinous pattern was present in more the 50% of the tumor region, in accordance with standard definitions endorsed by the World Health Organization (17). These links between tumor morphology and molecular features also imply that the gene expression profile may depend on the tumor region sampled for RNA extraction. The sensitivity of gene-based classifiers to tumor sampling raised concerns regarding the stability of consensus molecular subtypes (18) and may partially explain the low proportion of biomarkers that reach clinical relevance (19).

It is evident that, while intra-tumoral heterogeneity is recognized as a major challenge, we still lack the practical tools for its characterization that would easily translate into a diagnostic and predictive model. In contrast with previous results, in our study we explored region- (morphotype-) based transcriptomics approach as a possible solution to this problem. This method offers a trade-off between whole-tumor profiling and spatial transcriptomics. It has a better signal resolution than whole-tumor profiling, since it selects tumor regions with more similar cellular buildup, and covers the whole transcriptome, but clearly has a much lower spatial resolution than true spatial transcriptomics. However, it represents a practical approach where several regions of interest can be stably identified, and their profiling could be easily integrated in the current molecular pathology diagnostic practice.

Building on our previous results (3), we based our study on a detailed exploration of the transcriptome of the six morphotypes identified earlier as associated to the molecular subtypes of CRC: complex tubular (CT), desmoplastic (DE), mucinous (MU), papillary (PP), serrated (SE) and solid/trabecular (TB), respectively. As reference, we also profiled several tumor-adjacent normal (NR) and supportive stroma (ST) regions. The present study was based on a single center cohort and was designed to achieve several goals: (i) identifying representative samples for each of the morphotypes, (ii) providing a comprehensive characterization of their transcriptomics landscape, and (iii) studying the intra-tumoral heterogeneity from the perspective of morphotype-resolved transcriptomics. We characterized the morphotypes from several transcriptional angles: basic molecular programs as captured by differential expression and pathway analyses, molecular tumor classifiers and prognostic gene signatures. At the same time, we looked for variations both across all tumors and across matched samples (within tumors). The emerging picture is of an unexpectedly high heterogeneity, with clear implications both from fundamental biological and practical perspectives, opening new avenues for biomarker design.

Results

Data

From n=111 unique cases of primary CRC tumors (stages: II: 59, III:32, IV:20), n=202 regions were macrodissected representing either tumor morphological regions (n=149), tumor-adjacent normal tissue (NR, n=17), supportive stroma (ST, n=8), or whole tumor (n=28), respectively (Supp. Fig. SF1), each case being represented by one or several regions. Among the tumor morphological regions, n=126 “core samples” were identified based on “morphological purity”, indicating regions containing at least 80% of a unique morphological pattern. The six morphotypes of interest (Figure 1) consisted of (in brackets the additional non-core samples) 41(+11) CT, 13(+2) DE, 18(+3) MU, 10(+2) PP, 33(+7) SE, and 9 TB samples, respectively. The distribution of associated main clinical parameters is given in Supplemental Table ST1. The only statistically significant associations found were between MU or TB and grade 3 tumors, and SE and lower grade tumors (p=0.019, Supplemental Table ST2), respectively.

Figure 1

To complement the results presented here, we created a web application https://morphogene.recetox.cz allowing the interrogation of gene expression in various morphological regions.

Morphotype cellular admixtures

The transcriptomic profile of solid tumor sample is a mixture of gene expression profiles of individual cell types and their specific programs, including cancer cells at different levels of differentiation, specific immune cells, or supportive fibroblasts. As a first step, we performed in-silico deconvolution of the expression profiles to identify the most prevalent cell types in each of the morphotypes and GSEA to score cell type-specific gene sets (see Materials and Methods) in each morphotype, and NR and ST regions (used as controls, Figure 2, Supplemental tables ST3-4).

Figure 2

The results from ESTIMATE indicated, as expected, a high stromal content for ST, DE and MU and a high epithelial tumor cell content for normal region and TB and SE morphotypes, respectively (Figure 2A). A more balanced situation was observed for CT and PP morphotypes (similar to NR). This agreed with the (stroma-related) “Isella signatures” (20) where ST, DE and MU were enriched in endothelial cells, CAFs and immune cells (Figure 2B). When investigating the categories of epithelial cells, the signatures of top of the normal colon crypt cells (21) and colon differentiated epithelial cells (22) were enriched solely in NR regions, while DE, MU, CT, SE and TB were depleted in these cell types (Figure 2B). On the other hand, MU, CT, PP and TB regions expressed genes specific for the basal crypt cells (21) and ST, DE and MU were enriched in signatures of intestinal stem cells. These observations are in perfect agreement with the definition of the morphotypes and confirm the proper selection of the samples. quanTIseq revealed that all tumor morphotypes were enriched in M1 macrophages (with maximal presence in MU and DE), while M2 macrophages, NK cells and myeloid dendritic cells where highly present in supporting stroma and tumor-adjacent normal regions (Figure 2C). Additionally, TB morphotype had the lowest scores for regulatory T cells (TREGs) and B cells.

Further we refined the morphotype cell admixtures by testing signatures of different cell types and their active programs as derived from single-cell sequencing studies. We evaluated more than 150 signatures of stromal, epithelial, and immune cell population (supplemental tables of (23)) and cancer associated fibroblasts (CAFs) (24,25) (see Supplemental table ST3 and ST4 for full signatures). Interestingly, the morphotypes differed in the signatures of CAFs subpopulations (Figure 2D). ST, MU, and DE had high GSEA scores of most of the CAFs subpopulations, while the rest (CT, PP, SE, and TB) had mostly negative scores, indicative of depletion of corresponding cell types. DE and MU were most strongly enriched in signatures of ECM-myCAF S1 – associated with immunosuppressive microenvironment and pro- metastatic functions (25) – and wound-healing myCAF S1 populations, while the adjacent stroma mainly showed signatures of normal fibroblasts, detox-iCAF S1 and IL-iCAF S1 populations, both characterized by detoxification and inflammatory signaling. NR regions were enriched only in normal fibroblasts and detox- iCAF S1. By exploring the signatures from (23), we observed even finer differences between morphotypes within all three cell type populations and their programs (Supp. Fig. SF2 A-C). For instance, CT, TB, and SE had enriched pS04 (ribosomal) and pS12 (proliferation) stromal cell signatures, in addition, CT and TB expressed pS05 (interferon-stimulated genes, ISGs) and pS21 (FOS, JUN) signatures. Also, NR had a specific enrichment in mitochondrial (pS09), metallothionein (pS16) and BMP-producing (pS17) fibroblasts. CT and TB resembled MU in expressing pS20 signature and, additionally, TB showed similar levels of pS13 (inflammatory) signature as MU and DE. ST regions and DE and MU morphotypes had significantly increased pS02 (Fibro.matrix/stem cell niche) signature. Full results for other cell types and programs are provided in the Supplemental Table ST4.

CRC morphotypes and molecular programs

The molecular programs and pathways represented in MSigDB were scored by performing GSEA on differentially expressed genes (DEGs) in all morphotypes (and NR and ST).

For the first analysis, the ordered lists of DEGs per morphotype were obtained by contrasting the individual expression profiles to the average profile of pooled samples (Supplemental Table ST5 contains all DEGs). This allowed the identification of all molecular programs significantly de-/activated in each morphotype (Table 1, Figure 3, Supplemental Table ST6). When considering only the hallmark signatures (H collection), the discriminative gradients between the morphotypes (and NR and ST) were along the EMT and TNF-α signaling axes at one end, and the MYC and E2F targets at the other end (Table 1, Figure 3A). Desmoplastic and mucinous shared active pathways involved in immune system response (TNF-α signaling via NF-κB, interferon gamma response, complement, IL2-STAT5 signaling), neoangiogenesis, and increased metastatic potential (EMT, coagulation, TGF-β, NF-kB, NOTCH, Apical junction). At the other end of the spectrum, CT and TB morphotypes had activated major pathways involved in proliferation processes (P53, MTORC 1, Myc targets, G2M checkpoint, Mitotic spindle, NOTCH signaling, Protein secretion). In contrast with CT, TB morphotype shared with MU and DE active TGF-β signaling, apoptosis, and most pathways involved in immune system response. PP and SE morphotypes had activated MYC and E2F targets, with PP morphotype exhibiting downregulation of the KRAS signaling and upregulation of the WNT-β catenin signaling.

Figure 3

Table 1.

We performed principal component analysis (PCA) of the GSEA scores of hallmark pathways. Their projection onto the first two principal components revealed a specific bi-dimensional clustering of the morphotypes and illustrated the gradient of changes between morphotypes (Figure 3B, Supplemental Figure SF3). At one end, MU and DE shared the same region in PCA space with positive coordinates on the axis defined, among others, by EMT, inflammatory response, and UV response. At the same time, they had opposite projections on the second axis of variation, defined by p53, unfolded protein response and cholesterol homeostasis. In contrast, SE and PP shared the same quadrant with negative coordinates on the first axis, but positive on the second axis. The CT and TB fell between the two previous groups with respect to the first axis of variation, while having similar activations of pathways defining the principal components. Overlaid on top of the transcriptomics layer, an additional gradient could be observed: epithelial cell differentiation. Indeed, while SE, PP and CT were well or moderately differentiated, TB, DE and MU had low or undifferentiated morphology.

Figure 3C shows heatmap of median expressions of all top 5 up- and down-regulated genes of each morphotype with respect to the average profile (full lists in Supplemental Tables ST5-6). A second analysis identified morphotype-specific processes and pathways by GSEA of differentially expressed genes between each morphotype and all other five (excluding ST and NR) (Supplemental Tables ST7-8).

Several macro dissected regions originated from the same section allowing for paired comparison of morphotypes. While the reduced number of such pairs (MU vs SE: 8 pairs, DE vs SE: 7, CT vs DE: 5, and CT vs MU: 5, respectively) impacted the statistical power, we were able to identify genes differentially expressed (after p-value adjustment) in all but MU vs SE, indicative of regional differences (Supplemental Tables ST9-10). The differences between gene expression signatures from the matched paired comparisons were in line with those from comparisons not accounting for sample pairing, indicating that the morphotype specific effect was dominating the contrasts (see Supplemental Figure SF4).

We also performed comparison between all pairs of morphotypes (Supplemental Tables ST11-12). This comparison shows that, despite similar content in terms of fibroblasts or epithelial cells (discussed above), there are still differences both in terms of differentially expressed genes (ST11) and activated molecular programs (ST12) between DE and MU, on one side, and CT, PP, SE, and TB. These results refine those presented above and allow an ordering of morphotypes in terms of relative activation of pathways. For example, KRAS signaling appears to be highest in PP, followed by CT.

Morphotypes and molecular subtypes

The molecular subtyping taxonomies of CRC were derived from datasets representing profiles of whole tumor sections, therefore aggregating the expression of many cell types. In our previous work (3), we associated molecular subtypes with morphotypes assessed on the whole tumor and hence we were interested to see how this observation translated to the case of macrodissected morphological regions. We predicted both the consensus (CMS)(8) and intrinsic (iCMS) (9) molecular subtypes.

All ST regions were predicted as CMS4, and 82.4% of NR regions as CMS3. For the morphotypes, the predictions were more distributed across subtypes: DE and MU were most often assigned to CMS4 (63.6% and 58.8%), PP, SE, and CT to CMS2 (62.5%, 41.7% and 41.9%) and TB to CMS1 (80%) (Figure 4A, Supplemental Figure SF7). More importantly, this heterogeneity was also observed intra-tumoral, with regions within the same tumor section being assigned to different subtypes (Figures 4A, 5, Supplemental Figure SF5A-B).

Figure 4

In contrast, intrinsic molecular subtypes (iCMS2/3) were much more stable, most of the time all the morphotypes within a tumor sharing the same iCMS label (Figure 4B, Supplemental Figure SF7) and agreeing with the whole-tumor assignment. NR, MU, TB, and ST regions were classified most of the time as iCMS3 (100%, 94,4%, 71.4%, 66.7%), while PP, CT and DE were predominantly classified as iCMS2 (77.8%, 70.3%, 66.7%). The serrated morphotype was almost equally assigned to each of the iCMSs (iCMS2: 58%, iCMS3:42%).

Prognostic and predictive gene-based signatures

The morphotypes generally differed in terms of score distributions, with two signatures reaching statistical significance (Kruskal-Wallis’s test: Eschrich p=0.0228, Jorissen p=0.00085, Figure 4C-D). A more pronounced variability was observed when comparing tumor regions to matched whole tumor, with amplitude of the differences (region vs whole tumor) larger than 50% of the whole tumor score in some cases (Figure 4C). Figure 5 shows a case study with three different morphological regions (CT, MU, SE) which manifest rather large deviations from the whole tumor-based risk scores for most of the prognostic signatures (see also Supplemental Figure SF5A-B).

Figure 5

The predicted resistance/sensitivity to different therapeutics varied across morphotypes: MU resistance to gefitinib; DE sensitivity to azaticidine, dasatinib, and aplidin, and resistance to tamoxifen and gefitinib; PP resistance to cantharidin, SE resistance to aplidin, CT sensitivity to alkylating agents (Supplemental Figure SF6B). The differences were observed even within tumor (Figure 5), with some of the supposedly sensitive tumors (whole tumor scoring) having regions of predicted resistance (Supplemental Figure SF6A- B).

Discussion

The analysis of the morphotypes from transcriptomics perspective is meant to bridge the histopathology and gene expression. The present exploratory study was motivated by our earlier observations linking the morphological aspects of CRC to the molecular subtypes (3). The original observations semi-quantitatively scored the morphotypes as primary or secondary dominant in the whole tumor section and showed that subtype A (corresponding to CMS3) was enriched in PP and SE morphologies, subtype B (corresponding to CMS2) in CT morphology, subtype C (corresponding to CMS1) in MU and TB morphologies, and, finally, subtype D (CMS4) in DE/stromal reaction (3). In contrast, here we focused on tumor regions rather than whole tumor, which also allowed the characterization of the intra-tumor heterogeneity.

The results show a whole landscape of changes at gene and pathway levels, with morphotypes residing on a continuum space of molecular descriptors. The analysis of hallmark pathways and selected signatures combined with in silico deconvolution of cellular admixtures served two purposes. First, to confirm that the samples exhibit known properties (e.g., TB, SE and PP have high tumor epithelial cell content, and DE and MU are enriched in fibroblasts; molecular EMT signature is high in MU and DE, but low in CT and SE, etc.), thus ensuring proper quality of the data. Second, it served to refine the characterization of the morphotypes and sketching their “molecular portraits”. The morphotypes investigated had a fluid characterization from a transcriptomics perspective, with many pairwise similarities and some striking differences. Even from a strict histopathology perspective, it was difficult, if not impossible, to clearly distinguish the separation between adjacent morphological regions therefore a certain degree of contamination between morphotypes was to be expected. Nevertheless, the enrichment in specific cell types and states allowed the identification of characteristic molecular features.

MU and DE morphotypes (previously associated with CMS1 and CMS4, (3)), as expected, exhibited high score of genes up-regulated in colon fibroblast TGF-β signaling pathway, genes associated with high tumor stromal content, CAFs and endothelial cells as well as pathways involved in immune system response. The detailed analysis of CAFs in fibroblast-rich regions (DE, MU and ST) based on signatures derived from single cell sequencing studies (23,25) revealed some finer differences: the supportive stroma (ST) region had the complete panel of fibroblast tested, while DE and MU most notably missed the “normal CAFs”. The main difference between DE and MU appeared to be that former was enriched in CAFs associated with inflammatory response (IL-iCAF), all the other CAFs being present at similar levels. The other morphotypes were either significantly depleted in fibroblast signatures or their GSEA scores were not statistically significant. Deconvolution of immune cell fractions by quantiSeq showed enrichment of DE and MU in M1 macrophages. Given the involvement of CAFs in modeling the tumor microenvironment through ECM remodeling, angiogenesis promotion and immune system regulation (26), our results support the idea of scoring separately the stromal component by either molecular or histopathology descriptors, in addition to tumor regions themselves. Even though DE and MU (and ST) also had the highest scores for the molecular EMT signature, our observations rather support the description of CMS4 as stromal/desmoplastic subtype than “true” mesenchymal, in agreement with (27). Further, the poor prognostic associated with CMS4 could be explained by the stromal component: both (28) and (29) agree that a high stromal invasion/desmoplastic reaction is prognostic of shorter time to relapse.

CT morphotype represents a classic adenocarcinoma and is one of the most common morphologies. In our previous study (3), this morphotype was associated mainly with subtype B (vastly overlapping with CMS2). TB morphotype seems to be mostly representative of higher-grade tumors and was associated with CMS1. In contrast to NR, CT and TB showed significant enrichment of signatures of normal colon basal cells. From the molecular perspective, CT together with TB had both activated major pathways involved in proliferation processes. TB, in addition, resembled MU and DE morphotypes by sharing active TGF-β signaling, apoptosis, and active immune system response. SE and PP morphologies may be indicative of a different oncologic pathway – the “serrated pathway” (30). The two morphologies share common features like well to moderately differentiated, with low stromal content and crypt structure still preserved. From a molecular perspective, we found that both SE and PP were both distributed similarly across molecular subtypes (both CMS and iCMS) and had similar activation of hallmark pathways: EMT, IL2/STAT5, IL6/STAT3, KRAS signaling all being down-regulated, while MYC targets being up-regulated. Among the hallmark pathways, androgen response, heme metabolisms and IL6/STAT3 (all silenced), appeared to be specific (and statistically significant) to SE and PP.

Given the relatively small sample size and similarities already observed between the morphotypes, it came as no surprise that the lists of differentially expressed genes, morphotype-specific, were generally short (for FDR σ; 0.15). Nevertheless, literature search of the genes on top of these lists showed importance of these genes in CRC development, progression, EMT transition or response to therapy. For CT, the top gene was PIP5K1B which was related to PI3K/AKT signaling and seems to be involved in colorectal cancer development (31). TB had the most differentially expressed genes (n=662) in comparison with all other morphotypes, with top genes including FBXO5 – prognostic of shorter time-to-relapse in various cancers (32), FLRT3 – a proapoptotic gene which, when overexpressed, inhibits EMT (33), SETSIP – gene coding chromatin-binding protein capable of participating in fibroblast reprogramming and differentiation into epithelial cells (34), E2F7 – up-regulated by p53 in response to DNA damage (35), CXCL14 (downregulated) - depending on the cell of origin can have both tumor suppressive or supporting role (36), SEMA5A (downregulated) gene – proposed as prognostic marker in CRC (37). Among top overexpressed genes specific to MU morphotype we found FGF7 (fibroblast growth factor 7) whose disrupted signaling was associated with deregulation of cell differentiation (38), and MUC2 (intestinal mucin) whose downregulation has been suggested a marker of adverse outcome (39). At the same time, MUC2 was also among the DE-specific genes, but downregulated, consistent with the observation that desmoplastic reaction is a marker of shorter relapse-free survival (40). Still in DE, we found as top overexpressed genes PIEZO2 – a paralog of PIEZO1 which is involved in colorectal cancer metastases (41), SLIT3 – a member of the Slit/Robo pathway, a major regulator of several oncogenic pathways and potential therapeutic target (42), and OLFML2B – a potential biomarker for resistance to MEK inhibitors (43). SE morphotype had only one specific gene overexpressed at FDR σ; 0.15, CCDC175. At the other end of the list, very interestingly, we found significantly downregulated gene for dihydropyrimidine dehydrogenase (DPYD) gene – the variants of which are predictive of 5-fluoruracil toxicity in adjuvant colon cancer treatment (44), GLIPR2 which participates in positive regulation of ERK1/2 cascade and EMT transition (45), or the HOXA9A gene, the overexpression of which was suggested to contribute to stem cell overpopulation responsible for development of CRC (46) or the GLI3 gene – that participates in sonic hedgehog (Shh)-Gli-mediated tumorigenesis and the loss of Gli3 signaling was shown to initiate cell growth inhibition in colon cancer cells, while sensitizing colon cancer cells to treatment with anti-cancer agents (5-FU and bevacizumab) (47). The only specific gene marker of the PP morphotype was the downregulation of MZP – myelin protein zero.

We also found significant differences between pairs of morphotypes, especially in terms of molecular signatures/programs. These results reinforce the observations above and show that they are robust to the proportion of fibroblasts and/or epithelial cells present in the compared morphotypes.

In our collection, several cases were represented by several regions and an additional whole-tumor profile. Taking advantage of these matched samples, we investigated several molecular classifiers from an intra-tumor variability perspective as well. The CMS classification was less stable than iCMS, with whole tumor CMS class differing from at least one of the constituent morphological regions in about 60% of cases (11 out of 18, excluding cases in which CMS class was not predicted) (see Figure 5). Additionally, we tested several prognostic and predictive expression-based classifiers/signatures. The goal was not to compare them in terms of their predictive capabilities (the experimental design did not allow for such an exercise), but rather to have a clear picture of the extent to which the various morphotypes “distract” these predictors. We found that all the prognostic signatures varied with the morphological regions with some striking cases in which the morphotype scores exceeded the corresponding whole tumor scores by more than 50%. This observation suggests that, in some cases, the whole tumor-based predictions were too optimistic, the models failing to recognize higher risk cases. While these signatures were derived from whole-tumor expression profiles, their variability across tumor indicates the need for precise tumor sampling strategies.

Our exploratory study has, inherently, several limitations. The selection of cases may not represent the proportions of various morphotypes found in general population of CRC patients. Our selection tried to cover as many scenarios as feasible with a limited number of samples. Also, the tumor heterogeneity in terms of morphotypes cannot be estimated from these data since a single tissue block per tumor was considered. The reduced sample size in some of the paired comparisons within same tumor calls for further external validation. However, our results pave the way to future studies addressing these questions and others related to optimizing the tumor sampling strategy, for example.

We have analyzed the gene expression profiles of six morphotypes (and two peritumoral regions), building a comprehensive molecular picture of their salient features. The observed heterogeneity, especially intra- tumoral (Figure 5), calls for a finer resolution of the tumor sampling in profiling studies. Until spatial transcriptomics becomes integrated in routine clinical practice, using the morphotypes for anchoring the expression profiles is a feasible approach. Our study already provides indications of the molecular programs one would expect to find de-/activated in these regions, thus helping in designing future experiments. The implications for molecular classifiers are clear: it is necessary to account for tumor morphology when designing new biomarkers. Given the sensitivity of many gene-based classifiers to the tumor and stroma proportions in the samples, there is a need to adjust these classifiers to control for their relative proportions. This can be achieved by different means, and we presented an approach based on morphotypes.

From a molecular pathology practice perspective, the molecular descriptors found to vary across morphotypes may help in patient stratification and provide hints for further, more targeted investigations. Several questions call for further investigation: (i) how much of a tumor needs to be embedded to achieve a precise molecular diagnostic? and (ii) what precise tumor region(s) are needed for a molecular diagnostic? The morphotypes selected here may need further refinement and achieving consensus among pathologists regarding their exact definition, a point that could potentially be addressed by automatic image analysis approaches.

Ideas and Speculation

Our analyses indicate that both prognostic and response to therapy signatures may predict more severe cases (shorter relapse free survival or resistance to therapy) when applied to subregions than to the whole tumor. This might be one of the reasons the said signatures may fail their real-world validation. Therefore, morphologically heterogeneous tumors need several sampling locations to provide a more sensible result. Sensitivity and costs analyses need to be performed to estimate the benefits of multi-regional sampling.

Further, the fact that we were able to identify specific molecular programs associated with the morphotypes calls for investigating the inverse problem as well, i.e. whether sufficiently discriminatory features could be extracted for estimating the proportions of the morphotypes from whole tumor profiles.

Materials and Methods

Sample preparation

111 colon cancers (unique patients) were identified in the tumor archive of the Masaryk Memorial Cancer Institute and were assessed by two expert pathologists. Morphological regions of interest, representing complex tubular (CT), desmoplastic (DE), mucinous (MU), papillary (PP), serrated (SE) and solid/trabecular (TB) morphologies, respectively (see Figure 1), were digitally marked in scanned whole slide images (at 20x magnification) and macrodissected for RNA extraction. Additionally, from several slides, tumor- adjacent normal (NR) and tumor-associated stroma (ST). Tumor samples with limited contamination of additional morphologies (<20%) were called “core samples” and used morphotype molecular characterization. The labelling of the regions was repeated after 1 year to ensure a stable assignment. For n=28 cases, whole-tumor regions were macrodissected from the histology section immediately adjacent to the section used for morphological regions. Standard clinical and histopathological variables were retrieved for most of the patients. The study was approved by the ethical committee of Masaryk Memorial Cancer Institute, number 2018/861/MOU.

RNA samples were hybridized on Clariom^TM D (human) microarrays (Thermo Fisher Scientific, Waltham, MA, USA) and data preprocessed using standard Bioconductor (RRID:SCR_006442) packages (see Supplemental Methods).

The data generated in this study are publicly available in ArrayExpress under accession number E-MTAB- 12599.

Bioinformatics analyses

For the identification of differentially expressed genes we used linear models (limma package v.3.52.2) with a cut-off for false discovery rate FDR=0.15. The pathways were scored in terms of enrichment in specific signatures using gene set enrichment analysis (GSEA) (48) as implemented in fgsea package (v.1.22.0). For scoring the signatures in individual samples, we used gene score variation analysis (GSVA) (49) implemented in GSVA package (v.1.44.1). MSigDB (RRID:SCR_016863) (all collections: H, C1-8; v.7.4.1) (50) was used as the main source for gene sets and pathways. Additional cell type-specific gene sets, some derived from whole tumor others from single-cell sequencing studies, representing (i) cancer associated fibroblasts (CAFs) (20,23–25) (ii) epithelial cells (21–23), and (iii) immune cells (20,23) were used (see Supplemental Table 3 for full list). The consensus molecular subtypes were predicted using CMSCaller (51) (v.2.0.1) and the intrinsic epithelial subtypes (9) using the signatures therein (P. Tsantoulis, personal communication, July 2022). The cellular mixture of various tumoral regions was explored computationally using quanTIseq (52) (for immune cells) and ESTIMATE (53) (for tumor purity/epithelial cells). The core samples were used for deriving the lists of differentially expressed genes, for gene set enrichment analyses and for in silico deconvolutions of cell populations. The analyses treating the samples independently, were applied to all samples, including non-core.

Ten different survival/prognostic genomic signatures (full list in Supplemental Table ST13) were computed per-sample as (weighted, when weights were provided) means of signature genes, and 29 sensitivity/resistance signatures selected from MSigDB/C2 were scored by GSVA.

All data analyses were performed in R 4.2. More detailed description of methods can be found in Supplemental material.

Acknowledgements

We would like to thank Dr. Petros Tsantoulis, Hôpitaux Universitaires de Genève, Switzerland, for providing us with a model for classification of intrinsic epithelial subtypes based on gene expression data.

The authors acknowledge funding from the Czech Science Foundation (GAČR) grant no. GA19-08646S. Also, the support of the Research Infrastructure RECETOX RI (No LM2018121) and Cetocoen Plus project (CZ.02.1.01/0.0/0.0/15_003/0000469) financed by the Ministry of Education, Youth and Sports of the Czech Republic (MEYS), the Teaming project (CETOCOEN Excellence 857560; CZ.02.1.01/0.0/0.0/17_043/0009632), and the project National Institute for Cancer Research (Programme EXCELES, ID Project No. LX22NPO5102) - Funded by the European Union - Next Generation EU is acknowledged. This work was supported from the European Union’s Horizon 2020 research and innovation program under grant agreement No 857560. This publication reflects only the author’s view, and the European Commission is not responsible for any use that may be made of the information it contains.

Competing interests

The authors declare no potential conflicts of interest.

Additional information

Funding

The authors acknowledge funding from the Czech Science Foundation (GAČR), grant no. GA19-08646S (EB, MH, VP).