Hsa-miR-31-3p Expression Is Linked to Progression-free Survival in Patients with KRAS Wild-type Metastatic Colorectal Cancer Treated with Anti-EGFR Therapy

With more than 1.2 million new cases in 2008, colorectal cancer is the third most common form of cancer worldwide (1, 2). Nearly 20% of patients are diagnosed at an advanced stage with metastatic disease (1) and over half will ultimately develop metachronous metastases. The recent development of targeted therapies has improved outcomes in patients with advanced colorectal cancer. Cetuximab and panitumumab, two monoclonal antibodies, which neutralize the extracellular domain of EGFR, have demonstrated effectiveness in terms of both response and survival when used as first, second or third-line chemotherapy. Although anti-EGFR agents are routinely used in clinical practice, their use is limited to 60% of patients who have KRAS wild-type tumors (2) because randomized-controlled studies have shown activating mutations of this oncogene were predictive of resistance to these agents (3–12).

Because response to anti-EGFR therapy is seen in less than 40% of patients with KRAS wild-type tumors (13), additional factors are needed to help patient selection for this therapy to avoid the prescription of ineffective, expensive, and potentially harmful treatments. Many studies have evaluated other potential biomarkers, typically using a candidate gene approach focused on different intracellular pathways downstream of EGFR or a genome-wide approach. At the present time, however, the routine use of these biomarkers is not warranted. Indeed, BRAF V600E mutation seems to be a marker of poor prognosis (9, 10, 14, 15). Furthermore, studies that examined PIK3CA mutations, loss of PTEN expression, EGFR gene copy number, and polymorphisms of fragment c γ receptors have reported conflicting results (14, 16–27). Finally, mRNA expression level of amphiregulin and epiregulin seems to be of limited clinical interest because of the inability to set an appropriate threshold limit (28, 29).

MicroRNAs (miRNA) are small 19 to 25 nucleotides, noncoding RNAs that negatively regulate 30% of gene expression posttranscriptionally by inhibiting translation and degrading mRNAs. miRNAs control biologic processes such as cell proliferation, differentiation, angiogenesis, and apoptosis (30), while also acting as oncomiR or suppressor miRs depending of their target mRNA. miRNAs deregulation, through expression modification and point mutations, is involved in the occurrence of many types of cancer, including colorectal cancer (31). These abnormalities are of major interest and could serve for cancer diagnosis, prognosis, and chemotherapy response prediction.

This study aimed to develop a miRNA expression–based model for predicting survival of patients with KRAS wild-type metastatic colorectal cancer (mCRC) treated with anti-EGFR monoclonal antibodies.

Four cohorts of patients with mCRC were included in this study (n = 132). Three were used as training sets and one as a validation set. All samples were obtained from studies that had previously received appropriate ethical committee approvals. Samples for the first training series (TS-Frozen1) were obtained from a retrospective collection of 33 patients, some of whom were reported on in a previously published series (13, 20). All patients were refractory to a FOLFOX or FOLFIRI regimen, refractory to these regimens plus anti-EGFR antibodies, or received panitumumab as a single agent.

Samples for the second training series (TS-Frozen2) were obtained from a prospective collection of patients treated by anti-EGFR antibodies alone or in combination with an irinotecan-based chemotherapy regimen. All of these patients were considered refractory to FOLFIRI regimen.

The third training set (TS-FFPE) was made up of samples from a retrospective series of pooled Belgian and Finnish patients with mCRC. Some of these patients were a part of previously published study (32). All patients in this training set were considered refractory to a 5-fluorouracil–based regimen combined with irinotecan and oxaliplatin.

