A Clinically Applicable Gene-Expression Classifier Reveals Intrinsic and Extrinsic Contributions to Consensus Molecular Subtypes in Primary and Metastatic Colon Cancer

Colorectal cancer is a leading cause of cancer-related mortality (1–4). Differences in disease progression and response to therapy are frequently observed between individuals, which has often been attributed to disease heterogeneity. Stratification methods have therefore been developed to improve the selection of patients who might benefit from adjuvant therapy, and to identify those who might be at risk for developing distant metastases or have a poor prognosis (3, 5, 6). The main limitation of these approaches is their reliance on histopathologic assessment or evaluation of mutational status or expression of a small number of targets that are unable to capture important molecular features that have the potential to guide therapeutic intervention decisions.

Recent studies have relied on transcriptome-wide data starting with high-quality RNA derived from fresh or frozen tissues to classify tumors into distinct molecular subtypes (7–10). To facilitate clinical translation of colorectal cancer subtyping, an international consortium that included a panel of expert research groups agreed on four consensus molecular subtypes (CMS) of colorectal cancer (11). The CMS1 group was defined by microsatellite instability (MSI), BRAF mutation, and immune infiltration, while CMS2 tumors were dubbed the “canonical subset” and were marked by WNT pathway activation (11, 12). CMS3 malignancies presented with epithelial features and apparent metabolic dysregulation (11, 12). Finally, CMS4 defined a mesenchymal subset that exhibited a prominent stromal component, angiogenesis and transforming growth factor-β (TGFβ) activation (11, 12). In addition, these molecular subsets also present with differences in clinical behavior, including clinical outcome and response to therapy (13), making CMS potentially relevant for diagnosis and disease management. Nevertheless, several unanswered questions hinder the clinical implementation of CMS. For instance, how reliable are CMS assignments in archival tissue, which often consists of low RNA quality samples? Is the heterogeneity identified in primary tumors also represented in metastatic disease? How stable is a given subtype during colorectal cancer progression? In addition, preclinical evaluation of therapeutic modalities aimed at each subtype rests on the availability of relevant in vivo models that faithfully recapitulate the key molecular features of primary and metastatic human disease. Here, we address some of these questions by introducing a customized, colorectal cancer–focused panel based on the NanoString platform that is ideally suited for transcriptional characterization of archival tissues. We demonstrate that the content of our panel can accurately classify colorectal cancer samples from several public data sets into CMS. Using this panel, we then focus on assessing the CMS and associated biological and clinical outcomes of two metastatic colorectal cancer clinical cohorts. We further evaluate the stability of CMS during disease progression by ascribing subtypes in a data set containing primary tumors and patient-matched metastatic disease. Surprisingly, we find only a moderate concordance of subtypes between pairs. We demonstrate that the lack of CMS concordance between matched pairs can be explained by differences in stromal content at the primary and metastatic sites. Finally, we develop a series of murine implantation models that recapitulate several key tumor cell–intrinsic markers of CMS biology identified in primary human colorectal cancers and metastases, making them suitable for in vivo testing of therapeutic modalities aimed each at molecular subtype.

Three clinical cohorts were utilized in this study. These included archival tissues from 105 patients in the DARECK trial, a phase II clinical trial that assessed the activity of anti-EGFR therapy in combination with FOLFIRI in second-line mCRC (NCT01652482), 860 FFPE samples from the Ph3 AVANT trial that examined the efficacy of bevacizumab plus chemotherapy as adjuvant treatment for stage III colorectal cancer (14), and 312 FFPE from 205 colorectal cancer patients obtained from commercial vendors (The MT Group, Inc. and TriStar Technology Group; Supplementary Tables S3 and S5). Paired primary tumors and metastases were available for 107 of 205 subjects from the procured cohort, and clinical information was provided by the commercial vendors (Supplementary Table S5). All samples and clinical information were collected and used after IRB review and approval and in accordance with clinical protocols. Tumors were verified in hematoxylin and eosin (H&E)–stained slides by pathologists. Tumor sections were macrodissected to maximize tumor content before DNA and RNA isolation. Genomic DNA and total RNA were purified using QIAamp DNA FFPE tissue kits (Qiagen) and High Pure FFPE RNA micro-kit (Roche Applied Sciences), respectively (15).

