Expression of two parental imprinted miRNAs improves the risk stratification of neuroblastoma patients

Age at diagnosis, stage, and MYCN amplification are the cornerstones of the risk-stratification score of neuroblastoma that enables defining patients at low- and high risk. Refinement of this stratification is needed to optimize standard treatment and to plan future clinical trials. We investigated whether two parental imprinted miRNAs (miR-487b and miR-516a-5p) may lead to a risk score with a better discrimination. Expression levels of maternal miR-487b and paternal miR-516a-5p were determined using quantitative RT-PCR both for 231 neuroblastoma tumors (derivation set) and 101 independent neuroblastoma tumors (validation set). Survival outcomes were overall survival (OS) and disease-free survival (DFS). Multivariable Cox models were developed from derivation set and their performance evaluated using Akaike's information criterion (AIC) (goodness-of-fit) and time-dependent area under curves (discrimination). The selected model was validated using internal and external validation. The prognostic model including current prognostic factors plus miR-487b, miR-516a-5p, and interaction between two miRNAs was selected. Performance of this model was better in terms of both predictive ability (smallest AIC) and discrimination power (AUC close to 0.70). This model identifies three risk groups: high (3), intermediate (2), and low (1). Hazard ratios (HR) across risk groups were HR2/1 = 6.3 (2.7–14.6), HR3/1 = 14.8 (7.2–30.2) for OS and HR2/1 = 2.8 (1.5–5.4), HR3/1 = 7.2 (3.9–13.4) for DFS. The rank between these three risk groups was maintained and validated when performing internal and external validation. Expression of maternal miR-487b and paternal miR-516a-5p improves the risk stratification. This better discrimination at diagnosis is of clinical utility both for current and future treatments of neuroblastoma patients.


Introduction
Neuroblastoma, a malignant solid tumor of neural crest origin, accounts for 10% of all childhood cancers and exhibits a high clinical heterogeneity [1]. The combina-tion of age at diagnosis, tumor stage, and MYCN amplification status constitutes the score of the standard algorithm at diagnosis that enables to stratify patients into low-and high-risk groups, for both overall survival (OS) and disease-free survival (DFS) [2,3]. The low-risk group (up to 90% survival at 5 years) consists of non-MYCNamplified tumors that are either localized (stages 1, 2, and 3) or metastatic in children <18 months old (stages 4 and 4S). The high-risk group (around 35% survival at 5 years) comprises all cases MYCN-amplified neuroblastoma regardless of stage and age, plus non-MYCN-amplified stage 4 tumors in children ≥18 months old. Additional prognostic markers of relapse for low-risk patients and of survival for high-risk patients are required to optimize the management of neuroblastoma and to improve outcome for the patients [1,4]. Several attempts had been made to sort out molecular tumor markers of prognostic significance based on nucleic acids machinery (DNA, RNA) or biochemical traits [5,6]. So far the pan-genomic algorithm (gain or loss of whole chromosomes and intrachromosomal alterations) has been shown to improve the standard risk stratification [7,8]. On the other hand, several reports have identified gene expression patterns to predict the outcome of neuroblastoma patients [2,9]. Recently, Valentijn and colleagues defined a 157-genes signature that predicts clinical outcome of neuroblastoma irrespective of MYCN amplification [10].
MicroRNAs, due to their current exhaustive identification and chemical stability, represent another option to identify putative prognostic markers. These small RNAs interfere with diverse biological functions through the post-transcriptional regulation of gene expression by acting on the stability or translational rate of mRNA [11]. The involvement of miRNAs in tumor initiation and progression is well established and offers new perspectives in the therapeutic management of neuroblastoma patients [12,13]. Some recent works have been published in neuroblastoma patients. Each study proposed different miRNA signatures associated to neuroblastoma prognosis but not compared to the current classification systems [14][15][16][17]. De Preter and colleagues have proposed a 25-miRNA signature able to identify a cohort of highrisk neuroblastoma patients at greater risk of poor outcome [18].
We previously showed the prognostic role of miR-487b beyond the current risk factors in neuroblastoma tumors [19]. This miRNA belongs to the largest imprinted miRNA cluster (C14MC) at the 14q32.31 locus. Besides neuroblastoma tumors, alterations of expression of C14MC are frequently found in many other human cancers and diseases [20,21]. In human, the second large imprinted miRNA cluster, called chromosome 19 microRNA cluster (C19MC), locates at the 19q13.4 locus [22,23]. Several studies showed overexpression of C19MC miRNAs as associated with poor prognosis of many tumors [24,25].
In this study, we extended the work of Gattolliat and coworkers [19] in evaluating the prognostic of the expres-sion of only two miRNAs, each representing the most significant prognostic marker of C14MC and C19MC, the two imprinted clusters. To this purpose, we evaluated the added prognostic information from miR-487b and miR-516a-5p to the current prognostic factors (age, stage, MYCN) on the survival outcomes.

