MicroRNA Classifier and Nomogram for Metastasis Prediction in Colon Cancer

miRNAs are endogenous noncoding RNAs with a sequence length of 19 to 25 nucleotides (miRBase v19), that are key posttranscriptional regulators of gene expression (1), and highly conserved between species (2). They have important regulatory functions in basic cellular processes, thereby affecting major biologic systems such as stemness and immunity (3). Depending on the level of complementarity of the miRNA to its target mRNA, the particular transcripts are targeted for transcriptional repression or cleaved for degradation. miRNAs are predicted to be involved in regulation of more than 60% of the protein-coding RNAs (4). With their frequent location at fragile sites and genomic regions involved in cancer (5), it is not surprising that deregulation of miRNA expression, by both genetic and epigenetic mechanisms, is linked to tumorigenesis (6, 7), including the formation of metastasis (8).

Treatment modalities of patients with colon cancer are still essentially based on the anatomic extent of disease at diagnosis (TNM staging method of the American Joint Committee on Cancer/UICC). Even under the current classification, patient groups still show large differences in clinical outcome (9–11), especially for patients with locally advanced, but not distantly metastasized tumors (TNM stage II–III). Clinically, current criteria for identification of high-risk patients in this group of colon cancer, which would require adjuvant therapy, include clinical and histologic criteria T4 stage, poor differentiation, tumor perforation, or insufficient lymphadenectomy (9, 12). Biomarkers that can predict the development of distant recurrent disease would help identifying those patients with cancer that may benefit from adjuvant therapy and/or close follow-up in addition to surgery. The use of such biomarkers in conjunction with the TNM staging method, by means of a nomogram, may contribute to an approach of personalized treatment based on individual tumor characteristics. miRNA profiles, identified by comparing expression of miRNAs between patients with and without metastases in the follow-up, have been constructed for several cancer types (13). As discussed by de Krijger and colleagues (14), a similar profile for patients with colon cancer is warranted as the development of distant metastases mainly accounts for the high lethality rates of colon cancer.

Single miRNAs have multiple target genes (15), and target genes can be regulated by different miRNA transcripts. The tissue-specific (16) and cell type–specific (17) nature of miRNA expression, combined with the high degree of heterogeneity of colon tumors with regard to the amount of stroma and immune infiltrate, requires selection of target cells by microdissection (18), like laser capture microdissection (LCM). To be able to investigate expression of the whole micronome, and the possibility of detecting novel miRNAs, a next-generation sequencing (NGS) approach is recommended (19). Previous studies reporting on miRNA classifier profiles for metastasis prediction in colon or colorectal cancer (20–24) did not use the optimal approach of identification of miRNAs by comparison of their expression patterns between patients with and without metachronous metastasis using NGS in microdissected tumors. For PCR-based relative quantification of miRNA expression, a verified standard for normalization of the amount of input material is lacking. We therefore also aimed to assess which miRNA(s) would be suitable for usage as reference in miRNA-specific real-time PCR.

In this bicentric international study, we studied miRNA expression in two independent retrospective patient cohorts, both containing microsatellite stable tumors. Tumor tissue from the discovery cohort (Dutch colorectal cancer cohort, N = 50) and the validation cohort (German colon cancer cohort, N = 43) were subjected to LCM. High-depth NGS was used to study the predictive value of miRNAs for metachronous metastasis risk in the discovery cohort. NGS data were technically validated with low-depth NGS and miRNA-specific stem–loop real-time quantitative PCRs (qRT-PCR). Results were validated in an independent cohort, using qRT-PCR. Multivariate survival analyses for distant recurrence-free period showed a HR of nearly 11 in the validation cohort. With our study, we identified and validated a miRNA expression classifier in the primary tumor, and proposed a nomogram capable of predicting metastasis risk in microsatellite stable colon cancer. In addition, we recommend a new standard for normalization of the amount of input material in PCR-based relative quantification of colon tissue–derived miRNA expression.