The last cohort of patients (validation set, VS) was comprised of patients from a randomized phase II trial sponsored by GERCOR (Groupe Cooperateur Multidisciplinaire en Oncologie). From the 65 patients eligible for this phase II study, 50 KRAS wild-type patients were considered for subsequent analysis. Among them, 45 patients had enough tumor material to be included (84%) with 26 having only formalin-fixed paraffin-embedded (FFPE) samples available, 6 with fresh-frozen (FF) samples only, and 13 having both FFPE and FF samples available. All patients in this cohort were treated with third-line therapy using a combination of irinotecan and panitumumab after a demonstration of progression to oxaliplatin- and irinotecan-based chemotherapy regimens (33). This validation cohort was divided in two (Fig. 1A and B). The first cohort (VS-Frozen) consisted of patients with FF samples (n = 19) with the second cohort (VS-FFPE) consisting of patients with FFPE samples without FF samples available (n = 26) or with FF samples (n = 39).

DNA and RNA extraction and mutation analysis was only performed on specimens with greater than 50% tumor component. DNAs were extracted from frozen matched tumor and nontumor tissue samples using the QIAamp DNA Mini Kit (Qiagen). Total RNAs were extracted from frozen tumors, cells or FFPE using the mirVanamiRNA Isolation Kit (Ambion), miRNeasy, and miRNeasy FFPE Extraction Kit (Qiagen), respectively (see Supplementary Methods).

Methods of KRAS and BRAF mutation detection are described in the Supplementary Methods (13, 34).

Global miRNA profiling was performed by labeling and hybridizing 750 ng of extracted RNA from 43 tumor frozen samples randomly chosen from TS-Frozen1 and TS-Frozen2 on Illumina Human v2 microRNA Expression Beadchips. Beadchips were scanned with the Illumina I-Scan Reader and data were imported into GenomeStudio data analysis software (Illumina), quantile-normalized and log₂-transformed. For expression analysis on TaqMan (Applied Biosystems) of the 11 selected miRNAs, see Supplementary Methods.

Three colorectal adenocarcinoma cell lines (HTB-37, CCL-222, and CCL220-1) from the American Type Culture Collection (ATCC) were selected and used in the 3 month following reception of the cells for cell culture studies. The authentication of the cells has been provided by ATCC using the short tandem repeat profiling method. All cells were transfected with miRVanamiRNA (Ambion) mimic negative control or hsa-miR-31-3p miRVanamiRNA mimic (see Supplementary Methods).

We developed a database integrating contemporary miRNA target predictions from six individual databases (PITA; picTar 5-way Targetscan microRNA.org; MicroCosmandmiRDB). These databases allowed for the determination of miRNAs cotargeted by candidate genes, taking into account the number of miRNA prediction databases that predicted each miRNA/target relationship and the rank of this prediction. This database was current as of November 2012, when our analysis was performed.

Prediction of survival risks used the miRNA expression data in a supervised principal component (SuperPC) approach analysis, a modification of the conventional perchloric acid as described by Bair and colleagues (35, 36). To evaluate the predictive value of the method, a leave-one-out cross-validation was performed and permutation analysis (n = 100) estimated the statistical significance of separation of the Kaplan–Meier survival curves. A prognostic index is provided for each patient with a log expression profile given a vector x; a high value of the prognostic index corresponding to a high value of hazard of death.

Survival statistical analysis was performed using the R packages “survival” and “rms.” Univariate and multivariate analyses used a Cox proportional regression hazard model generating an HR and a performance index (area under the curve, AUC). Nomograms were developed on the basis of Cox proportional regression hazard models, which predict the probability of progression-free survival (PFS) based on a training set of patients. Internal validation relied on 150 bootstrap resamples, then independent validation. The validation set series and whole series (training + validation series) were submitted to a multivariate Cox model analysis, including the Eastern Cooperative Oncology Group (ECOG) performance status and the level of lactate dehydrogenase (LDH) when available. These two variables were obtained at the time of the first administration of anti-EGFR therapy. The analysis was performed by STATA version 11. The patients were dichotomized in two groups, high and low risk, as determined above. The ECOG performance status was obtained in 108 patients. The LDH have been prospectively recorded for the patients included in the phase II PIMABI study.

