Development and Validation of a 6-miRNA Prognostic Signature in Spinal Chordoma

Background: Published data have suggested a critical role for microRNA (miRNA) expression in chordoma progression. However, most of these studies focus on single miRNA and no multi-miRNA prognostic signature has been currently established for chordoma. In this study, we sought to develop and validate a 6-miRNA risk score (miRscore) model for survival prediction.

Methods: Medline, Embase, and Google scholar searches (from inception to July 20, 2018) were conducted to identify candidate miRNAs with prognostic value as per predefined criteria. Quantitative RT-PCR was used to measure miRNA levels in 114 spinal chordoma (54 in the training and 60 in the validation cohort) and 20 control specimens. Subsequently, the miRscore was built based on miRNAs data.

Results: Literature searches identified six prognostic miRNAs (miR-574-3p, miR-1237-3p, miR-140-3p, miR-1, miR-155, and miR-1290) with differential expression in tumor tissues. Bioinformatical analysis revealed an important regulatory role for miR-574-3p/EGFR signaling in chordoma and showed that the target genes of these prognostic miRNAs were mainly enriched in transcription regulation, protein binding and cancer-related pathways. In both cohorts, the miRscore was associated with surrounding muscle invasion by tumor and/or other aggressive features. The miRscore model well predicted local recurrence-free survival and overall survival, which remained after adjusting for other relevant covariates. Further time-dependent receiver operating characteristics analysis in the two cohorts found that the miRscore classifier had stronger prognostic power than known clinical predictors and improved the ability of Enneking staging to predict outcomes. Importantly, recursive-partitioning analysis of both samples combined separated patients into four prognostically distinct risk subgroups for recurrence and survival (both P < 0.001).

Conclusions: These data suggest the miRscore as a useful prognostic stratification tool in spinal chordoma and may represent an important step toward future personalized treatment of patients.

Introduction

Chordoma is a very rare and locally aggressive malignant mesenchymal neoplasm, which has been considered to arise from remnants of the embryonic notochord (1, 2). Despite the great advance in oncology field, the current treatment of choice for chordoma is still limited. Chordoma is resistance to traditional chemotherapy or radiotherapy (3). Therefore, radical surgical excision constitutes the most important therapeutic approach for chordoma, which has been reported to provide the optimal long-term disease control when combined with postoperative adjuvant radiation (4, 5). However, even with initial maximal resection, chordoma still has high local recurrence after surgery and 40–50% of patients may even develop remote metastasis (6), which poses a large challenge for the clinical management of this tumor entity. Considering the dismal prognosis of chordoma, it is highly imperative to develop better treatment strategy aiming at improving patient survival.

Previous studies have shown that epigenetic dysfunction plays a key role in chordoma development and progression (7, 8). To be in detail, it has been demonstrated that expression of several microRNAs (miRNAs) in chordoma tissues exhibits close association with patient outcomes (8), supporting the use of miRNAs as prognostic biomarkers in chordoma. However, given that tumor initiation is driven by a complicated molecular event, single miRNA may not provide optimal information on prognosis when compared with multi-miRNA signatures (9). To overcome this, researchers have attempted to establish miRNA-based system by simultaneously integrating multiple miRNAs data in outcome prediction of cancer patients (10–13). Currently, tissue-derived miRNA models have been created for several human cancers with good predictive performance (9, 14–18), some even display obvious advantage over traditional tumor-node-metastasis classification system for survival prediction (14, 18). However, there is currently a lack of studies investigating the role of miRNA-based system for chordoma prognostic risk stratification.

Clinically, the Enneking system is commonly used to determine disease staging and guide subsequent treatment for spinal chordoma patients. However, clinical outcome varies even among patients with the same Enneking staging, implying that this system is not adequate for prognosis. Presently, there is no recognized prognostic tool for spinal chordoma in clinics. Recently, Karhade et al. (19) developed a machine learning model to provide a quantitative prediction of 5-year survival for patients with spinal chordoma using the Surveillance, Epidemiology, and End Results database. Despite large sample size, the authors only considered clinical factors in their model without simultaneously incorporating tumor molecular features, which may not offer accurate prognostic information. In the present study, we aimed to develop and validate a 6-miRNA risk score (miRscore) for spinal chordoma outcome prediction in two independent cohorts. We also attempted to establish a miRscore-based decision tree to facilitate patient-individualized therapy for neurosurgeons.

Materials and Methods

Patients

This study included two patient cohorts. The training cohort consisted of 54 spinal chordoma patients who underwent surgical resection at our institute between June 2002 and April 2015. This cohort has been previously reported in our studies and patients were followed-up until September 2015 (20, 21). The validation cohort recruited 23 patients from our institute (during November 2015–December 2018 period) and 37 patients from Xiangya hospital (during February 2006–October 2018 period). These patients received curative tumor excision and follow-up information was updated in April 2019. Patients selection fulfilled comparable baseline, treatment and follow-up characteristics to that of the training population (Supplementary Table 1) (22). Sample size in validation group was estimated in advance by a power calculation (23), and minimum of 45 patients were required to show difference with a power of 80%, β error of 0.2, and α error of 0.05. Our validation cohort had 60 patients in total, which was adequate for analysis. For all patients, clinicopathologic data were obtained from medical records and evaluated as we previously detailed (20, 21). In the training cohort, expression of PD-L1 and Ki-67, as well as the Immunoscore pattern (low: I0-I2, high: I3-I4), was directly obtained from our prior data (20, 21). In the validation cohort, immunohistochemistry was conducted to produce the above data, in which the same antibodies, dilution, staining procedure, and evaluation criteria were used as those in the training cohort to ensure comparability (20, 21).The primary outcome parameters included local relapse-free survival (LRFS) and overall survival (OS). LRFS was calculated from the date of surgery to the first local recurrence, which was considered by clinical and imaging findings and further confirmed by histologic examinations of the second surgery (24). Similarly, OS was defined as the time from tumor resection to death related to any cause. Observations were censored when a patient was tumor-free (LRFS analysis) or alive (OS analysis).