Materials and Methods
Tumor's collection, cryosections, and total RNA extraction Tumors samples of a retrospective cohort, partially previously published by Gattolliat and colleagues [19], were gathered, between 1987 and 2009, at the Gustave Roussy, Villejuif, France (n = 223) and at the Centre L eon B erard, Lyon (n = 8). This cohort was referred to as derivation set. An independent cohort of 101 patients from Centre L eon B erard followed between 1988 and 2010 was constructed for the external validation (validation set). All tumors samples were collected with the approval of the appropriate ethic committees, according to the national law on the protection of people taking part in biomedical research. Patients were staged according to the International Neuroblastoma Staging System [26]. Storage of primary tumor tissues, cryosections, and total RNA isolation were previously described [19]. The first and last cryosections were used to select tumor tissues with a malignant tumor cell content of ≥60% [1,27].

miRNA quantitative RT-PCR
The quantitative real-time polymerase chain reaction (qRT-PCR) and statistical analyses were done as previously described [19]. Briefly, qRT-PCR was performed according to the MIQE guidelines [28] and the quantification method was done according to Livak and Schmittgen [29]. RNU-44 was used as normalizer and IGR-N-835 cell line as calibrator [30]. The statistical analyses compared the miRNA expression level between low-and high risk using the twotailed Student's t-test. These analyses included a reanalysis of miR-487b previously identified [19] as well as an analysis of miR-516a-5p, the most significant miRNA marker from the C19MC cluster.

Statistical analysis
The median follow-up was calculated using the Schemper's method [31]. OS and DFS were defined as the time from the diagnosis to death and to the first occurrence of relapse or death, respectively. Survival curves were estimated using the Kaplan-Meier method and were compared with the log-rank test. Five multivariable Cox proportional hazard models were developed to assess individually (or not) the added prognostic information of two markers: one from C14MC and the second from C19MC. Model 1 included the three current prognostic factors and was later referred to as the reference model. Models 2 and 3 were built in adding one miRNA to the reference model. Model 4 was the reference model plus the two miRNAs. This model was extended including an interaction term between the two markers (model 5). All models were compared using the likelihood test ratio and Akaike's information criterion (AIC). The discriminant performance of models was measured using the timedependent area under the curve (AUC) using the Surviv-alROC R package [32]. This criterion allows evaluating whether a model may correctly classify a patient with an event and a patient without event at a time point; AUC close to 1 indicating a better classification. The performance of the final model was also evaluated in deriving risk groups from the prognostic index (PI) (sum of weighted regression coefficients estimated from derivation set). The discrimination of the observed survival curves of risk groups (Kaplan-Meier) and the comparison with final model-based predicted mean survival curves were studied. The final model was evaluated through internal and external validation. The first validation used a leaveone-out cross-validation that allowed to validate the development process of the prognostic model [33]. The external validation enabled evaluating the generalization of the results. Whatever the type of validation, we estimated the PI in relying on (1) reuse of derivation set for internal validation (PICV) and (2) the validation set for external validation (PIV). In this aim, using the regression coefficients estimated (1) from the estimation set for the internal validation and (2) from the derivation set for the external validation, we derived three risk groups in categorizing the continuous PI estimated previously using cut-points based on the distribution of (1) PICV for the internal validation and (2) PI across the individual in the derivation set. The choice of three groups was arbitrary, but when two were chosen (low-and high-risk groups), this was not sufficient to the routine clinical practice; and, choosing four, may be too large regarding the number of events [33]. We estimated the observed survival curves (Kaplan-Meier) into the three risk groups and evaluated whether they were well separated and reported the hazard ratios across these risk groups. The comparison with the predicted survival curves allowed evaluating the accuracy prediction of the model (calibration). Less importance was attached to this point since our aim was to define a risk score with good ability to separate patients with different prognoses. An update of the final prognostic model was performed on the pool of the derivation and validation set with adjustment by dataset, thus providing more accurate estimate of regression coefficients. For all analyses, we categorized continuous miRNA expression levels as measured by qRT-PCR into two classes (<1, ≥1), using as cutoff the expression level of IGR-N-835 cell line taken as one [19]. This cutoff was found to be very close to those computed for the miRNA expression levels that best discriminate patients for survival (overall and DFS) of patients (data not shown). All tests were two-sided with a nominal significance level of 5% and statistical analyses were carried out with SAS â version 9.3 software (SAS Institute, Inc., Cary, NC).