Gene and miRNA expression value comparison analyses were performed using the nonparametric test (Kruskal–Wallis tests) with the pairwise Wilcox test function in R using the stats library.

False-discovery rate–adjusted P values were calculated using the Benjamini and Hochberg procedure for multiple testing correction. The correction test function was used to calculate Pearson correlations between expression values together with matching P values. Statistical significance was set at P < 0.05 for all analyses.

Patient characteristics of the different sets are described in Table 1. Overall, 132 patients with mCRC were included in this study. All patients were wild-type KRAS based on screening of exons 2, 3, and 4. From training sets 1 and 2, 43 RNAs were extracted from frozen tumor samples and analyzed by Illumina Beadchips. Using SuperPC analysis, univariate association between each miRNA and PFS allowed us to select a set of 11 miRNAs, for which Cox score statistics were found significant (P < 0.01; Supplementary Table S2 and Supplementary Fig. S1).

All miRNAs were tested for in the first and second training sets using TaqMan probes. Only one miRNA, hsa-miR-31-3p, showed a significant association with PFS on TS-Frozen1 and TS-Frozen2 with HRs being 1.8; 95% confidence interval (CI), 1.1–2.9; P = 0.01 and 2.5; 95% CI, 1.3–4.5; P = 0.002, respectively. In addition, a Super-PC analysis combining Cox proportional hazard model and principal component analysis performed on TS-Frozen1 using the TaqMan expression level of miR-31–3p allowed us to define a prognostic score, computed by the formula: 0.178x − 1.36, where x is the TaqMan log-transformed expression of hsa-miR-31-3p, and a cutoff value of −0.031 dividing patients into two groups with a high and low risk of progression (Supplementary Fig. S2). The HR was 2.6; 95% CI, 1.15–5.8; P = 0.021; the PFS of high- and low-risk patients was 13 versus 31.4 weeks, respectively. This prognostic score was calculated and the cutoff value was applied on TS-Frozen2. The PFS of high- and low-risk patients was 9 versus 35.3 weeks, respectively (HR, 4.1; 95% CI, 1.3–13.2; P = 0.018; Fig. 2A). The model was further validated for patients with FF sample in the PIMABI phase II study (validation set frozen, n = 19). The PFS of high- and low-risk patients was 23 versus 48.8 weeks, respectively (HR, 4.1; 95% CI, 1.1–15.3; P = 0.02; Fig. 2B).

To be independent of a dichotomized variable, we built a nomogram based on Cox proportional hazards regression modeling patients PFS probability with the log-transformed hsa-miR-31-3p expression on the TS-Frozen1 (Fig. 3). The BRAF status, a known prognostic marker in mCRC, was added into the model. The nomogram was then tested on TS-Frozen2 and VS-Frozen. The nomogram demonstrated the ability to predict PFS with an AUC superior to 0.70, whatever the time threshold tested (Supplementary Fig. S3).

In anticipation of the use of this marker in clinical practice, we evaluated its performance on RNA extracted from FFPE tissues. To validate the prognostic value of the expression of hsa-miR-31-3p on FFPE samples, we used two sets of patients in which FFPE tissue samples were available (TS-FFPE, n = 35 and VS-FFPE, n = 26). Patients from TS-FFPE were divided into two groups according to expression of hsa-miR-31-3p, as was done for frozen samples using a Cox model associated with a principal component analysis to calculate a new prognostic score (0.093x − 0.407, where x is the FFPE TaqMan log-transformed expression of hsa-miR-31-3p) and a new cutoff value (−0.058). An HR for PFS of 2.79; 95% CI, 1.12–6.96; P = 0.027 was observed between the high- and low-risk group of patients. This threshold was then validated on VS-FFPE from the PIMABI trial for patients in which only FFPE samples were available. An HR for PFS of 5.1; 95% CI, 1.48–14.44; P = 0.01 was observed.