Tissue Samples

Formalin-fixed paraffin-embedded (FFPE) tumor blocks of the 114 chordoma patients and 20 nucleus pulposus blocks (as controls) of disc herniated patients who underwent surgery during the same period were retrieved from the Department of Pathology. Control patients were randomly selected and matched for age and sex to cases from the training and validation cohort (25). Diagnosis of chordoma was made on hematoxylin and eosin-stained tissue sections by two neuropathologists per a previously published criterion (1). The study protocol was approved by our hospital ethical committee, and informed consent was obtained from each patient for publication of this study.

Selection of Prognostic miRNAs in Spinal Chordoma

We searched MEDLINE, EMBASE and Google scholar from the inception to July 20, 2018 to identify English language studies addressing the prognostic significance of miRNAs in chordoma. The keywords combinations used were (“miR” or “miRNA” or “microRNA”) and (“chordoma” or “chordomas”). Bibliographies of the included studies and existing systematic reviews were also manually checked for any additional relevant citations missed by the above search term. Specifically, the eligible studies should provide data analyzing the relationship between miRNA expression level and clinical outcomes (such as OS, LRFS, event-free survival or progression-free survival) of patients. In addition to this restriction, only literatures in which the studied miRNA was initially selected on the basis of microarray data either from authors’ group or previously published reports were considered for further evaluation. Studies were excluded based on the following criteria: studies focused on skull base chordoma, studies not involving human, studies with fewer than 20 patients, case studies, review studies, and in vitro and/or in vivo studies with no analysis of patient outcomes.

RNA Extraction and Quantitative RT−PCR

RNA extraction and quantitative RT−PCR were performed as we previously described (14, 26). Briefly, total RNA was isolated from 114 FFPE chordoma specimens and 20 FFPE nucleus pulposus samples using the mirVanaTM RNA isolation kit (Applied Biosystems, CA, USA) according to the manufacturer’s instructions. Total RNA (10 ng) was reversely transcribed into cDNA using the miRNA Reverse Transcription Kit (Qiagen, Dusseldorf, Germany). Quantitative PCR reactions were then conducted using a TaqMan^® Universal PCR Master Mix (Applied Biosystems, CA, USA). The relative expression of miRNAs was calculated by the 2^-ΔΔCT method. Each sample was analyzed in triplicate and U6 snRNA was used as the normalization control. As miR-155 and miR-1 have two different mature isoforms (miR-155-3p/5p and miR-1-3p/5p), the primers for these two miRNAs were determined as previously described (27, 28). The specific primers used are listed in Supplementary Table 2.

Bioinformatic Analysis of Prognostic miRNAs

