miR-17-5p-CXCL14 axis related transcriptome profile and clinical outcome in diffuse gliomas

ABSTRACT

GBM tissues are comprised of not only tumor cells but also tumor-associated nontumor cells, such as stromal cells and immune cells, which dilute the purity of glioma cells and function in glioma biology. However, the roles of miRNAs in modulating glioma purity are not clarified. In total, 838 glioma samples with transcriptome data, including 537 RNAseq data from TCGA project and 301 microarray data from Chinese Glioma Genome Atlas (CGGA project), were recruited into our investigation. Tumor purity, molecular subtypes and IDH status were also available. R language was employed as the main tool for statistical analysis and graphical work. Screening miRNA profiling and paired TCGA samples’ transcriptome data demonstrates that miR-17-5p expression harbors the most significant positive correlation with glioma purity among all miRNAs. CXCL14 shows robust negative correlation with miR-17-5p expression in TCGA and CGGA dataset. miR-17-5p directly targets CXCL14 and functions as a tumor-suppressor of GBM. CXCL14 showed lower expression in proneural subtype and may contribute as a potential marker for proneural subtype in glioma. Genes markedly correlated with CXCL14 are involved in essential functions associated with anti-tumor immune process. CXCL14 has a strong correlation with immune(T cells, Monocytic lineage and Neutrophils) and Fibroblasts within glioma environment. miR-17-5p and CXCL14 exhibited predictive values for high-grade glioma(HGG) patients: Higher miR-17-5p indicated significantly longer survival while lower CXCL14 indicated longer survival. Our results highlight the importance of the miR-17-5p-CXCL14 axis in regulating key steps of anti-tumor immune process and may serve as potential targets of immune treatments for gliomas.

KEYWORDS:

• Glioma
• immune response
• CXCL14
• miR-17-5p

Introduction

Glioblastoma multiform(GBM), a grade IV glioma is the most common and lethal type of tumor in the central nervous system (CNS) for which there is no available curative treatment.^1-3 The incidence of GBM is approximately 3 new cases per 100 000 annually in Europe and North America.⁴ Similar to other types of cancers, increasing age is relevant with higher incidence of GBM, with a median age of 64 years at diagnosis.⁴ Despite the continual accumulation in the wealth of knowledge of GBM pathogenesis and therapeutic treatments, patient prognosis is poor for newly diagnosed GBM, with a median survival of 12–15 months postdiagnosis.^5,6 Typical characters of glioblastoma include persistent proliferation, strong invasive capacity and aggressive angiogenesis within the tumour core.⁷ Malignant solid tumour tissues are comprised of not only tumour cells but also tumour-associated normal epithelial and stromal cells, immune cells and endothelium cells.⁸ Stromal cells have been proven to play important roles in tumour growth, disease progression and drug resistance. Stromal cells within GBM environment also could stimulate human GBM cell growth through the activation of extracellular signaling.⁹ Meanwhile, tumor-associated macrophages (TAMs), either of peripheral origin or representing brain-intrinsic microglia are recruited to the glioma environment and could facilitate immune functions. In this manner, TAMs could support stroma for GBM cell expansion and invasion.¹⁰ These nontumor cells play important roles in glioma development and progression.⁷

microRNAs (miRNAs) are an abundant class of endogenously expressed 18–25 nucleotide noncoding RNAs.^11,12 These molecules can bind to the target messenger RNA and induce translational silencing or degradation (mRNA), thereby inhibiting gene expression.¹³ Various studies have begun to map out the expression profiling and roles of miRNAs in GBM, aiming to garner new insight for treating this insidious disease. However, current studies pay little attention to the miRNAs expression pattern within GBM environments under various purities.

Here, transcriptomic data of TCGA (the United States) and CGGA(China) were collected for analysis. The ESTIMATE algorithm⁸ were employed to calculate tumor purity. The landscape of differently expressed miRNAs under increasing purities was systematically depicted. miR-17-5p expression showed the most significant positive correlation with glioma purity among all miRNAs.