We then added patients from VS-Frozen for whom we also had FFPE samples (n = 13). The PFS HR between high- and low-risk patients was of 2.44; 95% CI, 1.1–5.4; P = 0.028. When we added 6 patients to this cohort for whom only FF tissues were available to test whole series of PIMABI, the PFS HR between high- and low-risk patients was of 2.46; 95% CI, 1.22–4.94; P = 0.012 (Fig. 4). Finally, we observed a significant inverse correlation between the percentage variation of the RECIST criteria and hsa-miR-31–3p expression (r² = 0.49; P = 0.0035). This result was confirmed by the significant association between patient risk status and the percentage of variation of RECIST criteria (P = 0.02; Kruskal–Wallis rank test, Fig. 3). Multivariate analyses were performed, including the BRAF status (i.e., mutated vs. nonmutated), the LDH level and the ECOG performance status at time of inclusion in PIMABI phase II study in RAS wild-type population. The status of patients determined by the level of expression of the hsa-miR-31-3p remains significant. The PFS HR and overall survival (OS) HR between high- and low-risk patients were 2.6; 95% CI, 1.1–5.8; P = 0.023 and 3.2; 95% CI, 1.4–7.5; P = 0.008.

Finally, when we pooled the different series (n = 132) taking into account the classification based on the expression of hsa-miR-31-3p (Supplementary Fig. S4). The OR of nonresponse in the high-risk patients was 4.9; 95% CI, 2–12.5; P < 0.001). We performed a multivariate analysis, including the ECOG performance status and the BRAF status in RAS wild-type patients (n = 104). The PFS HR and OS HR between high- and low-risk patients were 2.3; 95% CI, 1.4–3.8; P < 0.001 and 1.9; 95% CI, 1.2–3.1; P = 0.008.

Three colorectal cancer cell lines were transfected with hsa-miR-31-3p mimic or with a mimic control. The transfection efficacy was demonstrated by a 1,500 times average increase of hsa-miR-31-3p levels. Expression profile analysis of transfected cells allowed us to identify 47 genes significantly downregulated (fold change < 0.7; P < 0.05), and 27 genes significantly upregulated by hsa-miR-31-3p (fold change < 1.3; P < 0.05; Supplementary Table S3).

As the role of a miRNA includes degradation of its transcript target, we studied whether our database predicted the 47 downregulated genes as a hsa-miR-31-3p putative target. Twenty-five of the genes were predicted to be putative direct targets of hsa-miR-31-3p and displayed a good rank in the prediction database. This number and the ranking of genes were significant (P < 0.0001 for both by the permutation test).Twenty-six of the upregulated genes were not predicted to be direct targets of hsa-miR-31-3p. Although one of the upregulated genes was predicted to be a direct target, it was the last ranked target in the prediction database.

The 25 putative direct and 27 indirect target genes were validated on qRT-PCR. Out of these 52 genes, 45 displayed an expression level comparable with the level obtained in the array. When the expression of these genes was analyzed in 39 tumor samples, we found a correlation between hsa-miR-31-3p expression pattern and PFS for DBNDD2 (P = 0.02) and EPB41L4B (P = 0.009). Interestingly, both genes displayed a negative correlation with hsa-miR-31-3p expression levels: DBNDD2 (−0.5; P = 0.001) and EPB41L4B (−0.3; P = 0.04; Supplementary Fig. S5).

By using an assumption-free approach to identify miRNAs associated with PFS in patients with KRAS wild-type mCRC treated by anti-EGFR therapy, we established for the first time a link between hsa-miR-31-3p expression and the risk of progression. All patients in our cohorts were refractory to the chemotherapy regimen used with anti-EGFR therapy. We were able to define a robust threshold by dividing patients in two groups according to their risk of progression after anti-EGFR therapy. Our results independently confirm those of a Finnish group (32) who demonstrated a significant differential expression of hsa-miR-31-3p between patients with stable and progressive disease, even though we used a part of their patient set in TS-FFPE to determine the threshold of the hsa-miR-31-3p in FFPE tissues.

