Heterogeneity of colorectal carcinoma (CRC) represents a major hurdle towards personalized medicine. Efforts based on whole tumor profiling demonstrated that the CRC molecular subtypes were associated with specific tumor morphological patterns representing tumor subregions. We hypothesize that whole- tumor molecular descriptors depend on the morphological heterogeneity with significant impact on current molecular predictors.
We investigated intra-tumor heterogeneity by morphology-guided transcriptomics to better understand the links between gene expression and tumor morphology represented by six morphological patterns (morphotypes): complex tubular, desmoplastic, mucinous, papillary, serrated, and solid/trabecular. Whole-transcriptome profiling by microarrays of 202 tumor regions (morphotypes, tumor-adjacent normal tissue, supportive stroma, and matched whole tumors) from 111 stage II-IV CRCs identified morphotype-specific gene expression profiles and molecular programs and differences in their cellular buildup. The proportion of cell types (fibroblasts, epithelial and immune cells) and differentiation of epithelial cells were the main drivers of the observed disparities with activation of EMT and TNF-α signaling in contrast to MYC and E2F targets signaling, defining major gradients of changes at molecular level. Several gene expression-based (including single-cell) classifiers, prognostic and predictive signatures were examined to study their behavior across morphotypes. Most exhibited important morphotype-dependent variability within same tumor sections, with regional predictions often contradicting the whole-tumor classification.
The results show that morphotype-based tumor sampling allows the detection of molecular features that would otherwise be distilled in whole tumor profile, while maintaining histopathology context for their interpretation. This represents a practical approach at improving the reproducibility of expression profiling and, by consequence, of gene-based classifiers.
This study presents a valuable finding on the putative molecular patterns underlying characteristic morphological regions observed in colorectal cancer (CRC). The authors provide a morphological framework through which clinicians might improve the performance of molecular signatures and consequently predict the clinical response of patients with better accuracy. The evidence supporting the claims of the authors is solid. The work will be of interest to clinicians and cancer biologists working in the field of CRC.
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.
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.
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).
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.
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).
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).
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).
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
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 ClariomTM 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.
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.
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.
The authors declare no potential conflicts of interest.
The authors acknowledge funding from the Czech Science Foundation (GAČR), grant no. GA19-08646S (EB, MH, VP).
- 1.Amin MB, American Joint Committee on Cancer, American Cancer Society, editors. AJCC cancer staging manual. Eight edition / editor-in-chief, Mahul B. Amin, MD, FCAP ; editors, Stephen B. Edge, MD, FACS [and 16 others] ; Donna M. Gress, RHIT, CTR-Technical editor ; Laura R. Meyer, CAPM- Managing editor. Chicago IL: American Joint Committee on Cancer, Springer; 2017. 1024 p.
- 2.Colon cancer molecular subtypes identified by expression profiling and associated to stroma, mucinous type and different clinical behaviorBMC Cancer 12
- 3.Gene expression patterns unveil a new level of molecular heterogeneity in colorectal cancer: Gene expression heterogeneity in colorectal cancerJ Pathol 231:63–76
- 4.Gene Expression Classification of Colon Cancer into Molecular Subtypes: Characterization, Validation, and Prognostic ValuePLoS Med 10
- 5.Poor-prognosis colon cancer is defined by a molecularly distinct subtype and develops from serrated precursor lesionsNat Med 19:614–8
- 6.A colorectal cancer classification system that associates cellular phenotype and responses to therapyNat Med 19:619–25
- 7.Colorectal cancer intrinsic subtypes predict chemotherapy benefit, deficient mismatch repair and epithelial-to-mesenchymal transition: Molecular subtypes in colorectal cancerInt J Cancer 134:552–62
- 8.The consensus molecular subtypes of colorectal cancerNat Med 21:1350–6
- 9.Single-cell and bulk transcriptome sequencing identifies two epithelial tumor cell states and refines the consensus molecular classification of colorectal cancerNat Genet 54:963–75
- 10.Comprehensive molecular characterization of human colon and rectal cancerNature 487
- 11.The single-cell sequencing: new developments and medical applicationsCell & Bioscience 9
- 12.Exploring tissue architecture using spatial transcriptomicsNature
- 13.Classification of colorectal cancer based on correlation of clinical, morphological and molecular featuresHistopathology 50:113–30
- 14.Experiments in molecular subtype recognition based on histopathology imagesIn: 2016 IEEE 13th International Symposium on Biomedical Imaging (ISBI) [Internet]
- 15.Image-based surrogate biomarkers for molecular subtypes of colorectal cancerBioinformatics 33:2002–9
- 16.Molecular pathological classification of colorectal cancerVirchows Arch 469:125–34
- 17., International Agency for Research on Cancer, editors. WHO classification of tumours of the digestive system. 4th ed. Lyon: International Agency for Research on Cancer; 2010. 417 p. (World Health Organization classification of tumours)
- 18.Challenging the Cancer Molecular Stratification Dogma: Intratumoral Heterogeneity Undermines Consensus Molecular Subtypes and Potential Diagnostic Value in Colorectal CancerClin Cancer Res
- 19.Standardising RNA profiling based biomarker application in cancer—The need for robust control of technical variablesBiochimica et Biophysica Acta (BBA) - Reviews on Cancer 1868:258–72
- 20.Stromal contribution to the colorectal cancer transcriptomeNat Genet 47:312–9
- 21.Gene expression patterns of human colon tops and basal crypts and BMP antagonists as intestinal stem cell niche factorsProceedings of the National Academy of Sciences 104:15418–23
- 22.The Intestinal Stem Cell Signature Identifies Colorectal Cancer Stem Cells and Predicts Disease RelapseCell Stem Cell 8:511–24
- 23.Spatially organized multicellular immune hubs in human colorectal cancerCell 184:4734–4752
- 24.Refining colorectal cancer classification and clinical stratification through a single-cell atlasGenome Biology 23
- 25.Single-Cell Analysis Reveals Fibroblast Clusters Linked to Immunotherapy Resistance in CancerCancer Discovery 10:1330–51
- 26.Cancer-associated fibroblasts: Key players in shaping the tumor immune microenvironmentImmunological Reviews 302:241–58
- 27.Comment on “Identification of EMT-related high- risk stage II colorectal cancer and characterisation of metastasis-related genes”cited
- 28.Histological phenotypic subtypes predict recurrence risk and response to adjuvant chemotherapy in patients with stage III colorectal cancerThe Journal of Pathology: Clinical Research 6:283–96
- 29.Clinical Value of Consensus Molecular Subtypes in Colorectal Cancer: A Systematic Review and Meta-AnalysisJNCI: Journal of the National Cancer Institute 114:503–16
- 30.The Molecular Hallmarks of the Serrated Pathway in Colorectal CancerCancers 11
- 31.Circular RNA PIP5K1A promotes colon cancer development through inhibiting miR-1273aWorld Journal of Gastroenterology 25:5300–9
- 32.Prognostic Significance and Immunological Role of FBXO5 in Human Cancers: A Systematic Pan-Cancer Analysis
- 33.TGF-β-Induced FLRT3 Attenuation Is Essential for Cancer- Associated Fibroblast–Mediated Epithelial–Mesenchymal Transition in Colorectal CancerMolecular Cancer Research 20:1247–59
- 34.Direct reprogramming of fibroblasts into endothelial cells capable of angiogenesis and reendothelialization in tissue- engineered vesselsProceedings of the National Academy of Sciences 109:13793–8
- 35.E2F7, a novel target, is up-regulated by p53 and mediates DNA damage-dependent transcriptional repressionGenes Dev 26:1533–45
- 36.The multifarious roles of the chemokine CXCL14 in cancer progression and immune responsesMolecular Carcinogenesis 59:794–806
- 37.A Combined ULBP2 and SEMA5A Expression Signature as a Prognostic and Predictive Biomarker for Colon CancerJournal of Cancer 8:1113–22
- 38.Fibroblast growth factor 7 signalling is disrupted in colorectal cancer and is a potential marker of field cancerisationJournal of Gastrointestinal Oncology [Internet
- 39.MUC1, MUC2, MUC5AC, and MUC6 in colorectal cancer: expression profiles and clinical significanceVirchows Arch 469
- 40.Prognostic value of desmoplastic reaction characterisation in stage II colon cancer: prospective validation in a Phase 3 study (SACURA Trial)British Journal of Cancer 124:1088–97
- 41.The function of Piezo1 in colon cancer metastasis and its potential regulatory mechanismJ Cancer Res Clin Oncol 146:1139–52
- 42.Slit/Robo pathway: a promising therapeutic target for cancerDrug Discovery Today 20:156–64
- 43.Pan-Cancer Analysis of OLFML2B Expression and Its Association With Prognosis and Immune Infiltration
- 44.DPYD Variants as Predictors of 5- fluorouracil Toxicity in Adjuvant Colon Cancer Treatment (NCCTG N0147)JNCI: Journal of the National Cancer Institute 106
- 45.Abrogation of Gli3 expression suppresses the growth of colon cancer cells via activation of p53Experimental Cell Research 318:539–49
- 46.HOXA9 Overexpression Contributes to Stem Cell Overpopulation That Drives Development and Growth of Colorectal CancerInternational Journal of Molecular Sciences 23
- 47.A DNA Repair Pathway–Focused Score for Prediction of Outcomes in Ovarian Cancer Treated With Platinum-Based ChemotherapyJNCI: Journal of the National Cancer Institute 104:670–81
- 48.Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profilesProceedings of the National Academy of Sciences 102:15545–50
- 49.GSVA: gene set variation analysis for microarray and RNA-Seq dataBMC Bioinformatics 14
- 50.The Molecular Signatures Database Hallmark Gene Set CollectionCell Systems 1:417–25
- 51.CMScaller: an R package for consensus molecular subtyping of colorectal cancer pre-clinical modelsSci Rep 7
- 52.Molecular and pharmacological modulators of the tumor immune contexture revealed by deconvolution of RNA-seq dataGenome Medicine 11
- 53.Inferring tumour purity and stromal and immune cell admixture from expression dataNat Commun 4