Description of the derivation and validation sets
Main characteristics of the two tumor sets used in our study are reported in Table 1, and individual patient data are presented in Table S1. The derivation set of 231 neuroblastoma patients was followed until March 2012. The median follow-up was 7.7 years (95% confidence interval [CI] = [6. 8-8.7]). According to the standard algorithm of disease outcome based on age at diagnosis, stage, and MYCN amplification, 65% (n = 151) patients were classified at low risk and 35% (n = 80) at high risk in the derivation set. The observed number of deaths was 69 (30%) and the number of events (death or relapse) was 90 (39%). In the external validation set (n = 101), the median follow-up of 6.7 years (5.4-8.0) was of the same order of magnitude of the derivation set. If patients were younger (P = 0.05) with higher number of MYCN amplification (P = 0.04) in the validation set, no significant difference into the standard risk stratification between the two datasets was observed (P = 0.38). There is no significant difference between derivation and validation set in terms of OS (P = 0.06) contrary to DFS (P = 0.04) (Fig. S1 reporting the log-rank test). These last two comparisons were performed after controlling for stage, age, and MYCN amplification but considering the analysis of retrospective data some important omitted confounders; for example, the difference of treatments in the institutions where these archival cohorts were obtained, could explain the observed difference in term of DFS. Figure 1A shows the updated survival curves reported in Gattolliat and coworkers with a significant log-rank test (P < 0.0001 for both OS and DFS). Indeed, within the cluster C14MC, miR-487b was the miRNA whose the expression was the most associated with both OS (P = 0.0018) and DFS (P = 0.0003) (Figs. 2 and 3, Table S2). The update of prognostic effect of miR-487b on OS and DFS reported in Gattolliat and colleagues is presented in column "model 2" of Table 2.
Strong expression of miR-516a-5p from the C19MC is associated with high-risk neuroblastoma Previous miRNA microarray profiling of a screening cohort of 13 tumors (eight high-risk and five low-risk) revealed that miRNAs from C19MC were globally overexpressed in high-risk neuroblastomas [19], in agreement with observations from other high-grade tumors of various cancers [25]. Among the microRNA of the C19MC cluster, miR-516a-5p was the only microRNA whose expression was significantly associated with survival. These results, confirmed by qRT-PCR in the derivation set (n = 231), showed that miR-516a-5p expression levels were higher in high-risk tumors (P = 0.0091, Fig. S4). Kaplan-Meier analyses showed that patients with high expression levels of miR-516a-5p showed a shorter OS (P = 0.0058) but not a shorter DFS (P = 0.1087), compared to patients with low expression levels of this marker (Fig. 1B). In terms of risk of outcomes for miR-516a-5p, OS was marginally associated (HR = 1.851, 95% CI = [0.940-3.642], P = 0.0748), but DFS was not (HR = 1.239, 95% CI = [0.741-2.071], P = 0.4144) when controlling for the current prognostic factors, that is, age, stage, and MYCN (column "model 3" in Table 2). These hazard ratios were in the opposite direction compared to those of miR-487b. In addition, these two miRNAs were not significantly correlated: the percentage of high expression of miR-516a-5p was 70.9% and 75.3% in low and high expression of miR-487b, respectively (P = 0.482).
miR-487b and miR-516a-5p add prognostic information to the current prognostic factors Figure 1C showed a global significant difference (P < 0.0001 for OS and DFS) between the four survival curves defined by the combination of miR-487b and miR-516a-5p. In particular, we observed that the difference between miR-487b-based survival curves was not significant in patients with low miR-516a-5p expression (P = 0.8282 and P = 0.3652 for OS and DFS, respectively), while it was significant in patients with high miR-516a-5p expression (P < 0.0001 for OS and DFS). This suggests a modifier effect of miR-516a-5p on the prognostic value of miR-487b. The gain obtained with these two miRNAs may also be observed in the two subgroups defined by the standard risk stratification (low-and high risk). Within each risk group, expression of these two miRNAs allowed yet to identify patients with different prognoses (Fig. S5). The analysis of multivariable Cox models (reference model, on one hand, models including miR-487b and miR-516a-5p separately and together, on the other hand) measuring the adjusted associations of each prognostic factor with the hazard of survival outcomes is reported in Table 2. Each model corresponds to a separate column in the table. The models including both miR-487b and miR-516a-5p for OS (column "model 4") and including only miR-487b for DFS (column "model 2") fit the data considerably better than the reference model: likelihood ratio test P = 0.002 with lowest AIC = 629.002 for OS and P = 0.007 with lowest AIC = 862.645 for DFS compared to the reference model. In these two models, high expression of miR-487b is associated with a 69% decrease in the hazard of death (HR = 0.31, 95% CI = [0.13-0.75]) and 54% decrease in the hazard of disease-related event (HR = 0.46, 95% CI = [0.25-0.84]). High expression of miR-516a-5p is associated with a 98% increase in the hazard of death (HR = 1.98, 95% CI = [1.00-3.92]). In "model 5," a significant interaction term (P = 0.013 and 0.0005 for OS and DFS, respectively) leads to a better predictive ability representing by the lowest AIC (AIC = 624.851 and 853.590 for OS and DFS, respectively) compared to all other models. In Table 3 which reformulates model 5 reported in Table 2, miR-516a-5p was a significant predictor of OS (HR = 2.77, 95% CI = [1.25-6.15],

Disease-free surviving (%)
High miR-516a-5p/High miR-487b High miR-516a-5p/Low miR-487b Low miR-516a-5p/Low miR-487b Low miR-516a-5p/High miR-487b High miR-516a-5p/High miR-487b High miR-516a-5p/Low miR-487b Low miR-516a-5p/Low miR-487b Low miR-516a-5p/High miR-487b "miR-516a-5p high /miR-487b high ", "miR-516a-5p high /miR-487b low ", "miR-516a-5p low /miR-487b high ", and "miR-516a-5p low /miR-487b low ". In addition, we evaluated the discriminant ability of the two models retained above, that is, including both miR-516a-5p and miR-487b (AE interaction term) in comparing their time-dependent AUCs with that of the reference model. Table 4 represents the time-dependent AUC within the first 5 years (see Fig. S6 for the time-dependent AUCs from five models until 10 years). The addition of the two miRNAs (with or without interaction) gave better risk discrimination than did the reference model, as shown by larger AUC leading to a relative increase which ranges from 8% to 62%. This suggests the potential of these prognostic models to identify groups of patients with different risks of event. Although model 5 had slightly lower discriminant ability than model 4, it was considered as the final model given it showed the lowest AIC criterion (overall performance). From it, three risk groups with significant different prognoses were identified: worse, intermediate, and good prognosis in the groups 3, 2, and 1, respectively (Fig. S7).

Internal and external validations provide a prognostic model including two miRNAs plus interaction term
The order of risk groups derived from the derivation set was maintained in the internal and external validation. The categorization of PICV (Prognostic Index for Inter-  and 95% confidence interval (95% CI) of prognostic risk sets using group 1 as baseline were reported. The different risk groups were defined as follows: 1 = "miR-516a-5p high /miR-487b high ", 2 = "miR-516a-5p low /miR-487b high and low ", 3 = "miR-516a-5p high /miR-487b low ". nal Validation) according to its median and third quartile indicated good separation of three groups ( Fig. 2A) . This model also fitted well with similarity between the observed (Kaplan-Meier) and predicted survival curves (calibration attribute) ( Fig. 2A). We confirmed the three prognostic groups in the validation set (Fig. 2B) with a similar order between them even if the discrimination was less marked between the intermediate (group 2) and good (group 1) prognostic groups and more particularly for DFS (see Table S3 for PI equation and the cutoff defining the risk groups). It may be possibly explained by a better DFS in patients of the validation set making difficult to separate patients into risk sets of different prognoses. Another reason is the limited number of events in the validation set ( Table 1). The HRs of prognostic risk sets were HR 2/1 = 3.493, 95% CI = (0.636-19.164) and HR 3/1 = 11.629, 95% CI = (3.781-35.767) for OS; and HR 2/1 = 1.249, 95% CI = (0.438-3.562) and HR 3/1 = 6.553, 95% CI = (2.471-17.375) for DFS. In term of calibration, the inaccuracy prediction between observed and predicted survival curves may result by substantial different baseline survival function in the two populations (derivation and validation sets). This last point is, however, less important regarding the aim focusing on a risk score able to well-discriminate neuroblastoma patients. We also provided an update of regression coefficients of model 5 (adjusted by dataset) in carrying out a pooling of the two datasets. Column "Combined" of Table 3 shows similar regression coefficients.

Discussion
We have identified two parental imprinted miRNAs (miR-487b and miR-516a-5p) able to improve the current risk stratification of neuroblastoma based on age, stage, and MYCN amplification. The addition of these two markers with these three current prognostic factors improves the predictive ability and allows identifying new groups of patients with different prognoses both in terms of DFS and OS. On a basic point of view, these two markers belong to the two imprinted clusters, one maternally imprinted, C14MC, and the other paternally imprinted, C19MC, both clusters located at fragile sites of the human genome. Given the monoallelic expression of these clusters in normal cells, one can assume that any variation in expression of microRNAs of these loci (increase or decrease) may be instrumental for oncogenesis. As a matter of fact, high-risk neuroblastomas show differential underexpression of miR-487b sharply contrasting to the overexpression of miR-516a-5p. Very recently, the role of miR-487b as tumor suppressor gene was assigned: indeed, a loss of miR-487b expression was shown during the carcinogenesis of airway epithelial cells [34]. But so far regarding function of miR-516a-5p, nothing at all is known. Therefore, it would be attractive to consider miR-516a-5p as an actor playing an oncogenic role.
However, in considering the combinations of miR-487b/miR-516a-5p expressions in various risk subtypes that our study delineates herein, it seems to be more appropriate to speculate for miR-516a-5p a role of a finetuning regulator of neuroblastoma oncogenesis. Besides neuroblastoma tumors, alterations of expression of either C14MC [20,21] or C19MC [25] are frequently found in many other human cancers. In relationship with neural crest origin and more particularly in neuroectodermal brain tumors of children, overexpression of a set of microRNAs from the C19MC (related to an amplification of the locus) hallmarked aggressive primitive tumors [35]. As demonstrated over the last 5 years, miRNAs have been identified as interesting prognostic markers in neuroblastoma [14][15][16][17]. One important advantage of miRNAs is their high stability as compared to mRNA that limits the clinical routine application of gene expression [36,37]. De Preter and colleagues proposed a global signature based on 25 miRNAs [18]. These authors evaluated the discriminant ability of this signature and obtained AUC values of 0.74 for OS and 0.72 for DFS with an endpoint to 3 years. We obtained similar results with AUC at 3 years estimated to 0.71 and 0.67 for OS and DFS, respectively, using a parsimonious model into account six parameters to estimate (age, stage, MYCN, miR-487b, miR-516a-5p, and interaction between the two markers). Noteworthy, miR-487b and miR-516a-5p do not belong to the signature of De Preter and coworkers. However, no functional studies warrant identity of target genes, so far. Whether or not our two miRNAs and this signature would regulate similar key genes involved in pathways related to neuroblastoma oncogenesis remain to be determined. The discriminant ability of our model, estimated from the derivation set, was also confirmed by an external validation set of 101 patients, albeit at a lesser extent between intermediate-and low-risk groups. This may be likely explained by a relative homogeneity of these patients in terms of prognosis. This evidence allows foreseeing the generalization of the model.
In a daily clinical practice, the quantification of two miRNAs expression is realistic and feasible in an easy and rapid way using quantitative RT-PCR; not to mention this analysis can be done with a low number of cells such as that resulting from fine needle aspirations. Our results were obtained from tumor tissues with a malignant tumor cell content of ≥60%. In this respect, similarly to global genomic approaches (array-comparative genomic hybridization (CGH) and genes expression panel), the routine application of qRT-PCR analysis in clinical practice might be considerably limited by stromal cells con-  Each model was controlling for the current prognostic factors. Model 1: age + stage + mycn; Model 4: Model 1 + miR-487b + miR-516a-5p; Model 5: Model 1 + miR-487b + miR-516a-5p + miR-487b 9 miR-516a-5p. 1 For each time point, the status of patients was defined as 1 if T ≤ t and 0, T > t and compute the AUC at this time. and pan-genomic alterations constitute an algorithm that could improve the risk stratification of neuroblastoma remains to be investigated. Further studies on larger external cohorts are needed to definitively assess the role of these two markers in neuroblastoma and to investigate whether the association between the two miRNA from the C14MC and C19MC in the tumor progression of neuroblastoma may be extended to other specific miRNA from these clusters in various human cancers. Also and importantly, future studies will ascertain the impact of the combined miR-487b/miR-516a-5p expression on the genomic stratification neuroblastoma through CGH-array [7,38]. Such findings would open to a better treatment (therapeutic escalation or des-escalation).
In conclusion, our findings demonstrate that two markers, miR-487b and miR-516a-5p, enhance the accuracy of the current risk score, thus defining a better discrimination of neuroblastoma patients at diagnosis. Clinical use of this model may provide better therapeutic management of the disease by pediatric clinicians.

Supporting Information
Additional Supporting Information may be found in the online version of this article: Figure S1. Kaplan-Meier curves comparing overall survival (respectively, disease-free survival) in the derivation and validation sets with the log-rank test. Figure S2. Kaplan-Meier curves of expression of eight miRNAs of the 14q32.31 cluster-miR-380-5p, miR-494, miR-376c, miR-487b, miR-154, miR-134, miR-409-3p, and miR-410-for overall survival of the derivation set (n = 231). miRNAs expression levels were converted into discrete variables by discriminating the samples into two classes (high expression in blue and low expression in red), under or over the cutoff defined as the expression level of IGR-N-835 cell line taken as one. Figure S3. Kaplan-Meier curves of expression of eight miRNAs of the 14q32.31 cluster-miR-380-5p, miR-494, miR-376c, miR-487b, miR-154, miR-134, miR-409-3p, and miR-410-for disease-free survival of the derivation set (n = 231). miRNAs expression levels were converted into discrete variables by discriminating the samples into two classes (high expression in blue and low expression in red), under or over the cutoff defined as the expression level of IGR-N-835 cell line taken as one. Figure S4. Scatter-plots showing the statistical analysis of miR-516a-5p expression level relative to the expression level of IGR-N-835 cell line in the derivation set (n = 231) stratified in two groups, low-risk (blue, n = 151) versus high-risk (red, n = 80). Two-tailed Student's t-test indicates significance between the two subgroups. Black bars inside the plots indicate the mean value of miRNA expression levels. Figure S5. Kaplan-Meier survival curves for overall survival and disease-free survival between four groups defined by the combination of miR-487b and miR-516a-5p in (A) low and (B) high standard risk groups of the derivation set (n = 231). Figure S6. Area under the curve (AUC) of different prognostic models applied to the derivation set (n = 231). Figure S7. Comparing observed (Kaplan-Meier) survival curves for risk groups (solid lines) and model-based predicted mean survival curves (dashed lines) (A) for overall survival and (B) for disease-free survival from the derivation set. Below each figure, the hazard ratio and 95% confidence interval into brackets of prognostic risk sets using group 1 as baseline were reported. Table S1. Individual data of the 332 neuroblastoma patients (derivation and validation sets) and the expression levels of the two selected miRNAs using qRT-PCR. Table S2. Cox proportional hazard models with backward selection of miRNAs from the C14MC cluster. Table S3. Prognostic index (PI) for overall survival and disease-free survival estimated from the derivation set (n = 231) and cut-points used to define the risk groups in the validation set. The cut-points correspond to the median and third quartile of the distribution of the PI across patients from the derivation set.