Immunotherapy research in cancers has sparked huge scientific interest and made great progress in the past few years.^14,15 Chemokines and their receptors have crucial roles in inflammatory human diseases, and efforts have been made to check the possibility of combining chemokine-based therapies with current cancer immunotherapies.^16,17 CC and CXC are two major chemokine classes, among which there are CXCL12 [also known as stromal cell-derived factor 1] and CXCL14 (originally designated BRAK).¹⁸ A clinical trial of combinatorial immunotherapy with the CXCR4 peptide antagonist is now undergoing since anti-CXCR4 or anti-CXCL12 antibodies could reduce tumour growth and prevent tumour extravasation in preclinical models.^17,19–23 CCL2 blocking synergistically promotes the cancer vaccine response in animal models of lung cancer and mesothelioma.²⁴ A phase Ib trial of a CCR2 antagonist was also recently initiated for patients with non-resectable pancreatic cancer.¹⁷ Based on current findings, directly targeting both pro-tumour and antitumour chemokine–chemokine receptor signaling pathways in combination with other immunotherapies promises to increase the efficacy of the immunotherapies for cancer treatment. Here, we evaluated the expression status of CXCL14 which was validated as a direct target of miR-17-5p and related biological process by analyzing mRNA profiling data across two databases including a total of 838 glioma samples. To our best knowledge, this is the largest and comprehensive study characterizing CXCL14 expression in whole grade glioma tumor masses.

Results

Landscape of differentially expressed miRNAs in association with glioma purity

As reported previously, the ESTIMATE algorithm⁷ was performed to assess glioma purity using TCGA transcriptome data, as well as ESTIMATE stromal and immune scores, which all exhibited moderate correlation with morphology evaluation. Then, we filtered paired GBM samples from TCGA miRNA microarrays dataset and calculated the correlation of all individual miRNAs with the tumor purity that was evaluated through transcriptome data. A total of 534 miRNAs were analyzed by Spearman correlation analysis and demonstrated 5 miRNAs (Spearman |R|> 0.5), which indicated a highly positive or negative correlation with glioma purity in the TCGA dataset (Supplementary Table 1 and Figure 1A). High purity GBMs were more likely to belong to proneural subtype (Figure 1A). To the contrast, low purity GBMs were enriched in mesenchymal subtype (Figure 1A). With regard to genomic alterations, the IDH mutation tended to indicate higher glioma purity although it did not reach significance (Figure 1A). As shown in Figure 1B-D, a robust positive correlation between the miR-17-5p expression level with glioma purity and a negative correlation with ESTIMATE immune and stromal score were determined, respectively. Meanwhile, we also validated the correlation level of miR-17-5p expression with glioma purity using the CGGA dataset. miR-17-5p expression also tended to positively correlate with GBM purity although the P value was not significant (Figure 1E-G).

Figure 1. Landscape of miRNAs, clinical and molecular characteristics in association with glioma purity.

A, TCGA RNAseq sets were arranged in an increasing order of GBM purity. The relationship between GBM purity and patients’ characteristics was evaluated (a, The distribution of GBM purity among several groups was assessed using one-way ANOVA.b, The distribution of GBM purity was assessed using the Student t test between two groups.) B-D, Correlation of miR-17-5p expression with GBM purity(B), ESTIMATE immune score(C) and ESTIMATE stromal score(D) in the TCGA cohort of GBM patients using Spearman analysis. E-G, Correlation of miR-17-5p expression with GBM purity(E), ESTIMATE immune score(F) and ESTIMATE stromal score(G) in the CGGA cohort of GBM patients using Spearman analysis.

Mir-17-5p expression is upregulated in glioblastoma and a potential marker for proneural subtype

In TCGA cohort, glioblastoma showed the higher miR-17-5p expression when compared to normal brain tissues (Student t test, P value< 0.0001; Figure 2A). When we examined miR-17-5p expression in CGGA cohort, we found that glioblastoma showed the highest miR-17-5p expression when compared to grade II and grade III glioma (Figure 2B). These findings suggested that as glioma neoplasia progresses, there is a trend towards increased expression of miR-17-5p.

Figure 2. miR-17-5p expression was upregulated in glioblastoma and downregulated in proneural subtype.

A, miR-17-5p expression in TCGA dataset. B, miR-17-5p expression in CGGA dataset. C-D, miR-17-5p expression in molecular subtypes. E-F, predictive value for proneural subtype.

To investigate the molecular relevance between miR-17-5p and GBMs, we determined the distribution of miR-17-5p expression in different molecular subtypes defined by TCGA dataset. As shown in Figure 2C, miR-17-5p was significantly upregulated in proneural subtype than other subtypes in TCGA dataset. In CGGA dataset, it also showed apparent trend although not significant (Figure 2D). This result enlightened us that miR-17-5p may serve as a biomarker for proneural subtype. ROC curves for miR-17-5p expression and proneural subtype in GBMs were analyzed and under curve area is 76.51.1% and 63.24% in TCGA and CGGA dataset, respectively (Figure 2E and F).