Animal studies were approved by Genentech's Institutional Animal Care and Use Committee and adhere to the NRC Guidelines for the Care and Use of Laboratory Animals. Female NSG mice were purchased from The Jackson Laboratory (8–12 weeks old; stock # 005557).

Mice were anesthetized using isoflurane, injected i.p. with 200 μL of 25 mg/mL d-luciferin (Invitrogen, L2912) and imaged on the Photon Imager (Biospace Lab). During image acquisition, animals received anesthesia from a nose cone delivery system, while their body temperatures were maintained on a thermostatically controlled platform. Photon counts per minute per cm² of observational area were calculated using M3 Vision software (Biospace Lab).

Custom-designed colorectal cancer NanoString Code Sets were used to measure gene expression starting with 250 ng of total RNA from FFPE samples (NanoString Technologies). These largely overlapped panels consisted of 802, 829, and 828 genes for procured, AVANT and DARECK data sets, respectively, and represented key colorectal cancer biology (Supplementary Fig. S1A; Supplementary Table S1). Positive and negative control probes were also included for hybridization efficiency and background calculations. Gene expression was quantified using the NanoString nCounter Analysis System, and raw counts were generated by nSolver software (NanoString). Raw counts were log₂ transformed, and sample-wise normalization was performed by centering each sample on the mean log₂ expression of three housekeeping genes: SP2, TMEM55B, and VPS33B. Subsequently, expression of each gene was mean centered and converted into z-scores. Normalized data for the complete cohort are provided in Supplementary Table S4.

Seven previously published microarray-based colorectal cancer data sets (AMC: GSE33113, ref. 16; NKIAZ: GSE35896, ref. 17; GSE13067 and GSE13294, ref. 18; GSE20916, ref. 19; GSE23878, ref. 20; and GSE37892, ref. 21) with a total of 636 samples were used for training. Testing was performed on 1,088 samples from two data sets (MVRM, ref. 9: GSE17536, GSE17537, GSE14333; CIT, ref. 10: GSE39582). Microarray data for all experiments were obtained in the form of normalized expression values on the log₂ scale and batch corrected using ComBat (22). Subsequently, expression of each gene was mean centered and converted into z-scores. CMS classification for all samples was obtained from Guinney and colleagues (11) based on the consensus calls across the six individual subtyping schemes and was assumed to be the ground truth for evaluation of our NanoString-based classifier. To enable classification of samples that were assayed on the colorectal cancer gene-expression panel, we restricted the microarray data only to genes on the panel. These expression data in conjunction with the CMS classification by Guinney and colleagues (11) were used to train a multinomial classifier based on a generalized linear model via penalized maximum likelihood (23) implemented in the “glmnet” R package. We used misclassification error as the criterion for 10-fold cross-validation and explored the elastic-net penalty (α) from 0 to 1 in increments of 0.1. The lowest misclassification error was achieved for α = 0.2 (Supplementary Fig. S1A), which was used to train the final classifier. Reclassification of AMC, MVRM, and CIT samples as well as novel classification of our NanoString data was performed for the value of the regularization parameter (λ) that gave a minimum mean cross-validated error. R code to perform the CMS classification is available at: https://github.com/rpiskol/CMString.