Patients for the discovery phase were selected from a colorectal cancer patient population (N = 200) who underwent surgery between 2001 and 2012 at the Leiden University Medical Center (Leiden, the Netherlands), with available snap-frozen primary tumor material. Microsatellite stability status was determined as described previously (25, 26). For the discovery study cohort, we included only histologically proven TNM stage II–III first primary colorectal carcinomas of the microsatellite stable phenotype, of patients who did not receive preoperative treatment. We selected 24 patients with distant metachronous metastases and matched these with patients without metastases. Matching criteria were TNM stage, age at surgery, gender, and location of the tumor. The discovery study cohort consisted of 50 patients (Table 1). All patients were under clinical follow-up surveillance according to the current Dutch National Guidelines.

For validation, patients were selected from a colon cancer patient population (N = 242), in part described previously (27), who underwent surgery between 1987 and 2006 at the Department of Surgery, Klinikum rechts der Isar, Technische Universität München (TUM, Munich, Germany). All patients gave written, informed consent, with approval of the Ethics committee of the Medical Faculty of TUM (#1926/07). Included tumors were histologically proven TNM stage II–III first primary colon carcinomas with microsatellite stable phenotype, with available snap-frozen primary tumor material. We selected the patients with metachronous distant metastases (N = 21) and matched these with patients without metastases. Matching criteria were similar as described above. The validation study population consisted of 43 patients (Table 1). All patients received the recommended postoperative follow-up care, according to the national German guidelines (S3 guidelines by the German Cancer Society DGK/AWMF).

All samples were coded, according to National ethical guidelines (“Code for Proper Secondary Use of Human Tissue,” Dutch Federation of Medical Scientific Societies). This study was performed according to the REMARK guidelines (NCI-EORTC; ref. 28).

Snap-frozen tumor tissues were collected for each patient. LCM was performed to select cancerous epithelial cells. Hematoxylin staining, LCM, and subsequent sample lysis was completed within 1 hour. Sections (8 μm thick) mounted on LCM-suitable slides (Zeiss) were processed on the PALM MicroBeam (Zeiss). A minimum of 3-mm² cancer cell area of the sections was collected. RNA was extracted using the Just-a-Tube sample total RNA/microRNA purification kit (Charm Biotech) for the discovery cohort or TRIzol Reagent (Life Technologies) for the validation cohort, according to manufacturers' recommendations. Total RNA quantity and miRNA quality were assessed with the NanoDrop 1000 (Thermo Scientific) and RNA 6000 Pico total RNA-Chip (Agilent Technologies), respectively.

Tumor tissue libraries for NGS were prepared as described by Buermans and colleagues (29). One universal forward primer and 14 reverse primers, with different barcodes, were available for the amplification PCR. A unique barcode was assigned to each sample within a sample pool. The presence of cDNA was assessed with a high-sensitivity DNA Chip (HS-DNA Chip; Agilent Technologies). In addition, a set of qRT-PCRs served to show whether miRNAs were present. Samples with at least 90% of desirable sequence lengths, with calculated insert sizes of 17–25 nucleotides, were randomly pooled into 4 pools in equimolar volumes of 2 nmol/L and subsequently reassessed with HS-DNA Chip.

High-depth NGS was performed on the four miRNA pools, loaded into 4 different lanes of a Illumina HiSeq2000 v3-flowcell with v3-reagents (Illumina). The sequence run type was a single-end read of 50 cycles with an added sample-index read for pooled-sample identification. In addition, low-depth NGS was performed using the Ion Torrent platform with an Ion314 chip (Life Technologies) to sequence a pooled subset of 13 of the 50 samples, and used for technical validation. Data from both NGS runs were extracted and processed into raw count data per miRNA as described previously (29, 30), with minor adaptations. In summary, the 3′-library adapter was removed from the sequences, and remaining reads of at least 15 bases were aligned to the human Hg19 reference genome using the Bowtie short-read aligner allowing for up to three loci with two mismatches. Crossmapping correction was applied to the alignments as described by de Hoon and colleagues (31). Aligned sequence tags were annotated on the basis of the miRBase v19 release resulting in separate expression values for the 5′ and 3′ mature miRNAs from the hairpin precursors. Data are publicly available (GEO number: GSE63119).

Data analysis was performed in R version 3.0.1. Library sizes were defined as the total number of annotated miRNA reads per sample. Normalization started with transforming raw miRNA count data to counts per million (cpm) using the library size consisting of only annotated miRNA isoforms (miRBase v19). Subsequently, quantile normalization (R-package limma) was performed, according to the recommendations of Garmire and colleagues (32, 33). All libraries passed quality control (R-package arrayQualityMetrics) before and after normalization. Low-expression miRNAs, which are those without expression of at least 1 cpm in at least 10 of 50 samples, were excluded from analysis. Subsequently library sizes were reset to the remaining number of miRNA reads.

Differential expression analysis was performed using the nonparametric Wilcoxon rank sum test (R-package stats) to restrict the influence of expression outliers on the results, and the number of false positive miRNAs. The number of normalized reads is representative for the expression of a miRNA in that sequence library. Multiple testing correction of P values was performed with false discovery rate (fdr).

Technical validation of the HiSeq NGS data with the Ion Torrent NGS data was performed by determining three different correlation coefficients. The intraclass correlation coefficient (ICC; absolute agreement) as a measure of intra- and intersamples variability and the Pearson product–moment (Pearson r) correlation coefficient as a measure of the linear correlation between the two datasets per sample were determined on cpm data of all detected miRNAs. In addition, the Spearman rank correlation coefficient (Spearman ρ) on rank data was calculated as a measure of statistical dependence between the two datasets per sample. Included were ranked data of miRNAs with more than 2 raw counts in both NGS datasets. A significant P value (P < 0.05) indicated correlation for the miRNA expression profiles per patient for the two datasets.

Identification of reference miRNAs for technical validation and validation of results by qPCR were determined using the coefficient of variation (CV) per miRNA, with a low CV indicating the miRNAs with the best linear relation between the number of isoform-specific miRNA reads and the total number of reads across the samples. To calculate the CV per miRNA, the SD was divided by the average number of reads.

Eight different genes from the top ten list of the differential expression analysis and two reference miRNAs were validated by qRT-PCR. miRNA-specific cDNA was prepared with the TaqMan miRNA Reverse Transcription Kit (Applied Biosystems) combined with miRNA assays (Applied Biosystems) in 96-well plates. The real-time PCR of the different TaqMan miRNA assays (Applied Biosystems) were run on the CFX-384 machine (Bio-Rad), using the TaqMan Universal Master Mix II (Applied Biosystems). Thresholds for calculation of the cycle threshold (C_t) were set separately for each miRNA. Data are available upon request. Triplicate measurements were averaged for further analyses. The average of miR16 C_t value and miR26a C_t value per sample served as reference for relative expression calculations (formula |$2^{\Delta {C_{\rm t}}$|⁠). The Pearson's r, as a measure of the linear correlation between the two datasets per sample, was determined per miRNA on NGS and qRT-PCR data of all patients. In addition, the Spearman's ρ on rank data, as a measure of statistical dependence between the two datasets, were calculated per miRNA on NGS and qRT-PCR data of all patients.

Relative expression of different miRNAs, as determined by qRT-PCR, was used as predictor in survival analyses with the statistical package SPSS 20.0 for Windows (SPSS Inc). Patients were dichotomized per miRNA based on the average of mean and median as determined in the discovery cohort. Receiver operating characteristic (ROC) curves for the individual miRNAs were used to assess the specificity and sensitivity of the chosen cut-off values. The relationship between dichotomized relative miRNA expression and known clinicopathologic parameters was investigated using the Pearson χ² test. The Cox proportional hazards model was used to analyze the association between relative miRNA expression with distant recurrence-free period (DRFP), defined as the time from surgery until the diagnosis of the first distant recurrence. Differences in DRFP between patient groups are presented as HRs. Cumulative risk curves were calculated, accounting for death by all causes as competing risk for metastasis, and served to visualize the associations (34). The most important covariates for colon cancer, being age at surgery, TNM stage, and gender, were included in multivariate analyses irrespective of statistical significance. The high expression group served as reference for all survival analyses. All tests were two-tailed and P values <0.05 were considered statistically significant. A nomogram for prediction of distant metastasis in colon cancer was constructed (R-package rms) using the model obtained from the multivariate survival analysis of the discovery cohort. Calibration plots (R-package rms) and ROC curves were used to assess the prediction performance and accuracy of the nomogram. For comparison, we performed the survival analyses for the validation cohort after building the survival model by using the dichotomization cutoff as determined on the discovery cohort.

miRNA gene target prediction was performed by combining results from three different target prediction databases and a database containing validated targets. Results from individual predicted miRNA target searches of the TargetScan, DIANA-microT-CDS, and the miRDB were combined with the results of the validated miRNA target database miRTarBase. Putative target genes were assessed for their associated Gene Ontology (GO) terms using the online functional annotation tool DAVID. In addition, a NCBI gene and PubMed search was performed for target gene associations with cancer (metastasis) using the HGNC symbol.

High-depth multiplex-NGS of miRNAs from tumor cells of 50 patients in the discovery cohort, randomly divided into four sample pools, yielded variable raw library sizes with an average of >7 million reads. The total number of reads aligned to miRNA sequences is denoted the miRNA library size of a particular tumor specimen. Spreading of the amount of sequence reads (e.g., library sizes) across samples is shown in Supplementary File S1A. Analysis revealed expression of 2,707 from 7,160 annotated miRNA isoforms (miRBase v19) in at least 1 of 50 tumor specimens. After removal of low-expression miRNAs, 667 miRNAs remained eligible for statistical analyses. miR21-5p was most abundant in our dataset, representing an average of almost 10% of the miRNA sequence reads per sample. The spread of percentages of the total miRNA library sizes represented by the three most abundantly expressed miRNAs is shown as box plots in Supplementary File S1B, per miRNA.

Low-depth NGS with the Ion Torrent was used to sequence 13 of the 50 tumor specimens for technical validation. Sequence depths ranged from 50 to 1,500 reads per sample. Analysis revealed expression of 325 annotated miRNA isoforms in at least one of the 13 tumor samples. All these miRNAs were also detected in the high-depth HiSeq NGS analysis. Correlation analysis of paired samples in the HiSeq and Ion Torrent datasets confirmed reproducibility of miRNA expression across the different NGS platforms, even with considerable differences in sequence depth between the two platforms (Supplementary File S2A).

The expression of the 667 miRNAs detected by HiSeq NGS were assessed for differences between metastasized and nonmetastasized primary tumors. Differential expression analysis identified 32 miRNAs with nonadjusted P values of 0.02 or lower (Table 2, Supplementary File S3). From these miRNA isoforms, 19 (59%) showed lower expression and 13 (41%) showed higher expression in metastasized tumors compared with nonmetastasized tumors. This top list indicated that both the 3p and 5p variant of miR378 were differentially expressed. In addition, for miR1, miR135a, miR365, and miR550a two different transcripts were differentially expressed. Eight of the top ten listed miRNAs were selected for technical validation with PCR.

miRNA expression patterns were assessed to identify reference miRNAs for (technical) validation using qRT-PCR. Those miRNAs with the best linear relation between the number of isoform-specific miRNA reads (or miRNA expression) and the total number of miRNA reads (or miRNA library size) across the samples were identified by calculating the CV per miRNA. A low CV indicates low SD combined with high average expression. Calculated CVs and data of subsequent linear modeling of the two miRNAs with lowest CV, being miR16-1-5p and miR26a1-5p, are shown in Supplementary File S4A. The linear relation of miRNA expression and miRNA library size, and the residuals plots were plotted for graphical inspection (Supplementary File S4B). By selecting two miRNAs as combined reference, we aimed to produce a more robust reference, as compared with a single miRNA. For normalization of the amount of input material in the validation phase, we identified the combination of miR16-1-5p and miR26a1-5p as the most suitable reference. Validation of this combination of miRNAs was performed using online available miRNA sequencing data from the cancer genome atlas (TCGA) data portal (Supplementary File S4C). Results endorse our findings. In addition, normalization of our Hiseq NGS data with the combination of these two miRNAs resulted in similar results in differential expression analysis compared with the selected normalization of NGS data as described in the Materials in Methods section (data not shown). These results indicate the capability/suitability of this combination of miRNAs for normalization purposes of quantitative PCR of colon tissue–derived miRNA expression.

For technical validation in the discovery cohort and validation of results in the validation cohort, qRT-PCR was performed for 8 top ten listed miRNAs identified in the differential expression analysis of HiSeq NGS data (Table 2). Expression was calculated relative to the average expression of miR16-1-5p and miR26a1-5p. Tested miRNAs were technically validated with Pearson's r and Spearman's ρ correlation analysis (Supplementary File S2B). Selection of variables that can provide robust prediction is an important aspect for survival model building. miRNAs varying strongly between cohorts can by definition not serve as such robust predictors. Therefore, selection of miRNAs for survival model building was based on miRNAs showing both similar expression patterns between NGS and qRT-PCR (discovery cohort), and similar relative expression ranges in the discovery and validation cohort as assessed by qRT-PCR. Four of the assessed miRNAs, namely miR25-3p, miR339-5p, miR378-3p, and miR550a-1-5p, showed similar expression patterns between the NGS and qRT-PCR data in the discovery cohort. For miR25-3p and miR339-5p relative expression ranges were similar between the tumors of the discovery and the validation cohort. Differential expression results, as identified by HiSeq NGS, were supported by qRT-PCR results (Fig. 1). In both patient cohorts, a lower expression of either miRNA was associated with the occurrence of a metachronous metastasis.

Survival model building was performed on the qRT-PCR data of the discovery cohort using univariate and multivariate Cox regression analyses on miR25-3p and miR339-5p together with predetermined covariates. Afterwards, for comparison, we performed the survival analyses in the validation cohort. Except for survival status and distant metastasis status, neither individual expression of miR25-3p nor miR339-5p expression, or their combination was associated with any other clinicopathologic parameter (Supplementary File S5). Survival analyses were performed to assess differences in DRFP between low and high relative expression groups for the two individual miRNAs. Cut-off values for group classification were determined in the discovery cohort at 0.279 and 0.0143 for miR25-3p and miR339-5p, respectively. ROC curves for the individual miRNAs were used to assess the specificity and sensitivity of the chosen cut-off values in the discovery set, and confirmed them to be optimal (Supplementary File S6). Univariate analyses showed an association between high relative miRNA expression and shorter DRFP for either miRNA in both the discovery and the validation cohort (Table 3). Histograms showing the relation between height of miRNA expression and number of patients with a distant metastasis at the time point 5 years from surgery (Fig. 2A), and cumulative incidence curves for DRFP (Fig. 2B) were plotted. Multivariate analyses identified both miR25-3p and miR339-5p as independent prognostic factor for DRFP (Table 3). HRs were consistently more pronounced in the validation cohort compared with the discovery cohort.

Combination of biomarkers may better reflect tumor aggressiveness and might increase the clinical discriminative and prognostic value. Therefore, expression of miR25-3p and miR339-5p were combined into two groups, being high expression of both miRNAs and low expression of one or both miRNAs. Cumulative incidence curves illustrating univariate survival analyses for DRFP were plotted (Fig. 2C). The 5-year survival percentage for DRFP in the discovery cohort and the validation cohort was, respectively, 70% and 93% for patients with a high miRNA expression for both miRNAs compared with, respectively, 39% and 31% for patients with low expression for one or two of the miRNAs. Multivariate analyses identified combined expression levels of miR25-3p and miR339-5p as independent prognostic factor for DRFP (Table 3). Patients with a high expression of both miRNAs were shown to have an 11 times lower change of developing a distant metastasis in the validation cohort based on multivariate survival analyses.

Using the multivariate survival model combining miR25-3p and miR339-5p for DRFP estimated on the discovery cohort (Table 3), we constructed a nomogram to quantitatively predict the probability of metachronous distant metastasis in new cancer patients (Fig. 3A), at 3 years and 5 years after surgery. The nomogram was based on primary tumors of TNM stage II–III colorectal cancer patients with a microsatellite stable phenotype. The miRNA-based classifier and the predetermined most important nonsubjective clinicopathologic risk factors were integrated in the nomogram. The calibration plot showed that the prediction performance of the nomogram was good compared with the ideal model (Fig. 3B). The predictive accuracy of the nomogram was assessed in an independent patient cohort with a ROC analysis (Fig. 3C). The concordance index (Fig. 3C) for the multivariate survival model (Table 3), indicated good probability prediction for the prognostic multivariate survival model in the independent validation cohort.

miRNAs have numerous target genes and exert tissue-specific transcriptional regulation. To assess possible target genes for miR25-3p and miR339-5p, gene target prediction was performed (Supplementary File S7A) by combining results from three well-established prediction tools (TargetScan, DIANA-microT-CDS, and miRDB) and a database containing validated targets (miRTarBase). For miR25-3p and miR339-5p we found, respectively, 38 and 36 overlapping target genes predicted in all three target prediction databases, indicating poor agreement between the target prediction databases. GO term enrichment was performed for these overlapping targets (Supplementary File S7B). In addition, only 1 or 2 of the overlapping predicted target genes for miR25-3p or miR339-5p was validated according to the miRTarBase. The number of common target genes for miR25-3p and miR339-5p, predicted by the same database(s), was 66 (Supplementary File S7C). We assessed the associated GO terms of the putative target genes using online annotation tool DAVID, and performed a NCBI gene and PubMed search using the HGNC symbols (Supplementary File S7). Annotation analysis of putative target genes suggests involvement of these genes in cellular metabolic processes, indicating a possible role for the identified miRNAs in cellular metabolism. In conclusion, functional studies are needed to unravel the actual target genes of miR25-3p and miR339-5p in colon cancer cells.

In this study, we identified the combination of miR25-3p and miR339-5p expression levels in tumor cells as independent prognostic factor for the occurrence of distant metastasis in TNM stage II–III colon cancer with a microsatellite stable phenotype. We proposed, and retrospectively validated, in an independent patient cohort, a nomogram for prediction of distant metastasis for this group of patients. In addition, we recommend combination of miR16-5p and miR26a-5p as standard for normalization in qPCR of colon cancer tissue–derived miRNA expression.

Multiple different processes are involved in carcinogenesis, with many key genes/proteins. The hallmarks of cancer (35) represent properties that a cell needs to attain to become and sustain itself as cancer cell. These hallmarks direct the cellular pathways involved in carcinogenesis, and particularly metastasis formation. Using expression of miRNAs for prognostication of clinical outcome for cancer has advantage over mRNAs, as miRNAs are considered the key posttranscriptional regulators of gene expression, with numerous simultaneous target genes (15) and involvement in various cellular pathways. In contrast to mRNA, these posttranscriptional master regulators, are highly conserved between species (2).

Several miRNA profiles for prediction of metastasis risk in patients with colon cancer have been proposed before (20–24). Recently, a six miRNA-based classifier, differing from the combination of miRNAs reported here, was published (24). Zhang and colleagues identified a six miRNA-based classifier as independent prognostic factor for DRFP in TNM stage II colon cancer patients. However, the discriminative value of these miRNAs could not be confirmed with our NGS data, as individual miRNAs were neither differentially expressed between patients with and without metachronous metastasis nor prognostic for distant metastasis risk. Differences in identified miRNAs might be due to the substantial dissimilarities of the respective study designs and patient cohorts. Given the large heterogeneity of colon cancers, with regard to the amount of stroma and immune infiltrate, it is essential to use microdissected cancer cells and high-depth sequencing of the micronome. A combination of these approaches with the specific assessment of differential expression between patients presenting with or without metachronous metastasis allows prediction of an accurate prognostic (multivariate) model for metastasis risk. Moreover, it is crucial to subsequently validate the model in an independent patient cohort. Exclusion of patients with low percentages of cancerous epithelial cells leads to biased patient cohorts, as high stroma content of tumors is itself associated with worse clinical outcome (36). Inclusion of these patients requires selection of target cells by microdissection (18), due to the tissue-type and cell-specific nature of miRNA expression (16, 17). To obtain the most uniform cell population, we used LCM, which is a laborious technique most likely not applicable for routine clinical use. Alternatively, manual microdissection can be used, which is already standard procedure for assessment of PCR-based biomarkers in many diagnostics laboratories. In addition, introduction of automated microdissection using robotics is upcoming. To the best of our knowledge, this is the first study to combine LCM with high-depth sequencing of miRNAs, aiming to address differential expression between patients presenting with and without distant metachronous metastasis for accurate prognostication.

A relatively small number of miRNAs account for the majority of raw sequence reads in the miRNA libraries, as seen before (22). We hypothesized that these miRNAs, being miR21-5p, miR192-5p, and miR19b-2-3p, are important in colon cells, and probably in colon carcinogenesis, but not necessary in the process of metastasis. miRNA-21-5p was previously reported to show increased expression in tumor compared with normal colon mucosa (37), and was linked to early onset carcinogenesis (38). Two independent manuscripts reported that increased expression of miR21, predominantly in the stroma compartment and not the cancer cells, is associated with clinical outcome in TNM stage II colon cancer patients (39, 40). These findings might explain why we did not see an association between miR21 expression and metastasis risk, unlike others (24). Decreased expression of miR192 was identified in tumor compared with normal colon tissue (41) and linked to early onset colon carcinogenesis (38). Upregulation of miR19b was associated with colon cancer (42) and with 5-fluorouracil resistance in colon cancer cells (43).

On the basis of the linear relation between miRNA expression and miRNA library size, we identified combination of miR16-5p and miR26a-5p as most suitable reference for qPCR. Both miRNAs have been recognized before as suitable reference miRNAs (44) in a TaqMan microarray-based approach. While our analysis was limited to TNM stage II–III colon cancer patients, Chang and colleagues validated the miRNAs in a colorectal cancer cohort of TNM stage I–IV patients. These according results confirm the utility of these miRNAs as reference for qPCR in cancer tissues of colorectal origin.

As discussed previously (45), function of miR25 can be oncogenic or tumor repressive depending on the tissue type. Generally, expression of miR25 is higher in cancer cells of epithelial origin, and in tumor stroma, compared with normal colon mucosa (45, 46). We showed that expression of miR25-3p is lower in patients with colon cancer with distant metachronous metastases compared with those without distant metastases. In prostate cancer cell lines, a similar expression pattern was observed (47), with increase in expression of carcinoma cell lines compared with immortalized normal cells, and subsequently lower expression in cell lines from prostate carcinoma metastasized to distant locations compared with those of primary tumors (47). Taken together, these data suggest that in colon cancer, expression of miR25-3p, after initial increase, might need to be decreased to some extent to enable the initiation of the distant metastasis process. Indeed, Li and colleagues showed that overexpression of miR25 suppressed migratory ability of colon cancer cells in vitro and decreased tumorigenic potential in a colon cancer xenografts (45). In addition, Baffa and colleagues reported decreased expression of miR25 in metastases compared with their associated primary colon tumor (48). Expression of miR339-5p has been linked to the process of distant metastases in breast cancer (49). Similar to the association we observed for miR339-5p, higher expression was associated with lower metastasis rates, and in addition reduced cancer cell growth and better patient survival in breast cancer (49). Also for prostate cancer, lower expression of miR339-5p was observed in metastatic cancer (50).

To date, treatment allocation of patients with colon cancer is based on the TNM classification system. Under current classification, patient groups still show large differences in clinical outcome (9–11). Moreover, clinically defined risk criteria for adjuvant therapy in locally restricted colon cancer are currently considered to be of insufficient quality (9, 12). The development of distant recurrent disease mainly accounts for colon cancer–related mortality. The use of the miRNA combination reported here, capable of predicting metastasis risk in colon cancer, in conjunction with the TNM staging method, may contribute to an approach of personalized treatment based on individual tumor characteristics. In conclusion, in this bicentric international study, we retrospectively identified and validated a miRNA classifier in the primary tumor using two independent patient sets, and propose a nomogram capable of predicting metastasis risk in microsatellite stable TNM stage II–III colon cancer.

No potential conflicts of interest were disclosed.

Conception and design: I.J. Goossens-Beumer, C.J.H. van de Velde, P.J.K. Kuppen

Development of methodology: I.J. Goossens-Beumer, R.S. Derr, H.P.J. Buermans, P.J.K. Kuppen

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): I.J. Goossens-Beumer, R.S. Derr, H.P.J. Buermans, H. Morreau, U. Nitsche, K.-P. Janssen, C.J.H. van de Velde, P.J.K. Kuppen

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): I.J. Goossens-Beumer, H.P.J. Buermans, J. Goeman, S. Böhringer, U. Nitsche

Writing, review, and/or revision of the manuscript: I.J. Goossens-Beumer, R.S. Derr, H.P.J. Buermans, S. Böhringer, H. Morreau, U. Nitsche, K.-P. Janssen, C.J.H. van de Velde, P.J.K. Kuppen

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): I.J. Goossens-Beumer, H. Morreau, P.J.K. Kuppen

Study supervision: C.J.H. van de Velde, P.J.K. Kuppen

The authors thank Sabine Leis and Anne Benard for their help in the collection and processing of frozen tissue samples of the validation cohort.

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.