Mir-17-5p directly targets CXCL14

To understand the potential role and mechanism of miR-17-5p in gliomas, we adopted the bioinformatic algorithms and mRNA profiling of glioma patients to identify potential target genes of miR-17-5p. A list of genes that were negatively correlated with miR-17-5p expression (Spearman r < −0.35) in both TCGA and CGGA glioma patients was produced. Overlapping this list with the candidates generated by 5 prediction algorithms (miRanda, RNAhybrid, miRWalk, RNA22 and Targetscan), we narrowed the miR-17-5p target candidates down to 2 genes (Figure 3A). To explore whether miR-17-5p directly targets these 4 genes, 3′ -UTR sequences containing putative binding sites of WT or mut were cloned into the pMIR-REPORT vector. miR-17-5p transfected U251 and N3 cells indicated a significant reduction of luciferase activities when cells were cotransfected with reporter plasmid (CXCL14–WT) (Figure 3B). However, miR-17-5p overexpression did not affect the luciferase activities of CXCL14 3′ -UTR mut reporter (Figure 3B). Meanwhile, miR-17-5p transfected cells indicated no significant change when cells were cotransfected with the reporter plasmids (GPR132–WT) (Figure 3C). Then, Western blotting analysis was conducted to determine whether CXCL14 expression was indeed regulated by miR-17-5p at the protein level. The expression of CXCL14 protein was downregulated in miR-17-5p treated cells, but increased in cells transfected with anti-miR-17-5p inhibitor (Figure 3D).

Figure 3. MiR-17-5p directly targets CXCL14.

A, Venn diagram indicating 4 genes computationally predicted to be putative miR-17-5p target by 5 prediction algorithms and negatively associated with miR-17-5p in both TCGA and CGGA dataset(Left). Sites of miR-17-5p seed matches in the CXCL14 3ʹ-UTR and a schematic diagram of the mutated CXCL14 3ʹ-UTR sequence(Right). B, Cells co-transfected with the luciferase/CXCL14 3ʹ-UTR reporter vector (WT or Mutant) and either miR-NC or miR-17-5p. C, Cells co-transfected with the luciferase reporter vector (WT or Mutant) and either miR-NC or miR-17-5p. D, CXCL14 protein expression of cells transfected with miR-17-5p mimics or miR-17-5p inhibitors E-F, CXCL14 expression in TCGA (E) dataset and CGGA (F) dataset according to IDH status. G-H, CXCL14 expression in molecular subtypes in TCGA (G) dataset and CGGA (H) data set. I-J, CXCL14 expression predictive value for proneural subtype. K-L, Survival analysis for CXCL14.

CXCL14 is downregulated in IDH mutant glioma and a potential marker for proneural molecular subtype

Subsequently, the CXCL14 expression was evaluated in gliomas with various purities. As expected, CXCL14 expression levels were significantly higher in glioma patients with low tumor purity than in those with high tumor purity according to the TCGA and CGGA dataset (Supplementary Figure 1). IDH1/2 mutations are the most common mutations in glioma. The sub-classifier of IDH mutation demonstrated that IDH mutant-type glioma indicated universally lower expression of CXCL14 than those of IDH wild-type(WT) glioma, across different grades, though no statistical significance was observed in one group. This reminded us that CXCL14-related biological process were more prevalent in IDH wild-type gliomas (Figure 3E-F). To analyze the molecular relevance between CXCL14 and high-grade gliomas(HGGs), we depicted the distribution of CXCL14 expression in different molecular subtypes. As shown in Figure 3G-H, CXCL14 was significantly downregulated in proneural subtype than other subtypes in both TCGA and CGGA dataset. ROC curves for CXCL14 expression and proneural subtype in all HGGs were performed and under curve areas are 89.08% and 87.17% in TCGA and CGGA dataset, respectively (Figure 3I-J). We analyzed the possible relationship between CXCL14 expression and prognosis in glioma patients using the TCGA and CGGA dataset. Up-regulation of CXCL14 expression correlated with a worse survival rate (Figure 3K; p = 0.02809) in TCGA patients with GBMs.Meanwhile, high expression of CXCL14 also showed decreased survival rate (Figure 3L; p < 0.0001) in patients with HGGs from CGGA dataset. Univariate and multivariate Cox proportional hazards analysis incorporating the age and molecular subtype factors were also performed (Supplementary Table 2 and 3).