The status of mutation hotspots covering 13 oncogenes was determined by allele-specific PCR, as described previously (24, 25). MSI status was assessed using MSI Analysis System (cat # MD1641, Promega) with five mononucleotide repeat markers (BAT-25, BAT-26, NR-21, NR-24, and MONO-27) and two pentanucleotide markers for sample identification and contamination (Penta C and Penta D). Samples with one or altered makers were classified as MSI-H, and samples with no altered markers were considered MSI-stable. Penta C and Penta D markers were used to confirm that paired primary and metastatic tumors were from the same patients. The matched identity of 20 paired primary tumors and metastases was further confirmed by a custom, targeted colorectal cancer panel using the Illumina platform (data not shown).

Differential gene-expression analysis was performed using the “limma” package in R (26). To identify genes that were specific to a certain subtype, we contrasted that subtype to all others. We performed Quantitative Set Analysis for Gene Expression (QuSAGE; ref. 27) to identify colorectal cancer subtypes' relevant features. For that purpose, we contrasted primary tumors and metastases from each subtype with samples from all other subtypes. Subsequently, we repeated the analysis using either primary tumors or metastases alone to understand whether the same patterns are observed in both primary tumors and metastases. To determine statistical significance of gene set activity, false discovery rate (FDR)–adjusted P values were obtained by comparing the probability distribution function of log-fold changes in a given gene set to a baseline value of zero using a one-sided test. Only significant results (P < 0.1) are shown. The gene sets used are shown in Supplementary Table S9.

RNA-seq libraries were prepared using 200 ng of FFPE tissue-derived RNA (depending on RNA quality) with TruSeq RNA Access kit (Illumina) by Cofactor Genomics. The libraries were sequenced on Illumina NextSeq500 sequencers to obtain an average of 44 million paired-end reads (75 bp) per sample. Gene-expression quantification from the resulting fastq files was performed using Salmon version 0.6.0 (28) and the human reference transcriptome (NCBI Build 38, annotation release 106). The obtained counts were scaled by library size, quantile-normalized, and precision weights were calculated using the “voom” R package (29). Normalized counts served as the basis for CMS classification using the SSP classifier developed by Guinney and colleagues (11).

BaseScope assays were performed at Advanced Cell Diagnostics. Double Z-probes for AXIN2 (ACD #702191), negative control (bacterial dapB, ACD # 701011), and positive control (PPIB, ACD #701031) were used. In brief, FFPE tissue sections were first treated with heat and protease digestion for target retrieval. Probe hybridization, preamplifier and amplifier steps were performed sequentially, and signals were developed by red chromogenic substrates. Samples were counterstained with Gill's Hematoxylin and analyzed by HALO image analysis software. RNA in situ hybridization (ISH) signals were binned by dots per cell (0 = 0 dot/cell, 1 = 1–3 dots/cell, 2 = 4–9 dots/cell, 3 = 10–15 dots/cell, 4 ≥ 15 dots/cell). A composite H-score for each sample was computed based on the percentage of cells in each bin, as follows: H-score = (0 × % of cells in bin 0) + (1 × % of cells in bin 1) + (2 × % of cells in bin 2) + (3 × % of cells in bin 3) + (4 × % of cells in bin 4). The H-score ranged from 0 to 400. H-scores for tumor and stroma were scored separately.

DLD1, Ls180, and RKO cell lines were obtained from ATCC (CCL-221, CL187, and CRL-2577, respectively). MDST8 was obtained from ECACC (99011801). All of the cell lines were maintained in RPMI-1640 media and supplemented with 10% heat-inactivated FBS (HyClone, SH3007003HI), 1 × GlutaMAX (Gibco, 35050-061), and 1 × penicillin–streptomycin (Gibco, 15140-122). All cell lines were utilized before passage 20 and treated in exponential growth phase at 50% to 75% confluence. All stocks are tested for Mycoplasma prior to and after cells are cryopreserved. Two methods are used to avoid false-positive/negative results: Lonza Mycoalert (Lonza, #LT07-701) and Stratagene Mycosensor (Strategene, #302108).

The lumen implantation procedure has previously been described (30). Briefly, mice were anesthetized by isoflurane inhalation and injected intraperitoneally (i.p.) with buprenorphine at 0.05 to 0.1 mg/kg. A blunt-ended hemostat (Micro-Mosquito, No. 13010-12, Fine Science Tools) was inserted ∼1 cm into the anus. The hemostat was angled toward the mucosa and opened slightly such that a single mucosal fold could be clasped by closing the hemostat to the first notch. The hemostat was retracted from the anus, exposing the clasped exteriorized mucosa. A 10 μL cell suspension containing 50,000 cells admixed with 50% matrigel (Corning) was directly injected into the colonic mucosae. After reversing the prolapse, the hemostat was then released.

Mice were restrained without anesthesia in a conical acrylic restrainer with a heating element to dilate blood vessels. Each animal was injected intravenously in the tail vein with 1.6 mL of the solution containing 50,000 cells in PBS in a single dose administered as a bolus intravenous injection (tail vein) over 4 to 5 seconds (8 seconds maximum). Animals were observed continually for adverse clinical signs for at least 20 minutes after dose.

Total RNA from isolated tumor tissues was extracted using the RNeasy Micro-Kit (Qiagen). Five replicate samples were collected for each treatment condition. The concentration of RNA was determined using NanoDrop 8000 (Thermo Scientific) and the integrity of RNA was determined by Fragment Analyzer (Advanced Analytical Technologies). One microgram of total RNA was used as an input material for library preparation using TruSeq Stranded Total RNA Library Prep Kit (Illumina). The size of the libraries was confirmed using High Sensitivity D1K screen tape (Agilent Technologies), and their concentration was determined by the qPCR-based method using the Library quantification kit (KAPA). The libraries were multiplexed and sequenced on Illumina HiSeq4000 (Illumina). We obtained on average 39 million single-end RNA-seq reads (50 bp) per sample.

Sequencing reads were processed using a XenoFilter pipeline to obtain human- and mouse-specific gene-expression quantifications. Reads were first aligned to ribosomal RNA sequences to remove ribosomal reads. The remaining reads were aligned to the human reference genome (NCBI Build 38) using GSNAP (31) version “2013-10-10,” allowing a maximum of two mismatches per 50 base pair sequence (parameters: “-M 2 -n 10 -B 2 -i 1 -N 1 -w 200000 -E 1 –pairmax-rna = 200000 –clip-overlap”). The same procedure was used to map all reads to the mouse reference genome (NCBI Build 38). Reads that were multimappers in either of the two organisms were removed. Subsequently, concordant and unique reads with ≤ 3 mismatches to the human genome and >3 mismatches to the mouse genome were called human-specific. The reverse mismatch criteria were used to determine mouse-specific reads. Lastly, gene-expression levels were quantified by calculating the number of reads mapped to the exons of each RefSeq gene using the HTSeqGenie R package. Transcript annotation was based on the RefSeq database NCBI Annotation Releases 104 and 106 for mouse and human, respectively.

CMS assignments were performed using the NanoString-based classifier that was developed in this work. The CMS assignment in vitro was described by Linnekamp and colleagues (32).

Total RNA from tumors was isolated using an RNeasy Mini Kit (Qiagen, 74106). cDNA was synthesized using SuperScript IV VILO Master Mix (Invitrogen, 11756050). Quantitative PCR on quadruple samples per condition was carried out in a ViiA7 machine (Applied Biosystems) using the Comparative CT program. Expression of the interested genes was determined relative to the housekeeping gene ACTB. Primer/probe information is listed below: ACTB: Thermo Fisher Scientific Hs01060665_g1; AXIN2: Thermo Fisher Scientific Hs00610344_m1; CDX2: Thermo Fisher Scientific Hs01078080_m1; FRMD6: Thermo Fisher Scientific Hs00378563_m1; ZEB1: Thermo Fisher Scientific Hs01566408_m1.

Overall survival was defined from the time of primary tumor diagnosis to death or the last day of follow-up (censor). P values for Kaplan–Meier analyses were derived using the log-rank test. Kruskal–Wallis test was used to analyze the differences in RNA ISH signals, tumor and stroma contents between CMS2 and CMS4 samples. The hazard ratios between CMSs are calculated by a univariate Cox proportional hazards model in R.

RNA-seq data for human colorectal cancer lines have been deposited to the NCBI GEO database (accession GSE125382).

To enable colorectal cancer subtyping of samples from archival clinical material, we assembled a custom gene-expression panel using the NanoString technology that captures key biological processes involved in colorectal cancer pathogenesis (Supplementary Fig. S1A; Supplementary Table S1; refs. 8–13, 16, 33, 34). The NanoString platform is a well-established methodology for transcriptional profiling of RNA derived from formalin-fixed paraffin-embedded (FFPE) samples (35–38). To generate a NanoString-based classifier, we used nine cohorts of colorectal cancer patients (8, 11, 16) that were partitioned into training (636 patients) and validation sets (1,088 patients; see Material and Methods and Supplementary Table S2). In these cohorts, CMS was previously assigned using the network consensus calls from Guinney and colleagues (8, 11), which we used as a “gold standard” (8, 11). We restricted the expression data to genes that were only present on our NanoString panel, trained a 322-gene classifier that reliably categorized the training set into all four CMS and validated this classifier in the CIT (10) and MVRM (8) data sets (Fig. 1A and B; Supplementary Fig. S1B; Supplementary Table S3). The NanoString-based classifier demonstrated a concordance with the gold-standard transcriptome-derived classifier that reached 90.5% for the CIT and 93.8% for the MVRM cohorts (Supplementary Fig. S1C and S1D; Supplementary Table S3). The classification generated by our NanoString panel retained key CMS clinical features that have been previously described such as the poor outcome of the mesenchymal subset CMS4 (Fig. 1C and D) relative to the other subtypes. To gain additional confidence in our classifier, we subjected 46 high-quality FFPE samples obtained from a procured colorectal cancer cohort to both NanoString and whole-transcriptome RNA-seq and assigned them to CMS using the NanoString and the Guinney classifiers (11), respectively (Supplementary Table S6). In 35 of 41 (85.3%) cases, we observed concordance in CMS assignments made using the NanoString-based classifier and the Guinney classifier (Supplementary Fig. S2A). Furthermore, we detected a high correlation between NanoString- and RNA-seq–based expression values for the majority of the 322 genes present in the classifier [median correlation (rho) = 0.65; Supplementary Fig. S2B; Supplementary Table S7]. These results support the ability of our NanoString classifier to correctly assign CMS onto archival mCRC samples.

To directly evaluate the clinical utility of our NanoString classifier, we classified FFPE colorectal cancers from patients treated with MEHD7945A (39) as part of the DARECK study, a phase II clinical trial that assessed the activity of anti-EGFR therapy in mCRC (NCT01652482; Supplementary Table S3). This cohort comprised 63 of 104 (62%) stage IV patients, RAS wild-type primary tumors, and 41 of 104 (38%) patients who were initially diagnosed with stage I–III colorectal cancer that subsequently metastasized (Supplementary Table S3). All four CMS were detected in that cohort (Fig. 2A), while 4% could not be assigned to a subtype with high confidence. Notably, the prevalence of CMS3 was low in DARECK (3%; Fig. 2A) in comparison with previously published sets (11). Further molecular analysis confirmed that BRAF mutations were enriched in CMS1, as were tumors originating from the right side of the colon (Fig. 2B), both well-described CMS1 attributes (11). Next, we took advantage of the available clinical follow-up to assess the overall survival (OS) and progression-free survival (PFS) for patients in each subtype. Although the mesenchymal CMS4 subtype has been consistently associated with the worst outcome in early-stage colorectal cancer (Fig. 1C and D; ref. 11), we found that patients with CMS1 tumors had significantly worse OS and PFS probabilities than the other CMS (Fig. 2C; Supplementary Tables S11 and S12). As DARECK consisted mainly of metastatic colorectal cancer, we hypothesized that the unfavorable outcome for the CMS1 group in this cohort could be similar to the previously described poor survival after relapse (SAR) from that subset (11). Indeed, like the aggressive CMS4 group, CMS1 trended toward worse SAR within relapsed patients from the combined CIT + MVRM data set, although this did not reach statistical significance due to the small numbers of CMS1 tumors (Supplementary Fig. S3A; Supplementary Tables S11 and S12). For additional validation, we performed NanoString gene expression analysis and subsequent subtype classification of patients from the AVANT trial, a randomized Ph3 study that examined the efficacy of bevacizumab plus chemotherapy as adjuvant treatment for stage III colorectal cancer (14). We analyzed 860 FFPE samples from AVANT and found again that the CMS1 group trended toward worse SAR compared with CMS2 tumors (Supplementary Fig. S3B; Supplementary Tables S11 and S12). These results further support the clinical utility of our NanoString-based classification of colorectal cancer and highlight the aggressive nature of relapsed CMS1 tumors.

Having demonstrated the clinical applicability of NanoString-based CMS classification on archival tissues from clinical trials, we decided to evaluate the stability of CMS during metastatic progression. This is especially relevant since most treatment decisions for patients presenting with metastatic disease are made based on molecular characteristics of their primary resected tumors. We procured tumor samples from a cohort of patients with mCRC at the initial time of diagnosis or presented with stage I–III colorectal cancer that later progressed to metastatic disease (Fig. 3A; Supplementary Fig. S2C; Supplementary Tables S4 and S5). This cohort consisted of samples from 182 primary tumors and 130 metastases. Overall, we found a highly similar proportion of CMS2 and CMS4 in primary tumors and metastases and observed a slightly higher proportion of CMS1 in 22 of 130 (16.9%) metastases compared with 17 of 182 (9.34%) primary colorectal cancers (Fig. 3A and B). We could not assign CMS to 9.34% and 7.7% of the primary and metastatic samples from that cohort, respectively. Intriguingly, we found a substantially lower fraction of CMS3 in 1 of 130 (<1%) in metastases compared with 19 of 182 (11%) in primary tumors. The lack of CMS3 patients was seen despite a high prevalence of KRAS mutations in these samples, which is known to strongly associate with that subtype (11). To further corroborate the observed paucity of CMS3 in metastatic disease, we assigned CMS in a publicly available cohort that comprised liver metastasis (40) and primary tumors, in which gene-expression profiles were generated using Affymetrix microarrays. Classification of these metastatic samples using either the SSP classifier (11) or our NanoString classifier revealed a depletion of CMS3 in liver metastasis when compared with primary tumors, supporting our initial observation (data not shown). Similar to the results obtained from the DARECK trial, CMS1 tumors in the procured cohort presented with worse OS (Fig. 3C). This further reinforces the poor outcome of CMS1 in metastatic disease. Finally, we found a moderate association between CMS and the site of metastasis. For example, CMS1 were mainly found to disseminate to the liver and lung while non-liver/lung metastasis were enriched in CMS4 tumors (Supplementary Fig. S2C). These preliminary data potentially reflect the different tropism of distinct CMS, although they require further validation given the small numbers of samples, potential CMS misclassification, and possible population selection bias in our data set.

We next compared the biological attributes of the two most prevalent CMS2 and CMS4 subsets between primary tumors and metastasis. First, we evaluated the desmoplastic stroma content, a defining feature of the aggressive CMS4 tumors. Histologic examination of tumor/stroma content revealed a significantly higher percentage of tumor cells and conversely a lower fraction of stromal cells in CMS2 compared with CMS4 in primary malignancies (Fig. 3D). This tumor/stroma ratio was maintained in metastases. To bolster our histologic data, we examined all genes on our panel that were differentially expressed between CMS1, 2, and 4 in primary tumors and, separately, in metastases. We observed a significant overlap of differentially expressed genes between primaries and metastases for all CMS groups except CMS1 (Supplementary Fig. S4; Supplementary Table S8). Next, we examined gene pathway activity with QuSAGE using pathway gene signatures and well-characterized cellular processes (Supplementary Table S9; ref. 11). Consistent with previous reports (11, 12), we detected expression of immune markers, including PDL1, in CMS1 tumors (Fig. 3E). In line with our histologic data, CMS2 primary tumors displayed upregulation of MYC as well as WNT activation signatures (Fig. 3E), while CMS4 primary tumors showed elevated expression levels of genes related to stromal biology, TGFB, EMT, and angiogenesis (Fig. 3E). Interestingly, primary tumors from the CMS4 group also exhibited upregulation of PDL1 and an immune infiltration signature (Fig. 3E). More specifically, CMS4 primary tumors displayed an increase in the expression of genes involved in T-cell regulation as well as TH17 and TH1 activation, which were not shared with CMS1 malignancies (Fig. 3E). Importantly, identical patterns of up/downregulation of CMS-related genes in primary tumors were also observed in metastases (Fig. 3E). This demonstrates that the biological properties of CMS2 and 4, such as WNT activity and stroma content, respectively, are similar in colon tumors located at the primary or metastatic site.

To test the concordance and stability of CMS during metastatic progression, we used the subset of 71 patients from the procured data set in which both primary tumor and metastasis was obtained. These 71 pairs were confirmed to originate from the same patients by short tandem repeat (STR) analysis (see Materials and Methods). We observed a concordance of CMS assignments in 37 of 62 (60%) primary/metastasis pairs (Fig. 4A). The sizable number of discordant cases directly suggests that CMS identification in primary tumors cannot be readily used to infer the disease subtype in the corresponding metastatic lesion.

To probe the underlying cause of CMS discordance between primary tumors and metastases, we performed additional molecular analysis on a subset of paired samples. First, we evaluated whether clonal drift and selection had occurred and could explain the shift in subtype identity. Next-generation sequencing of a panel of colorectal cancer relevant mutations revealed an overall strong overlap in allelic variants between matched primary and metastasis pairs (Supplementary Fig. S5A), although the number of private mutations between primary tumors and metastases trended to be higher in the discordant compared with concordant cases (Supplementary Fig. S5B). Discordant pairs had, in general, lower classification probabilities than the concordant pairs (Supplementary Fig. S6A), but this could not fully explain CMS switching. As an example, CMS4 tumors that were classified as CMS2 in the metastasis had high probabilities of CMS assignment (Supplementary Fig. S6B). This suggests that the lack of concordance is not due to a high misclassification error rate alone, but, rather, it could be related to a biological CMS switch. To investigate the biological nature of the CMS switch, we studied gene-expression changes that were associated with discordant cases. We found that discordance was strongly correlated with a change in gene-expression signatures known to define a given subset. For instance, CMS2 primary tumors that shifted to CMS4 in metastasis displayed a shift in epithelial gene expression toward EMT, stromal infiltration and TGF beta signaling, while the reciprocal was observed for CMS4 to CMS2 shifters (Fig. 4B; Supplementary Fig. S6C). Furthermore, tumor epithelial content was high in CMS2 concordant cases but decreased in liver metastasis that shifted from CMS2 in the primary to CMS4 (Fig. 4C).

Because these data suggest that changes in environment composition (e.g., stromal content) could explain the discordance observed, we decided to evaluate whether more intrinsic tumor cell features were also changed in liver metastasis compared with their primary counterpart. First, we assessed the spatial distribution of AXIN2, a canonical WNT pathway target gene that is highly expressed in the epithelial subtype CMS2 and relatively lower in CMS4. ISH revealed intense AXIN2 staining in the epithelial compartment of primary CMS2 tumors and a much lower signal in epithelial cells from CMS4 tumors (Fig. 4D and E). This staining pattern was similar in primary colorectal cancers and in metastases, suggesting that WNT activity is at least partly tumor intrinsic in CMS2 (32). Second, we generated an experimentally tractable system, by implanting luciferase tagged human colorectal cancer lines encompassing all CMS groups (refs. 32, 41; Supplementary Fig. S7A) in vivo into immunocompromised animals (Fig. 4F). Of note, we could not recover luciferase tagged MDST8 after viral transduction, preventing our use of this line for live imaging. Orthotopic implantation of the lines in the colon through induced rectal prolapse or in the liver through a modified hydrodynamic tail-vein injection (Fig. 4F and G) yielded tumor growth at each site, albeit with variable take rates (Fig. 4G; Supplementary Fig. S7B and S7C). Next, we isolated several primary colon tumors and multiple liver metastases from independent animals and subjected them to whole-transcriptome analysis. By restricting our analysis to human-specific reads we directly measured tumor intrinsic gene expression. Principal component analysis using all expressed genes revealed a strong separation of the different lines, which clustered together irrespective of tissue location (i.e., colon or liver; Supplementary Fig. S7D). Furthermore, CMS classification using our NanoString gene classifier showed a high concordance of CMS calls between primary and metastatic tumors (Fig. 4H). Interestingly, although Ls180 was initially classified in vitro as CMS3 by Linnekamp and colleagues (32), that line was consistently ascribed to CMS2 in both primary tumors and liver metastasis. These data suggest that tumor-intrinsic gene expression of CMS is maintained between primary tumors and liver metastasis in this setting. Altogether, we conclude that tumor-intrinsic features, including genetic alterations and tumor-specific gene expression, are maintained during colon cancer progression but changes in environment composition at the primary and metastatic sites can have an impact on accurate CMS classification.

Recent efforts have led to the stratification of colorectal cancer patients into distinct molecular subtypes (7, 8, 10, 42, 43). Such molecular classification has been achieved primarily through whole-transcriptome profiling, where the requirements for high-quality RNA or sufficient bulk tumor can hamper widespread clinical adoption. Therefore, there is a need for rapid and cost-effective alternatives to whole genome expression profiling in both the clinical setting and in translational research. Preferentially, such tools should enable the identification of patient subtypes to guide subsequent treatment, allow access to large patient cohorts from clinical trials for which only FFPE tissue is available, and retain sufficient biological information to inform on subtype-specific molecular processes.

With these criteria in mind, we have developed a customized colorectal cancer NanoString panel that is cost effective and ideally suited for gene-expression analysis of archival tissues. The portability of the method was highlighted by the successful classification of patients into distinct molecular subtypes in three independent clinical colorectal cancer cohorts, including metastatic resected specimens. In these data sets, we found that poor prognosis was consistently endowed within the CMS1 inflammatory-like subtype in the DARECK trial and vendor-procured metastatic cohorts, in contrast to the favorable clinical performance of this subtype in the nonmetastatic setting. We hypothesize that the poor prognosis portended by CMS1 in mCRC is consistent with the poor SAR of that subset in patients with nonmetastatic disease. Although the underlying molecular foundation for the poor outcome in metastatic CMS1 patients remains unknown, it may be explained, in part, by the enrichment of mutations in BRAF or the high prevalence of microsatellite instability, molecular features both previously associated with a dismal prognosis in relapsed colorectal cancer (44). Importantly, we provide evidence to suggest that the metabolic CMS3 is strongly underrepresented in advanced metastatic disease, even in data sets that have a selection bias toward KRAS-mutant colorectal cancer. RAS wild-type status in the DARECK trial may influence CMS distribution and, potentially, lead to underrepresentation of CMS3 in this cohort. It is tempting to speculate that the liver environment influences expression of key colorectal cancer metabolic genes that are required for accurate classification of that subset.

Our classifier demonstrated an overall 85% concordance with RNA-seq–based classification using the genes that comprise the Guinney classifier, which is in line with previous work comparing custom NanoString panels with other transcriptional platforms (45, 46). Although promising, the clinical implementation of our panel will require an improvement of CMS classification consistency, which will be aided by standardizing processing guidelines and optimization of the gene panel. Unfortunately, the lack of fresh frozen (FF) specimen from our cohort precluded a direct cross comparison of our panel on paired FF and FFPE tissues. This is an important limitation of the current study since gene-expression measurements can significantly vary between FF and FFPE samples. In addition, because the “gold standard” CMS assignment as defined by Guinney and colleagues were made on FF samples, further studies aimed at evaluating our NanoString panel in paired FF/FFPE should be considered in order to more critically assess the accuracy of our classifier.

Despite these limitations, we have further leveraged our NanoString panel to provide critical information on the stability of CMS during colorectal cancer metastatic progression. Although we observed that many of the biological attributes of primary colorectal cancer are intrinsic, we also demonstrated that CMS ascribed at the primary site does not always recapitulate in disseminated disease. We attributed this lack of concordance primarily to distinct environmental composition at primary and metastatic sites. More specifically, we showed that CMS2/4 discordance between matched pairs can be attributed, at least in part, to differences in the stromal composition of primary colorectal cancers and metastases. To support this conclusion, we further demonstrated that human cell lines implanted in the colon or liver of mice have high CMS concordance between the two sites, further suggesting that tumor-intrinsic gene expression is largely maintained during metastatic progression. Our observation has direct clinical impact as it cautions that the diagnosis made on a resected primary tumor may be inadequate to guide subtype-specific treatment of the corresponding metastasis.

Importantly, our efforts have also led to the generation of in vivo models that closely recapitulate key aspects of the human disease. These can be useful tools to further evaluate CMS-specific therapeutic modalities in vivo in relevant tissue sites and will hopefully aid at improving the clinical outcome of this disease.

R. Piskol holds ownership interest (including patents) in Roche. R. Patel holds ownership interest (including patents) in Roche. X. Qu holds ownership interest (including patents) in Roche. H. Koeppen holds ownership interest (including patents) in Roche. F. de Sauvage holds ownership interest (including patents) in Roche. M. R. Lackner holds ownership interest (including patents) in Roche. F. de Sousa e Melo holds ownership interest (including patents) in Roche. No potential conflicts of interest were disclosed by the other authors.

Conception and design: R. Piskol, L. Huw, D. Kim, M.R. Lackner, F. de Sousa e Melo, O. Kabbarah

Development of methodology: R. Piskol, L. Huw, D. Kim, R. Patel, F. de Sousa e Melo, O. Kabbarah

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): L. Huw, I. Sergin, Z. Modrusan, D. Kim, R. Tam, J. Burton, X. Qu, H. Koeppen, T. Sumiyoshi, F. de Sousa e Melo, O. Kabbarah

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): R. Piskol, L. Huw, I. Sergin, E. Penuel, H. Koeppen, M.R. Lackner, F. de Sousa e Melo, O. Kabbarah

Writing, review, and/or revision of the manuscript: R. Piskol, L. Huw, C. Kljin, R. Tam, E. Penuel, H. Koeppen, F. de Sauvage, M.R. Lackner, F. de Sousa e Melo, O. Kabbarah

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): N. Kljavin, R. Patel, O. Kabbarah

Study supervision: F. de Sauvage, F. de Sousa e Melo, O. Kabbarah

We would like to thank Gregoire Pau for the development of the XenoFilter procedure for the deconvolution of human and mouse reads and Matthew Chang for code review. We thank Meghna Das Takur for providing assistance with the AVANT data set. Weilan Ye and Aditya Murthi for critical comments on the manuscript.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.