The availability of multiple drugs for mCRC (e.g., panitumumab, cetuximab, and bevacizumab) underlines the need for better decision-making tools. The survival improvement of patients with mCRC often results from the use of these different drugs via a personalized approach or precise manner. The nomogram we developed, which predicts PFS on an individualized basis, provides such a decision-making tool. The accuracy of this nomogram on frozen samples was assessed by the calculation of the AUC that gives a value superior to 0.70, which can be considered as a good performance for a predictor. There is also a good correlation between hsa-miR-31-3p expression levels measured in frozen samples and those measured in FFPE samples, confirming the interest to use miRNAs as surrogate for wider gene expression signatures. Nevertheless, we needed to readjust the threshold of the hsa-miR-31-3p expression to divide the high- and low-risk progression groups for FFPE samples. Furthermore, although PFS can be considered a good surrogate marker of efficacy in well-designed randomized clinical trials, its value per se in retrospective series is unclear. We also investigated the impact of hsa-miR-31-3p with other critical variables such as response rate and OS. There was a significant link between response status (nonresponder vs. responder) and the classification in high- and low-risk patients according hsa-miR-31-3p expression as shown in the waterfall plot of the patients included in the PIMABI phase II study. We also validated this association in the whole series, the risk of nonresponse is five times more frequent in the high-risk patient group than in the low-risk patient group. Finally, we also showed a significant link between the OS and the hsa-miR-31-3p expression level in the whole series adjusted on the other prognostic variables (ECOG and BRAF mutational status). The high-risk patient group has a significant shorter survival than the patients belonging to the low-risk group.

Several miRNAs have been shown to be associated with response to chemotherapy in different types of cancers (37–39). To our knowledge, only one group has reported an association between the expression of a miRNA (hsa-miR-143) with PFS and cancer-specific survival in patients with mCRC treated by anti-EGFR alone or in combination with chemotherapy in refractory setting (40). The authors, however, failed to identify a link between the expression of this miRNA and response to anti-EGFR therapy, suggesting that this biomarker is more a prognostic marker than predictive, even if some relation with 5FU sensitivity have been described (41)

The specific role of the hsa-miR-31-3p has been investigated through in silico and in vitro approaches and the list of genes is not easy to interpret. The 25 genes putatively directly downregulated by the hsa-miR-31-3p do not give a clear picture to a regulatory mechanism of response to anti-EGFR therapy. A recent link between some specific regulation of miRNAs maturation and EGFR kinase activity has been emphasized and could enlighten our results (42). In hypoxic conditions the maturation of pre-miRNAs with long loop structures are dependent on the phosphorylation of AGO2 protein by EGFR. The phosphorylated AGO2 protein reduced its interaction with DICER protein leading to decrease significantly the loading of pre-miRNA or the expression of mature miRNA. Among the miRNAs regulated by this mechanism, the premature mir-31 is one of the most likely candidates. The absence of EGFR induced increased expression of mir-31 and inhibition of EGFR by tyrosine kinase inhibitors such as gefitinib. The high level of hsa-miR-31-3p expression in our patient's tumors could be the result of an absence of EGFR response to hypoxia (43) and thus reduced benefit from anti-EGFR therapy. Among the genes derepressed by the inhibition of maturation of miRNA by EGFR response to hypoxia (42), two genes were found in common with our list of downregulated genes by the hsa-miR-31-3p; AMFR and EPB41L4B, the latter being also associated with PFS.

We were unable to demonstrate the predictive role of hsa-miR-31-3p expression in our study because there was no control arm consisting of patients without anti-EGFR therapy. In an attempt to assess the predictive nature of hsa-miR-31-3p, we first tested its prognostic value in the PETACC3 cohort (44, 45) and did not find a correlation between hsa-miR-31-3p expression and survival after relapse (HR, 0.93; P = 0.48; data not shown). These data were based on 124 KRAS and NRAS wild-type patients. Furthermore, in a recent study by Shen and colleagues (42), the expression of phosphorylated AGO2 protein has been shown to be associated with worse prognosis in patients with breast cancer, suggesting that maturation inhibition of miRNAs controlled by EGFR is a marker of aggressiveness in at least breast cancer patients. Our results are exactly at the opposite, suggesting that high expression of hsa-miR-31-3p is rather a predictive than a prognostic marker.