CXCL14 related biological process

Genes significantly correlated with PD-1 expression (Pearson |R| > 0.4) were screened out in TCGA (429 genes) and CGGA (821 genes) datasets. We further analyzed these genes by Functional Annotation Tool online (DAVID, https://david.ncifcrf.gov/). These genes were mainly involved in angiogenesis, extracellular and immune response related biological process (Figure 4).

Figure 4. Gene ontology analysis for CXCL14 in gliomas.

A, Gene ontology analysis for CXCL14 in TCGA dataset. B, Gene ontology analysis for CXCL14 in CGGA dataset.

CXCL14 expression correlates with distinct immune cell populations infiltration and immune system-related metagenes

As indicated above, CXCL14 plays an important role in the immune response in glioma. We used the Microenvironment Cell Populations-counter method described by Etienne Becht et all²⁵ to evaluate the association between CXCL14 and immune cell populations from transcriptomic data. We detected a strong correlation of CXCL14 with T cells, Monocytic lineage, Neutrophilis and Fibroblasts cells (Figure 5A and B).To specifically analyze the role of CXCL14 in the glioma immune response, we evaluated its association with seven immune related metagenes.²⁶ As shown in Figure 5C and D, CXCL14 expression was positively associated with 5 of these metagenes expect IgG and Interferon.

Figure 5. Lower CXCL14 expression correlates with distinct immune cell populations infiltration and immune system-related metagenes.

A-B, Association between CXCL14 expression and immune cell populations in TCGA (A) and CGGA (B) datasets. Positive correlations are shown in blue and negative correlations in red color. Color depth and the size of the circle are proportional to the correlation coefficients. The numbers shown inside boxes are the correlation coefficients r. The corrplot display is symmetric about a diagonal line. For example, in Figure 5A, the r = -0.23(in 1^st column and 4^th row), corresponding to the red circle in 4^th column and 1^st row is the correlation coefficient of CXCL14 level with cytotoxic lymphocytes population in glioma. C-D, CXCL14-related immune response in TCGA (C) and CGGA (D) datasets.

Reduced mir-17-5p expression predicts poor survival in GBM patients

To further explore whether miR-17-5p has consistently inhibitory effect on tumor growth in vivo, we intracranially injected GBM cells stably expressing miR-17-5p or miR-NC into nude mice. As shown in Figure 6A, tumor growth was markedly slower in mice injected with primary GBM N3 cells overexpressing miR-17-5p than those overexpressing miR-NC. In agreement with the tumor growth, the mice in miR-17-5p overexpression mice showed significantly better survival than those from the miR-NC overexpression group (Figure 6B). IHC analysis also demonstrated that miR-17-5p-overexpressing tumours showed decreased CXCL14 staining (Figure 6C, right), a lower proliferation index as indicated by Ki-67 staining, and significantly lower levels of angiogenesis as indicated by CD31 expression (Figure 6C, left). To explore the association of miR-17-5p expression level with prognosis of GBM patients, Kaplan-Meier survival curve analysis with a log-rank test of GBM patients from TCGA dataset was performed. As shown in Figure 6D, GBM patients with lower miR-17-5p expression had a worse overall survival time in comparison with those with higher miR-17-5p expression (P < 0.0001). Meanwhile, the same result was also confirmed in HGGs patients form CGGA (P = 0.04969) validation sets (Figure 6E). Univariate and multivariate Cox proportional hazards analysis were used to further reveal the prognostic value of miR-17-5p in GBM patients. As shown in Table 1–2, lower miR-17-5p expression was a risk factor for GBM patients. After incorporating the age and molecular subtype factors, the analysis indicated that miR-17-5p value was not an independent prognostic marker for the overall survival of GBM patients.

Figure 6. Reduced miR-17-5p expression predicts poor survival in GBM patients.

A, Representative bioluminescence images of mice bearing intracranial pCDH- or pCDH miR-17-5p-transfected N3 cells on the days indicated (n = 10 each group). B, Kaplan–Meier survival curve of mice injected with bearing pCDH- or pCDH miR-17-5p-transfected N3 cells. C, Immunohistochemical analysis of paired sections of Ki-67, CD31 and CXCL14 expression in intracranial tumours. Scale bar = 100 μm (lower). D-E, Survival analysis for miR-17-5p in TCGA (D) and CGGA (E) datasets. F-G, Univariate and multivariate analysis of gliomas in TCGA (F) and CGGA (G) datasets.

Discussion

MicroRNAs, a subset of small regulatory RNAs, have been elucidated to play important roles in a wide variety of oncogenic activities, such as tumor growth, invasion, migration, immune response and angiogenesis.²⁷ Aberrant expression of miRNAs has been observed in various kinds of cancers, such as glioma and its aggressive glioblastoma subtype.²⁸ Tumor tissues is comprised of diverse mixture of tumor and nontumor cells within their microenvironments, and the nontumor cells, such as immune cells and stromal cells dilute the tumor purity. Tumor purity is routinely and previously analyzed by pathologists via visual evaluation, which could be influenced by the sensitivity of histopathology, interobserver bias, and variation in accuracy. Here, we employed ESTIMATE algorithm in both microarray- and RNAseq-based transcriptome profiles to systematically predict glioma purity of patients from TCGA and CGGA dataset. The main advantage of our study was the utilization of large sample-sized glioma cohorts and systematic analysis in multidimensional conditions to depict the miRNAs expression pattern within GBM environments under various purities. In this study, we identified an uncharacterized miR-17-5p to be lowly expressed in GBM patients with lower tumor purity and functionally required for the malignant phenotype. miR-17-5p overexpression inhibited proliferation of GBM cells by complementarily binding 3 ′ -UTR of CXCL14. Moreover, miR-17-5p expression predicted significantly better survival for HGG patients.Meanwhile, miR-17-5p expression is upregulated in proneural subtypes and may serve as a promising indicator of proneural subtype.

Chemokines are demonstrated as critical modulators of cell migration in development, homeostasis, and immune surveillance. Some chemokines are considered pro-inflammatory and can be induced during an immune response to recruit immune cells to a site of inflammation. In non-cancerous tissue, CXCL14 is predominantly expressed in the epithelium while it indicates an altered expression pattern during tumorigenesis. In some tumors, such as those of the prostate and pancreas, CXCL14 showed increased expression.^29-31 In contrast, CXCL14 expression is frequently decreased in cancer cells of the head and neck, breast, cervix, and kidneys.³² However, the role of CXCL14 in tumorigenesis of glioma are not clearly explained. Here, we found that CXCL14 gene is a target of miR-17-5p in glioma. CXCL14 oncogene is experimentally validated as the novel target of miR-17-5p not only in vitro, but also in vivo. Within the tumor environment, CXCL14-overexpressing fibroblasts could accelerate tumor growth and vascularization by inducing upregulation of FGF basic and VEGF-A, -B, and – C.³³ CXCL14 in tumors also could modulate tumor proliferation by attracting fewer dendritic cells, thereby regulating the initiation of an anti-tumor immune response in head and neck squamous cell carcinoma.³⁴ Here, we found that CXCL14 expression was downregulated in higher-purity gliomas. GO term analysis indicated that CXCL14 expression was involved in some extracellular, angiogenesis and immune response related biological process within glioma environment. GBM had been classified into four molecular subtypes: proneural, neural, classical, and mesenchymal.¹³ The proneural subtype is characterized of marker expression of neuroblasts or neurons, low rate of proliferation, and associated with better outcome.^35,36 Here, we found proneural subtype expressed significantly lowest level of CXCL14 while mesenchymal subtype tended to express highest level of CXCL14 among four subtypes. These results were consistent with previous study that mesenchymal glioblastomas are characterized of anti-tumor inflammatory responses than other subtypes. Moreover, CXCL14 expression predicted significantly worse survival for HGG patients.

The mechanisms of CXCL14 expression regulation are not well elucidated. Through screening transcriptomic data, we firstly determined the association between CXCL14 and immune cell populations. We found that recruitment and function of immune cells played important roles in this process. As reported previously, CXCL14 expression was involved in the aspects of neutrophil chemotaxis, cancer-associated fibroblasts progression and activation, and regulatory T Cells activation,^37-40 which have attracted the interest of targeting these cells in cancer therapy. In the present study, a strong correlation of CXCL14 expression with T cells, Monocytic lineage, Neutrophils and Fibroblasts cells within glioma environment in both TCGA and CGGA datasets were demonstrated.Futhermore, we also explored the association between CXCL14 expression and 7 immune related metagenes.²⁶ CXCL14 were positively associated with HCK, LCK, MHC-I, MHC-II and STAT1.These 5 metagenes positively correlating with CXCL14 were mainly involved in activities of T cells and antigen-presenting cells, which was in line with the analysis mentioned above.

In conclusion, our findings demonstrate that dysregulation of the miR-17-5p/CXCL14 signaling pathway in GBM promises to improve our knowledge of the biological basis of GBM development and tumor environment. It might also be a potential therapeutic target in immunotherapy of GBMs.

Patients and methods/materials and methods

Sample and data collection

For CGGA transcriptome data(http://www.cgga.org.cn), we have enrolled 301 samples ranging from WHO grade II to grade IV and generated via Agilent Whole Human Genome Array platform. For CGGA miRNA microarray dataset (http://www.cgga.org.cn), we have collected 198 samples ranging from WHO grade II to grade IV and generated via Infinium HumanMethylation27 BeadChip platform. Dataset of TCGA was downloaded from TCGA portal (https://tcga-data.nci.nih.gov/docs/publications/tcga/). For CGGA project, written informed consent was obtained from each participant. This study protocol was approved by the Ethics Committee of Nanjing Medical University, Nanjing, China. All the involvement with human subjects in the present study complies with the Declaration of Helsinki.

IDH mutation detection

For the subjects of CGGA, IDH1/2 mutation status was determined by the method of pyrosequencing as described previously. Meanwhile, IDH1/2 mutation data for the TCGA subjects were analyzed by whole exon sequencing or pyrosequencing.

Cell culture

Primary human GBM N3 cell (N3) was a gift from Beijing Neurosurgical Institute, Capital Medical University.N3 cells and U251 cells that were purchased from ATCC were incubated in Dulbecco’s modified Eagle’s medium (DMEM) supplemented with 10% fetal bovine serum (FBS). These cell lines were authenticated in January 2014 through the STR genotyping method (Genechem, Shanghai, China).

Luciferase assays and western blotting

U251 and N3 cells were seeded into a 24-well plate. After cultured overnight, cells were co-transfected with the wild-type or mutated pMIR-reporter luciferase vector plasmid, and transfected with miR-17-5p or miR-NC mimics. Luciferase assays were done 48 h after transfection via the Dual Luciferase Reporter Assay System (Promega, WI, USA). Western blot analysis were performed as previously described⁴¹ and anti-bodies against CXCL14(Abcam,1:1000) and GAPDH (1: 5,000, Santa Cruz) were used.

Xenograft tumour assay

A lentiviral pCDH empty vector(pCDH) and miR-17-5p expressing vector (pCDH miR-17-5p) were purchased from GeneChem Co. Ltd (Shanghai, China). Male BALB/c nude mice (6 weeks old) were used in this study. All procedures involving experimental animals were performed according to a protocol approved by Nanjing Medical University. For establishing intracranial xenografts, 5.0 × 10⁴ N3 cells stably expressing luciferase reporter were stereotactically implanted as previously described.⁵ The bioluminescence imaging was employed to determine the Fluc activity for mice. IHC detecting CXCL14, Ki-67 or CD31 (Abcam) expression was determined on tumour tissue from nude mouse xenografts through methods described previously.^5,42

Statistical analysis

In the present study, statistical analysis was mainly done via R language (https://www.r-project.org/) with several publicly available packages. Figures were mainly generated by R packages: ComplexHeatmap package was employed to depict the expression pattern shown in Figure 1A and Figure 4; pROC package was used to draw the receiver operating characteristic (ROC) curve for evaluating the quality or performance of diagnostic tests in Figure 2E-F and 3I-J; corrplot package was applied to determine correlation coefficient r of CXCL14 expression with various immune cell population (Figure 5A-B) or inflammatory metagenes(Figure C-D); Univariates and multivariates cox proportional hazard model were analyzed using coxph function provided by survival package. Meanwhile, The Kaplan–Meier survival curve was depicted to estimate survival distributions. P-values of < 0.05 (two-sided) were considered significant (single asterisks in the Figs).

Disclosure of interest

The authors report no conflict of interest

Acknowledgments

This work was supported by a grant from the National Natural Science Foundation of China (No. 81402056 and 81472362), National High Technology Research, Development Program of China (863) (No.2012AA02A508) and the National Key Research and Development Plan (No. 2016YFC0902500).

Supplementary data

Supplementary data can be accessed here.

Additional information

Funding

This work was supported by the National Natural Science Foundation of China [81402056 and 81472362];National Key Research and Development Plan [2016YFC0902500];National High Technology Research, Development Program of China (863) [2012AA02A508];