We used four online complementary computational databases (miRTarBase, http://mirtarbase.mbc.nctu.edu.tw/php/index.php; miRDB, http://www.mirdb.org/; miRWalk2.0, http://zmf.umm.uni-heidelberg.de/apps/zmf/mirwalk2/; TargetScan, http://www.targetscan.org/) to predict the target genes of the six prognostic miRNAs. Genes were included for analysis only when they were jointly predicted by algorithms of at least three databases. After this, we conducted gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis to explore the main function and signaling pathways in which the annotated genes were involved. Finally, we constructed the miRNA-gene networks according to the regulatory relationships between miRNA and genes by using Cytoscape version 3.5.1 (http://www.cytoscape.org/), which is an open-source software tool developed and updated by researchers around the world to visualize biomedical networks involving proteins, genes and other interactions. We also retrieved protein-protein interaction (PPI) networks from STRING database (version 11.0, http://string-db.org).

Construction of 6 miRNA-Based Risk Score (miRscore)

The miRscore was established as previously suggested (17, 29). First, the univariate Cox analysis was performed to screen for candidate miRNAs that were significantly correlated with OS of spinal chordoma patients. Then, candidate miRNAs were fitted into a multivariable Cox regression analysis with OS as the dependent variable to determine the contribution of each miRNA to the outcome prediction. Finally, the risk score for each patient was derived as a linear combination of the expression data for each candidate miRNAs and the corresponding multivariable Cox regression coefficients.

Statistical Analysis

In both cohorts, quantitative data were presented as mean ± standard deviation and analyzed by Student’s t test or One-Way ANOVA test. Categorical data were analyzed by Chi-square test or Wilcoxon’s rank sum test, when appropriate. In the training cohort, the X-tile software (version 3.6.1, https://medicine.yale.edu/lab/rimm/research/software/, Rimm’s lab, Yale School of Medicine, New Haven, CT, USA) was applied to determine the cutoff points with the most significant split (the minimum P value approach) for continuous variables in survival analysis by regarding OS as the outcome parameter (30). To avoid overstating significance of the data, the P value from the log-rank test above was subsequently corrected as previously suggested (31). These cutoff values were directly used for subsequent analyses in the validation data. In the two patient cohorts, the Kaplan-Meier curve by univariate log-rank test was used to compare difference in survival probabilities between different subgroups. Similarly, multivariate Cox regression analysis was also performed on these data sets to see whether the miRscore model could independently predict patient survival, after controlling for other clinical predictors that were previously reported in literature (32–34) and/or significant in our univariate survival analysis. Time-dependent receiver operating characteristic (ROC) curve was used in both cohorts to evaluate the sensitivity and specificity of variables for survival prediction. Recursive-partitioning analysis (RPA) for creation of a decision tree was performed using the combined training and validation data to provide reference for clinical decision-making and individualized therapy (14). All statistical analyses were performed using R version 3.5.1 (R Foundation for Statistical Computing, Vienna, Austria). All tests were two-sided and a P value ≤ 0.05 was considered statistically significant.

Results

Clinical Characteristics

A total of 114 patients (54 in the training and 60 in the validation cohort) were included. Patient characteristics were recently communicated by our group (22) and summarized in Supplementary Table 1. Briefly, there were 35 males and 19 females in the training cohort. Of them, 43 were primary tumors and 11 were recurrent tumors. A total of 36 patients received Enneking inappropriate (EI) resection, and 18 patients underwent Enneking appropriate (EA) resection. By contrast, the validation cohort comprised 42 males and 18 females. 47 patients had primary diseases and 13 had recurrent ones. Twenty-four patients had EI resection, and 36 had EA resection. There were no significant differences for baseline characteristics between patients in training and validation cohort (Supplementary Table 1). All chordoma cases belonged to classic pathology type.

Identification of Six Prognostic miRNAs Associated With Survival of Spinal Chordoma Patients

A total of 122 studies were retrieved and the full-texts were fully reviewed independently by two authors (WH and MXZ) for eligibility (Supplementary Figure 1). After literature reviewing, six miRNAs (miR-574-3p, miR-1237-3p, miR-140-3p, miR-1, miR-155, and miR-1290) were identified to have significant association with spinal chordoma prognosis (Supplementary Figure 1) (20, 26–28, 35, 36). For each single miRNA, only one study reported its significant impact on spinal chordoma prognosis. Specifically, miR-140-3p (high expression in chordoma tissues compared to controls) was shown to have negative prognostic implication, while the other five miRNAs (low expression in chordoma tissues compared to controls) displayed positive association with survival of patients. These miRNAs were all initially screened by miRNA microarray and showed significantly differential expression in chordoma tissues compared with nucleus pulposus tissues (as controls) as analyzed by quantitative RT−PCR assay (Figure 1).

FIGURE 1

Figure 1 Quantitative expression levels of the six prognostic miRNAs by RT-PCR in the training, validation and control groups. All miRNAs were differentially expressed in tumor tissues compared with nucleus pulposus tissues (as controls).

Bioinformatic Analysis of the Six Prognostic miRNAs

To get deeper insights into the biological function of the six prognostic miRNAs, we performed GO and KEGG analysis, and also generated miRNA-mRNA target regulatory and protein-protein interaction networks. It should be noted that miR-155 and miR-1 have two different mature isoforms, we therefore included all of them in computational analysis in order to obtain more comprehensive outcomes. A total of 731 nonoverlapped genes targeted by these prognostic miRNAs were identified (Supplementary Figure 2). The most significant GO terms corresponding to these predicted genes included regulation of gene transcription and protein binding (Figure 2 and Supplementary Table 3). Pathway enrichment analysis of target genes revealed 45 pathways, including FoxO signaling, pancreatic cancer, neutrophin signaling, transcriptional misregulation in cancer, axon guidance, pathway in cancer, and other cancer-associated signaling pathways (Figure 2 and Supplementary Table 4). miRNA-mRNA regulatory network showed a potentially functional link between these prognostic miRNAs (Supplementary Figure 3). PPI network indicated a tight connection formed by these proteins with MAPK1, MAPK8, and EGFR as the potential hub genes (Supplementary Figure 4). Further integrated analysis of both the miRNA-mRNA and PPI networks data identified an important regulatory role for miR-574-3p-mediated signaling pathways in chordoma, especially the EGFR pathway (Figure 3).

FIGURE 2

Figure 2 Results of GO and KEGG pathways annotation analysis of the predicted genes targeted by the six prognostic microRNAs.

FIGURE 3

Figure 3 Identification of an important regulatory role for miR-574-3p-mediated signaling pathways (especially the EGFR/RAC1 pathway) in spinal chordoma by integrated analysis of both the miRNA-mRNA and protein-protein interaction networks data. miRNA, microRNA.

Construction and Description of the miRscore Model

Cutoff points for the 6 prognostic miRNAs are depicted in Supplementary Figure 5. Patients were separated into low (≤cutoff) and high (>cutoff) groups. Univariate Kaplan-Meier analysis with the log-rank test showed that the 6 miRNAs were significantly associated with patients’ OS in both cohorts (Supplementary Figures 6 and 7). This was also the case for these miRNAs in the survival analysis of LRFS, except for miR-1237-3p and miR-574-3p in the training data (Supplementary Figures 8 and 9). Further multivariate Cox regression analysis including the 6 miRNAs data was conducted and coefficients for each miRNA were obtained (Figure 4A). miRscore was then calculated as a linear combination of the expression data for each miRNA (Figures 4B, C) and its corresponding coefficient as described above, given as follows:

$$\text{miRscore}=\Sigma \left[lo{g}_2\right(Expressio{n}_{miRNA}\text{) }\times \hspace{0.17em}\text{Coefficien}{\text{t}}_{miRNA}]$$

FIGURE 4

Figure 4 (A) Coefficient of each prognostic miRNA obtained by Cox regression analysis. Heatmap shows the expression data of the six miRNAs used for the miRscore construction in the training (B) and validation cohort (C). Columns represent patients who were sorted descendingly according to their miRscore levels. (D) Determination of the cutoff point (-5.83) for miRscore in prognosis analysis with overall survival as the outcome parameter. Distribution of the miRscore in the training (E) and validation patients (F). The x-axis represents patients. The intersection point of black and yellow dotted line represents the miRscore cutoff classifying patients into high and low subgroups. (G, H) Covariate tracks show the clinical attributes of each patient in the training and validation cohort, respectively. Red boxes indicate that the miRscore group is enriched for the indicated attribute (P < 0.05, Chi-square test). miRscore, microRNA risk score.

For the purpose of prognosis analysis, miRscore was also classified into low (≤-5.83) and high (>-5.83) groups according to its threshold value (Figure 4D). Figures 4E, F illustrates the distribution of miRscore data among all patients. The average value for the miRscore in the training and validation cohort was -4.16 ± 3.39 and -4.36 ± 5.62, respectively.

Association Between miRscore and Clinicopathologic Parameters

In both cohorts, high miRscore was correlated with surrounding muscle invasion by tumors (Figures 4G, H and Supplementary Tables 5 and 6). In addition, tumors with high miRscore seemed to be more likely to have positive PD-L1 expression and advanced Enneking staging in the training cohort, although this correlation was not statistically significant (Figures 4G, H and Supplementary Tables 5 and 6). However, in the validation cohort, significant positive association was observed between the miRscore and grade, Ki-67 index, tumor hemorrhage, and the Enneking staging system (Figures 4G, H and Supplementary Tables 5 and 6).

Association Between miRscore Classifier and Patient Outcome

In both cohorts, high miRscore appeared to correlate with a high risk of death (Figures 5A, B). However, low miRscore was linked to increased tumor relapse (Figures 5A, B). The Kaplan-Meier curve showed that patients harboring high miRscore experienced worse LRFS and OS than those with low miRscore (Figures 5A, B). Furthermore, multivariate Cox regression analysis adjusting for other clinical variables revealed that the miRscore was an independent predictor of both LRFS and OS (Figures 5C, D).

FIGURE 5

Figure 5 (A) Distribution of recurrence and survival status in chordoma patients stratified by the miRscore data in the training cohort (top). Red dotted line represents the miRscore cutoff dividing patients into high and low groups. Kaplan-Meier curve of LRFS and OS according to low or high miRscore group (middle). (B) Distribution of recurrence and survival status in chordoma patients stratified by the miRscore data in the training cohort (top). Red dotted line represents the miRscore cutoff dividing patients into high and low groups. Kaplan-Meier curve of LRFS and OS according to low or high miRscore group (middle). (C) Multivariate Cox regression model for LRFS and OS of spinal chordoma patients in the training cohort. The boxes indicate the hazard ratio, and the horizontal lines represent 95% confidence intervals. (D) Multivariate Cox regression model for LRFS and OS of spinal chordoma patients in the validation cohort. miRscore, microRNA risk score; LRFS, local recurrence-free survival; OS, overall survival.

Comparison of the miRscore Model With Six Prognostic miRNAs, the Immunoscore, Enneking Staging System, and Other Clinicopathologic Variables in Outcome Prediction

In both cohorts, time-dependent ROC analysis found that the miRscore model displayed comparable prognostic value to Immunoscore system in predicting survival (Figure 6). Moreover, the miRscore classifier provided stronger prognostic power than the Enneking staging and other classical clinical factors in LRFS and OS prediction (Figure 6). Importantly, combined miRscore and Enneking staging enhanced the ability of each alone for outcome prediction (Figure 6).

FIGURE 6

Figure 6 Time-dependent sensitivity and specificity derived AUCs show the predictive performance of miRscore plus Enneking staging, miRscore, Immunoscore, Enneking staging and other known clinicopathological factors for LRFS and OS prediction in the training (left) and validation samples (right) at follow-up years 1 to 3 (for LRFS) or 1-5 (for OS). AUC, area under the curve; miRscore, microRNA risk score; LRFS, local recurrence-free survival; OS, overall survival.

Establishment of miRscore-Based Decision Tree for Clinical Decision-Making

Using combined data, RPA analysis identified four different risk groups for both tumor recurrence and patient death, including “low-risk”, “low-intermediate-risk”, “high-intermediate-risk” and “high-risk” groups (Figures 7A, B). Specifically, in analysis of LRFS, the 6-miRNA signature represented the strongest predictor together with tumor hemorrhage, type of surgery, Ki-67 index and tumor invading into surrounding muscle tissues or not (Figure 7A). The worst prognostic group referred to patients who received EI tumor resection and had high miRscore as well as tumor muscle invasion, while low miRscore patients without tumor hemorrhage had the best prognosis. Kaplan-Meier analysis by the log-rank test showed that the four risk groups significantly differed with regard to LRFS (Figure 7C). Subgroup analysis detected significant difference in LRFS time between low-risk group and low-intermediate-risk (median, 45 vs. 27 mo, P = 0.002), high-intermediate-risk (median, 45 vs. 19 mo, P < 0.001), or high-risk group (median, 45 vs. 10 mo, P < 0.001), and between low-intermediate-risk group and high-intermediate-risk (median, 27 vs. 19 mo, P = 0.001) or high-risk group (median, 27 vs. 10 mo, P < 0.001), as well as between high-intermediate-risk group and high-risk group (median, 19 vs. 10 mo, P < 0.001).

FIGURE 7

Figure 7 RPA tree for tumor recurrence (A) and survival (B) defined four different risk groups in the pooled data set (n = 114), including “low-risk”, “low-intermediate-risk”, “high-intermediate-risk” and “high-risk”. Each node displays the predicted probability of recurrence or death (color code low to high: green-red), the number of events for total patients, and the percentage of observations. Kaplan-Meier curves with log-rank test for LRFS (C) and OS (D) according to the four identified risk groups. LRFS, local recurrence-free survival; OS, overall survival. RPA, recursive-partitioning analysis.

Similarly, the RPA tree showed that the miRscore signature was also the strongest parameter in terms of OS, together with type of surgery, tumor invasion of the surrounding muscle tissue, Ki-67 index and primary or recurrent tumors (Figure 7B). The worst prognostic group included high miRscore patients who received EI tumor resection and had recurrent disease on admission, while patients with no muscle invasion by tumor as well as low miRscore and low Ki-67 index experienced the best prognosis. Subsequent Kaplan-Meier curve revealed that the four risk groups exhibited significant difference in relation to OS (Figure 7D). Statistically significant differences in OS were seen between low-risk group and low-intermediate-risk (median, not estimable [NE] vs. 78 mo, P = 0.003), high-intermediate-risk (median, NE vs. 27 mo, P < 0.001), or high-risk group (median, NE vs. 20 mo, P < 0.001), and between low-intermediate-risk group and high-intermediate-risk (median, 78 vs. 27 mo, P < 0.001) or high-risk group (median, 78 vs. 20 mo, P < 0.001), as well as between high-intermediate-risk group and high-risk group (median, 27 vs. 20 mo, P = 0.002).

Discussion

In this study, we developed a miRNA-based signature to predict the clinical outcome of spinal chordoma patients. We found that the miRscore system was closely associated with tumor surrounding muscle invasion and other aggressive clinicopathological features. The miRscore model accurately reflected survival independent of known prognostic parameters. Moreover, this signature showed good predictive performance in comparison with the Immunoscore system and other clinical factors, which could also improve predictive accuracy of the Enneking staging. Importantly, the decision tree integrating the miRscore with clinical predictors effectively separated patients into four prognostically distinct subgroups. These data suggest the use of miRscore signature as a clinical patient stratification tool in spinal chordoma and may facilitate future individualized therapy of patients.

Previous studies have demonstrated that the deregulation of miRNAs expression plays a crucial role in chordoma initiation and progression (7, 8). In support of this, our study revealed that the 6-miRNA signature was significantly correlated with adverse clinicopathological features and could independently predict survival of patients. Currently, the precise mechanism underlying the effect of miRNA signature on chordoma outcome remains undetermined. Published data have shown that the six prognostic miRNAs used for signature construction could influence chordoma prognosis by mediating expression of various target genes. For instance, Duan and colleagues found that low miR-1 expression in chordoma tissues promoted the proliferation and invasion of tumor cells by upregulation of slug (37), while miR-155 led to aggressive chordoma phenotype and poor patient survival mainly via regulating SOCS1 and TP53INP1 expression (28). Similarly, Wang et al. (38) proved that miR-1290 could suppress tumor cell proliferation and invasion by targeting Robo1 in chordoma. Our previous observations suggested that miR-140-3p, miR-1237-3p, and miR-574-3p could control biological behavior of chordoma cells possibly by modulating receptor tyrosine kinase pathways and expression of MMP2 and PD-L1, respectively (20, 26, 36). Collectively, these findings hint that our miRNA signature impacts chordoma outcome through complicated molecular events. This speculation can be corroborated by our subsequent bioinformatic analysis, showing that target genes of the six miRNAs were involved in various GO terms and signaling pathways related to cancer.

In addition, our computational analysis also disclosed a key regulatory role for miR-574-3p/EGFR pathway in chordoma. We previously reported that miR-574-3p could influence immune profiles of chordoma possibly by targeting PD-L1 (20). Moreover, recent data have provided evidence for EGFR/Ras/Mek signaling to control PD-L1 expression in cancers (39, 40). Taken together, these outcomes further highlight the clinical significance of immune microenvironment in chordoma (21, 41, 42) and suggest a complete miR-574-3p/EGFR/Ras/Mek/PD-L1 signaling network involved in chordoma development. This finding deserves further investigation and may have implication for sole or combination immunotherapy approaches for chordoma treatment.

Our study found that the miRscore system could effectively portend chordoma outcome independent of known parameters. Moreover, this miRNA-based signature had better prognostic power than the traditional Enneking system and also possessed comparable predictive ability to that of the Immunoscore system. Importantly, combined miRscore and Enneking staging model significantly improved predictive accuracy of each alone in both cohorts. Altogether, these results indicate that adding miRNA expression data into chordoma prognostic grading system is necessary in order to produce a reliable disease model and then guide therapeutic optimization. Preceding reports have suggested that the Enneking system as an indicator of anatomic features of tumors is not adequate to inform chordoma prognosis (21). Further, it has been determined that the immune microenvironment reflective of tumor bioimmunological properties plays a pivotal role in chordoma progression (20, 21, 41). Besides these data, recent studies have implied that radiological findings on magnetic resonance imaging are also related to chordoma outcome (43, 44). As the miRscore signature represents different molecular features from those above, future prognostic model considering all these aspects may prove a more accurate and reliable tool for risk stratification in spinal chordoma.

Interesting, we also found that the miRscore-based decision tree could effectively split patients into four different risk subgroups, with significantly distinct survival. Of note, in this setting, the miRscore was the most important element among other clinical parameters for both disease recurrence and death. These data reiterate the clinical relevance of the 6-miRNA signature in spina chordoma and offer an approach for more detailed, clinically meaningful stratification of patients, which can be useful in identifying cases with advanced disease and then making treatment decisions timely.

Limitations

Additional studies are needed to optimize the miRscore calculation for clinical utility. Besides, it should be noted that parts of our validation datasets were from the same institute as the training cohort, which might compromise the generalization of our findings. We did this because this intervention could enable us to provide more reliable results by analyzing the current validation data as enough sample size was reached to obtain statistically significant difference. But before translating the miRscore model into clinical practice, prospective external validation studies with large sample sizes especially from a different institute are still needed to corroborate our current data in the future. In addition, published data have suggested that many other miRNAs may impact patient outcomes in spinal chordoma (7). Our future study will further evaluate the eligibility for adding other miRNAs (such as miR-219-5p and miR-31) into the miRscore model and assessed their effect on spinal chordoma prognosis (45, 46). Finally, our study is a correlation study in nature. More studies will be required to clarify how the miRscore signature affects clinical outcome of patients.

Conclusions

The presented study identified a six-miRNA molecular signature that was associated with aggressive clinicopathological features and survival of spinal chordoma patients. The miRscore model had good predictive performance and also compensated the deficiency of Enneking system for survival prediction. The integration of the miRscore with known prognostically clinical parameters in decision tree allowed the definition of four subgroups with significantly different survival. These data suggest the miRscore as a useful prognostic stratification tool in spinal chordoma and may represent an important step toward future personalized treatment of patients.

Data Availability Statement

All datasets presented in this study are included in the article/Supplementary Material.

Ethics Statement

The study protocol was approved by the Institutional Review Board at the First Affiliated Hospital, University of South China, Hunan, China. Written informed consent was obtained from each patient for publication of this study.

Author Contributions

All authors participated in data acquisition. WH, G-HL, JL, and M-XZ contributed to the conception and design of the study. WH, Y-GY, W-JW, Z-HO, BW, and M-XZ did the data analysis and interpretation. WH, X-LL, T-LZ, and X-BW performed the experiments. WH, Z-HO, T-LZ, X-BW, and M-XZ contributed to drafting and revision of the manuscript. All authors contributed to the article and approved the submitted version.

Funding

This work was supported by the National Natural Science Foundation of China (81802211 to X-BW and 81871821 to JL), Natural Science Foundation of Hunan Province (2019JJ50542 to T-LZ and 2018JJ3738 to X-BW), and Project for Clinical Research of Hunan Provincial Health Commission (20201978 to T-LZ and 20201956 to M-XZ).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank Dr. Yi Jiang and Dr. Xiao-Ling She from Department of Pathology, The Second Xiangya Hospital, Central South University for pathological analysis of the study.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fonc.2020.556902/full#supplementary-material

Abbreviations

miRNA, microRNA; miRscore, miRNA risk score; FFPE, formalin-fixed paraffin-embedded; LRFS, local relapse-free survival; OS, overall survival; GO, gene ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; ROC, receiver operating characteristic; RPA, Recursive-partitioning analysis; EA, Enneking appropriate; EI, Enneking inappropriate; PD-L1, programmed cell death-1 ligand 1; EGFR, epidermal growth factor receptor; PPI, protein-protein interaction; NE, not estimable.

References

1. Walcott BP, Nahed BV, Mohyeldin A, Coumans JV, Kahle KT, Ferreira MJ. Chordoma: current concepts, management, and future directions. Lancet Oncol (2012) 13(2):e69–76. doi: 10.1016/S1470-2045(11)70337-0

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Samson IR, Springfield DS, Suit HD, Mankin HJ. Operative treatment of sacrococcygeal chordoma. A review of twenty-one cases. J Bone Joint Surg Am (1993) 75(10):1476–84. doi: 10.2106/00004623-199310000-00008

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Colia V, Stacchiotti S. Medical treatment of advanced chordomas. Eur J Cancer (2017) 83:220–8. doi: 10.1016/j.ejca.2017.06.038

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Dial BL, Kerr DL, Lazarides AL, Catanzano AA, Green CL, Risoli T Jr, et al. The Role of Radiotherapy for Chordoma Patients Managed with Surgery: Analysis of the National Cancer Database. Spine (Phila Pa 1976) (2020) 45(12):E742–51. doi: 10.1097/BRS.0000000000003406

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Stacchiotti S, Gronchi A, Fossati P, Akiyama T, Alapetite C, Baumann M, et al. Best practices for the management of local-regional recurrent chordoma: a position paper by the Chordoma Global Consensus Group. Ann Oncol (2017) 28(6):1230–42. doi: 10.1093/annonc/mdx054

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Fuchs B, Dickey ID, Yaszemski MJ, Inwards CY, Sim FH. Operative management of sacral chordoma. J Bone Joint Surg Am (2005) 87(10):2211–6. doi: 10.2106/JBJS.D.02693

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Gill CM, Fowkes M, Shrivastava RK. Emerging Therapeutic Targets in Chordomas: A Review of the Literature in the Genomic Era. Neurosurgery (2020) 86(2):E118–23. doi: 10.1093/neuros/nyaa008

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Yu X, Li Z. Epigenetic deregulations in chordoma. Cell Prolif (2015) 48(5):497–502. doi: 10.1111/cpr.12204

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Schmidt L, Fredsøe J, Kristensen H, Strand SH, Rasmussen A, Høyer S, et al. Training and validation of a novel 4-miRNA ratio model (MiCaP) for prediction of postoperative outcome in prostate cancer patients. Ann Oncol (2018) 29(9):2003–9. doi: 10.1093/annonc/mdy243

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Li X, Yu X, He Y, Meng Y, Liang J, Huang L, et al. Integrated Analysis of MicroRNA (miRNA) and mRNA Profiles Reveals Reduced Correlation between MicroRNA and Target Gene in Cancer. BioMed Res Int (2018) 2018:1972606. doi: 10.1155/2018/1972606

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Falzone L, Scola L, Zanghi A, Biondi A, Di Cataldo A, Libra M, et al. Integrated analysis of colorectal cancer microRNA datasets: identification of microRNAs associated with tumor development. Aging (Albany NY) (2018) 10(5):1000–14. doi: 10.18632/aging.101444

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Falzone L, Lupo G, La Rosa GRM, Crimi S, Anfuso CD, Salemi R, et al. Identification of Novel MicroRNAs and Their Diagnostic and Prognostic Significance in Oral Cancer. Cancers (Basel) (2019) 11(5):610. doi: 10.3390/cancers11050610

CrossRef Full Text | Google Scholar

13. Candido S, Lupo G, Pennisi M, Basile MS, Anfuso CD, Petralia MC, et al. The analysis of miRNA expression profiling datasets reveals inverse microRNA patterns in glioblastoma and Alzheimer’s disease. Oncol Rep (2019) 42(3):911–22. doi: 10.3892/or.2019.7215

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Hess J, Unger K, Maihoefer C, Schüttrumpf L, Wintergerst L, Heider T, et al. A Five-MicroRNA Signature Predicts Survival and Disease Control of Patients with Head and Neck Cancer Negative for HPV Infection. Clin Cancer Res (2019) 25(5):1505–16. doi: 10.1158/1078-0432.CCR-18-0776

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Lim EL, Trinh DL, Ries RE, Wang J, Gerbing RB, Ma Y, et al. MicroRNA Expression-Based Model Indicates Event-Free Survival in Pediatric Acute Myeloid Leukemia. J Clin Oncol (2017) 35(35):3964–77. doi: 10.1200/JCO.2017.74.7451

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Bagnoli M, Canevari S, Califano D, Losito S, Maio MD, Raspagliesi F, et al. Development and validation of a microRNA-based signature (MiROvaR) to predict early relapse or progression of epithelial ovarian cancer: a cohort study. Lancet Oncol (2016) 17(8):1137–46. doi: 10.1016/S1470-2045(16)30108-5

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Bai F, Zhou H, Ma M, Guan C, Lyu J, Meng QH. A novel RNA sequencing-based miRNA signature predicts with recurrence and outcome of hepatocellular carcinoma. Mol Oncol (2018) 12(7):1125–37. doi: 10.1002/1878-0261.12315

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Yang Y, Qu A, Zhao R, Hua M, Zhang X, Dong Z, et al. Genome-wide identification of a novel miRNA-based signature to predict recurrence in patients with gastric cancer. Mol Oncol (2018) 12(12):2072–84. doi: 10.1002/1878-0261.12385

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Karhade AV, Thio Q, Ogink P, Kim J, Lozano-Calderon S, Raskin K, et al. Development of Machine Learning Algorithms for Prediction of 5-Year Spinal Chordoma Survival. World Neurosurg (2018) 119:e842–7. doi: 10.1016/j.wneu.2018.07.276

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Zou MX, Guo KM, Lv GH, Huang W, Li J, Wang XB, et al. Clinicopathologic implications of CD8(+)/Foxp3(+) ratio and miR-574-3p/PD-L1 axis in spinal chordoma patients. Cancer Immunol Immunother (2018) 67(2):209–24. doi: 10.1007/s00262-017-2080-1

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Zou MX, Lv GH, Wang XB, Huang W, Li J, Jiang Y, et al. Clinical Impact of the Immune Microenvironment in Spinal Chordoma: Immunoscore as an Independent Favorable Prognostic Factor. Neurosurgery (2019) 84(6):E318–33. doi: 10.1093/neuros/nyy274

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Zou MX, Pan Y, Huang W, Zhang TL, Escobar D, Wang WB, et al. A four-factor immune risk score signature predicts the clinical outcome of patients with spinal chordoma. Clin Transl Med (2020) 10(1):224–37. doi: 10.1002/ctm2.4

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Zhao S, Su Y, Duan J, Qiu Q, Ge X, Wang A, et al. Radiomics signature extracted from diffusion-weighted magnetic resonance imaging predicts outcomes in osteosarcoma. J Bone Oncol (2019) 19:100263. doi: 10.1016/j.jbo.2019.100263

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Meng T, Yin H, Li B, Li Z, Xu W, Zhou W, et al. Clinical features and prognostic factors of patients with chordoma in the spine: a retrospective analysis of 153 patients in a single center. Neuro Oncol (2015) 17(5):725–32. doi: 10.1093/neuonc/nou331

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Zou MX, Lv GH, Li J, She XL, Jiang Y. Upregulated human telomerase reverse transcriptase (hTERT) expression is associated with spinal chordoma growth, invasion and poor prognosis. Am J Transl Res (2016) 8(2):516–29.

PubMed Abstract | Google Scholar

26. Zou MX, Huang W, Wang XB, Li J, Lv GH, Wang B, et al. Reduced expression of miRNA-1237-3p associated with poor survival of spinal chordoma patients. Eur Spine J (2015) 24(8):1738–46. doi: 10.1007/s00586-015-3927-9

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Duan Z, Shen J, Yang X, Yang P, Osaka E, Choy E, et al. Prognostic significance of miRNA-1 (miR-1) expression in patients with chordoma. J Orthop Res (2014) 32(5):695–701. doi: 10.1002/jor.22589

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Osaka E, Kelly AD, Spentzos D, Choy E, Yang X, Shen JK, et al. MicroRNA-155 expression is independently predictive of outcome in chordoma. Oncotarget (2015) 6(11):9125–39. doi: 10.18632/oncotarget.3273

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Hou Y, Yu Z, Tam NL, Huang S, Sun C, Wang R, et al. Exosome-related lncRNAs as predictors of HCC patient survival: a prognostic model. Am J Transl Res (2018) 10(6):1648–62.

PubMed Abstract | Google Scholar

30. Camp RL, Dolled-Filhart M, Rimm DL. X-tile: a new bio-informatics tool for biomarker assessment and outcome-based cut-point optimization. Clin Cancer Res (2004) 10(21):7252–9. doi: 10.1158/1078-0432.CCR-04-0713

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Altman DG, Lausen B, Sauerbrei W, Schumacher M. Dangers of using “optimal” cutpoints in the evaluation of prognostic factors. J Natl Cancer Inst (1994) 86(11):829–35. doi: 10.1093/jnci/86.11.829

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Bakker SH, Jacobs WCH, Pondaag W, Gelderblom H, Nout RA, Dijkstra PDS, et al. Chordoma: a systematic review of the epidemiology and clinical prognostic factors predicting progression-free and overall survival. Eur Spine J (2018) 27(12):3043–58. doi: 10.1007/s00586-018-5764-0

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Bettegowda C, Yip S, Lo SL, Fisher CG, Boriani S, Rhines LD, et al. Spinal column chordoma: prognostic significance of clinical variables and T (brachyury) gene SNP rs2305089 for local recurrence and overall survival. Neuro Oncol (2017) 19(3):405–13. doi: 10.1093/neuonc/now156

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Zhou J, Sun J, Bai HX, Huang X, Zou Y, Tan X, et al. Prognostic Factors in Patients With Spinal Chordoma: An Integrative Analysis of 682 Patients. Neurosurgery (2017) 81(5):812–23. doi: 10.1093/neuros/nyx081

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Wang Y, Chen K, Chen H, Zhang K, Lu J, Mao H, et al. Low expression of miRNA-1290 associated with local invasion and recurrence in sacral chordoma. Int J Clin Exp Pathol (2017) 10(11):10934–40.

PubMed Abstract | Google Scholar

36. Zou MX, Huang W, Wang XB, Lv GH, Li J, Deng YW. Identification of miR-140-3p as a marker associated with poor prognosis in spinal chordoma. Int J Clin Exp Pathol (2014) 7(8):4877–85.

PubMed Abstract | Google Scholar

37. Osaka E, Yang X, Shen JK, Yang P, Feng Y, Mankin HJ, et al. MicroRNA-1 (miR-1) inhibits chordoma cell migration and invasion by targeting slug. J Orthop Res (2014) 32(8):1075–82. doi: 10.1002/jor.22632

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Wang B, Zhang K, Chen H, Lu J, Wu G, Yang H, et al. miR-1290 inhibits chordoma cell proliferation and invasion by targeting Robo1. Transl Cancer Res (2019) 8(2):542–51. doi: 10.21037/tcr.2019.03.13

CrossRef Full Text | Google Scholar

39. Shen X, Zhang L, Li J, Li Y, Wang Y, Xu ZX. Recent Findings in the Regulation of Programmed Death Ligand 1 Expression. Front Immunol (2019) 10:1337. doi: 10.3389/fimmu.2019.01337

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Cha JH, Chan LC, Li CW, Hsu JL, Hung MC. Mechanisms Controlling PD-L1 Expression in Cancer. Mol Cell (2019) 76(3):359–70. doi: 10.1016/j.molcel.2019.09.030

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Zhou J, Jiang Y, Zhang H, Chen L, Luo P, Li L, et al. Clinicopathological implications of TIM3(+) tumor-infiltrating lymphocytes and the miR-455-5p/Galectin-9 axis in skull base chordoma patients. Cancer Immunol Immunother (2019) 68(7):1157–69. doi: 10.1007/s00262-019-02349-1

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Thanindratarn P, Dean DC, Nelson SD, Hornicek FJ, Duan Z. Advances in immune checkpoint inhibitors for bone sarcoma therapy. J Bone Oncol (2019) 15:100221. doi: 10.1016/j.jbo.2019.100221

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Zuckerman SL, Amini B, Lee SH, Rao G, Tatsui CE, Rhines LD. Predictive Value of Preoperative Magnetic Resonance Imaging Findings for Survival and Local Recurrence in Patients Undergoing En Bloc Resection of Sacral Chordomas. Neurosurgery (2019) 85(6):834–42. doi: 10.1093/neuros/nyy578

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Wei W, Wang K, Liu Z, Tian K, Wang L, Du J, et al. Radiomic signature: A novel magnetic resonance imaging-based prognostic biomarker in patients with skull base chordoma. Radiother Oncol (2019) 141:239–46. doi: 10.1016/j.radonc.2019.10.002

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Wei W, Zhang Q, Wang Z, Yan B, Feng Y, Li P. miR-219-5p inhibits proliferation and clonogenicity in chordoma cells and is associated with tumor recurrence. Oncol Lett (2016) 12(6):4568–76. doi: 10.3892/ol.2016.5222

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Bayrak OF, Gulluoglu S, Aydemir E, Ture U, Acar H, Atalay B, et al. MicroRNA expression profiling reveals the potential function of microRNA-31 in chordomas. J Neurooncol (2013) 115(2):143–51. doi: 10.1007/s11060-013-1211-6

PubMed Abstract | CrossRef Full Text | Google Scholar