Altogether, our results suggest that the expression hsa-miR-31-3p could be of help in the decision to add anti-EGFR therapy to common chemotherapy regimens in the population of KRAS wild-type patients by giving a probability of PFS, which will be of interest with the development of several targeted therapy in mCRC.

R. Thiébaut, F. Liébaert, F. Rousseau, and B. Génin are employees of IntegraGen. J.-B. Bachet is a consultant/advisory board member for Amgen and Merck Serono. O. Bouché reports receiving speakers bureau honoraria from Amgen and is a consultant/advisory board member for Merck Serono. F. Bibeau reports receiving speakers bureau honoraria from and is a consultant/advisory board member for Amgen, Merck Serono, and Sanofi. S. Tejpar reports receiving speakers bureau honoraria from and is a consultant/advisory board member for Merck Serono. P. Österlund reports receiving a commercial research grant from Roche, speakers bureau honoraria from Prime Oncology, and is a consultant/advisory board member for Amgen, Bayer, Merck Serono, Roche, and Sanofi. T. André reports receiving speakers bureau honoraria from Amgen and is a consultant/advisory board member for Amgen and Merck Serono. P. Laurent-Puig is a consultant/advisory board member for Amgen, Merck Serono, and IntegraGen. No potential conflicts of interest were disclosed by the other authors.

Conception and design: G. Manceau, R. Thiébaut, F. Liébaert, F. Rousseau, B. Génin, P. Österlund, S. Knuutila, P. Laurent-Puig

Development of methodology: S. Imbeaud, R. Thiébaut, F. Rousseau, B. Génin, D. Le Corre, A. Didelot, S. Knuutila, P. Laurent-Puig

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): R. Thiébaut, K. Fontaine, F. Rousseau, J.-B. Bachet, B. Chibaudel, O. Bouché, B. Landi, F. Bibeau, K. Leroy, F. Penault-Llorca, J.-L. Van Laethem, S. Tejpar, N. Mosakhani, P. Österlund, R. Ristamäki, V. Sarhadi, S. Knuutila, V. Boige, T. André, P. Laurent-Puig

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): G. Manceau, S. Imbeaud, R. Thiébaut, B. Génin, M. Vincent, F. Penault-Llorca, J.-L. Van Laethem, S. Tejpar, S. Rossi, N. Mosakhani, P. Österlund, S. Knuutila, T. André, P. Laurent-Puig

Writing, review, and/or revision of the manuscript: G. Manceau, R. Thiébaut, F. Liébaert, B. Génin, M. Vincent, J.-B. Bachet, B. Chibaudel, O. Bouché, B. Landi, J.-L. Van Laethem, P. Demetter, S. Tejpar, S. Rossi, P. Österlund, S. Knuutila, V. Boige, T. André, P. Laurent-Puig

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): R. Thiébaut, J.-B. Bachet, B. Landi, P. Demetter, P. Laurent-Puig

Study supervision: G. Manceau, R. Thiébaut, F. Liébaert, T. André, P. Laurent-Puig

The authors thank Céline Vazart and Virginie Decaulne for their expertise and contribution associated with the extraction and quantification of miRNA for this study.

This work was funded by the Programme Hospitalier de Recherche Clinique (PHRC), Institut National de la Santé et de la Recherché Médicale (INSERM) through the Biointelligence program, and by IntegraGen through the Biomos program funded by Banque Publique d'Investissement, France. PETACC3 research on miRNAs is supported by the Foundation Fournier-Majoie, FFMI, Belgium. S. Tejpar is a senior clinical investigator of the Fund for Scientific Research Flanders and is supported by Belgian National Cancer Plan.

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.