Online ISSN: 1945-4589
Research Paper Volume 15, Issue 6 pp 2321—2346

Downregulated annexin A1 expression correlates with poor prognosis, metastasis, and immunosuppressive microenvironment in Ewing’s sarcoma

Xiangwen Shi1,2, *, , Yipeng Wu1,2, *, , Linmeng Tang3, *, , Haonan Ni1, , Yongqing Xu2, &, ,

  • 1 Kunming Medical University, Kunming, China
  • 2 Laboratory of Yunnan Traumatology and Orthopedics Clinical Medical Center, Yunnan Orthopedics and Sports Rehabilitation Clinical Medical Research Center, Department of Orthopedic Surgery, 920th Hospital of Joint Logistics Support Force of PLA, Kunming, China
  • 3 Bone and Joint Imaging Center, Department of Medical Imaging and Radiology, The First Affiliated Hospital of Hebei North University, Zhangjiakou, China
* Equal contribution

Received: January 15, 2023       Accepted: March 15, 2023       Published: March 28, 2023      

https://doi.org/10.18632/aging.204615
How to Cite

Copyright: © 2023 Shi et al. This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Abstract

Objective: Ewing’s sarcoma (ES) is a common bone malignancy in children and adolescents that severely affects the prognosis of patients. The aim of this study was to identify novel biomarkers and potential therapeutic targets for ES.

Methods: Highly prognosis-related hub genes were identified by independent prognostic analysis in the GSE17679 dataset. We then performed survival analysis, Cox regression analysis and clinical correlation analysis on the key gene and validated them with the GSE63157, GSE45544 and GSE73166 datasets. Differentially expressed genes (DEGs) were screened based on the high and low expression of key gene, Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses, Gene Set Enrichment Analysis (GSEA), and Gene Set Variation Analysis (GSVA) were performed to explore the underlying mechanisms of ES, and significant module genes were established based on protein-protein interaction (PPI) networks. Furthermore, the correlations between module genes and the immune microenvironment were analyzed and the correlations between the key gene and immune infiltration levels in sarcoma were investigated using TIMER and TISIDB. Finally, the expression levels of these key genes in ES cell lines (RD-ES and A673 cells) were further validated by real-time quantitative PCR (RT-qPCR). CCK-8 and EdU assays were performed to assess the effect of ANXA1 knockdown on RD-ES cell proliferation.

Results: ANXA1 was identified as a key gene for ES prognosis. The overall survival (OS) time of patients with low ANXA1 expression was shorter, and the expression level of ANXA1 in the metastatic group was significantly lower than that in the primary group (P<0.01). Additionally, the abundance of 12 immune cells in the ANXA1 low-expression group was significantly lower than that in the high-expression group (all P<0.05), which may be related to the inhibition of the immune microenvironment. A PPI network was constructed based on 96 DEGs to further identify the five ANXA1-related module genes (COL1A2, MMP9, VIM, S100A11 and S100A4). The expression levels of ANXA1, COL1A2, MMP9, VIM, S100A11 and S100A4 were significantly different between ES cell lines and mesenchymal stem cells after validation in two ES cell lines (all P<0.01). Among these genes, ANXA1, COL1A2, MMP9, VIM and S100A4 were significantly associated with the prognosis of ES patients (all P<0.05). Importantly, ANXA1 knockdown significantly promoted the proliferation of RD-ES cells, which may explain the susceptibility to ES metastasis in the ANXA1 low-expression group.

Conclusions: ANXA1 may serve as an independent prognostic biomarker for ES patients and is associated with metastasis and the immunosuppressive microenvironment in ES, which needs to be validated in further studies.

Introduction

Ewing’s sarcoma (ES) is the second most frequent and highly invasive bone tumor, accounting for 15% of all bone tumors in children and adolescents [1]. According to early data from the American Cancer Center, there are approximately three cases of ES per million people each year [2]. The incidence of ES is higher in white populations, and it is more common in males for unknown reasons [3, 4]. Despite being less than 1% of all human malignant tumors [5], sarcomas are all invasive, and the most common metastatic sites are the lung, bone, and bone marrow [6].

The standard treatment for ES includes a combination of surgery, local radiotherapy, and drug chemotherapy [7, 8]. After the introduction of chemotherapy, the survival rate of local ES patients increased from 10% to approximately 75%, but the metastasis of ES is still not optimistic [9, 10]. It has been reported that approximately 20-25% of patients with ES have metastases at the time of diagnosis, with a 5-year survival rate of less than 30% [11, 12]. In view of the above factors, new treatments are urgently needed to improve the prognosis of ES patients. In contrast to other sarcoma types, specific chromosome translocations such as EWS-TL11, are necessary conditions for Ewing’s sarcoma [13, 14]. Therefore, the pathological process controlled by the fusion protein determines the importance of molecular targeted therapy for Ewing’s sarcoma. Several previous studies evaluated indicators affecting ES survival and metastasis based on the SEER database and additional multicenter cohort data, and further developed nomogram models to predict the incidence of ES [15, 16]. Recently, the development of public databases and genomics has made it critical to seek potential prognostic biomarkers for ES, which would improve the prognosis of patients and guide their treatment.

In this study, we analyzed four datasets, GSE17679, GSE63157, GSE45544 and GSE73166, from the Gene Expression Omnibus (GEO). High prognosis-related genes in the GSE17679 dataset were screened by independent prognostic analysis, and survival analysis, Cox regression analysis and clinical correlation analysis were performed based on these genes. 88 ES samples from the GSE17679 dataset were used as the experimental group, and 107 ES samples from the GSE63157, GSE45544 and GSE73166 datasets were used as the validation group to validate the results of the survival analysis, independent prognostic analysis and clinical correlation analysis. Differential analysis, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis were performed with median expression values of hub gene. Additionally, a protein-protein interaction (PPI) network was constructed to identify the module genes associated with hub gene to further explore the relationship between these module genes and the ES immune microenvironment. Finally, the expression of module genes was validated in ES cell lines (RD-ES and A673 cells) and mesenchymal stem cells (MSCs), and survival analysis was conducted. Exploring the effect of ANXA1 knockdown on the proliferation of RD-ES cells by CCK8 and EdU assays.

Results

Screening of hub genes

First, a total of 64 genes were screened by Cox regression analysis of patients with ES based on OS time and survival status (P<0.001) (Supplementary Table 1). Next, 50 genes were further obtained from 88 samples screened by survival analysis according to OS time and survival status (P<0.001) (Supplementary Table 2). Finally, 374 genes were significantly associated with at least one clinical characteristic through clinical correlation analysis, of which 50 genes were significantly associated with two clinical features, and 324 genes were significantly associated with only one clinical feature (Supplementary Table 3). Based on survival (P<0.001), independent prognostic (P<0.001), and clinical correlation analyses (SigNum=2), ANXA1 was determined to be the final hub gene.

Survival analysis of ANXA1 and validation

A total of 88 samples in the training set were divided into a high-expression group and a low-expression group according to the median expression of ANXA1, and a survival curve was drawn. The results suggested that the OS time of patients with low-expression was significantly lower than that of patients with high-expression (P < 0.001), and the five-year survival rates of the ANXA1 low and high-expression groups were approximately 30% and 70%, respectively (Figure 1A). To validate the results of ANXA1 survival analysis, 85 samples from the GSE63157 dataset were divided into high and low-expression groups. According to the results, patients in the low-expression group had a shorter OS time than those in the high expression group (P=0.001), and the five-year survival rates of the ANXA1 low and high-expression groups were approximately 50% and 75%, respectively, which supported the conclusion of the experimental group (Figure 1B). To further analyze the prognostic value of ANXA1 in patients with different types of ES, survival analyses were performed in primary ES, recurrent ES, and metastatic ES patients from the GSE17679 dataset. The results showed that the OS time of the ANXA1 low expression group was significantly lower than that of the ANXA1 high expression group in primary (P<0.001) and metastatic ES patients (P=0.028), while there was no significant difference in recurrent ES patients (Supplementary Figure 1).

ANXA1-based survival curves, univariate and multivariate Cox regression analysis. (A) The results of Kaplan–Meier survival analysis based on the GSE17679 dataset suggested that patients in the ANXA1 low expression group had worse OS times than those in the low expression group (PB) The results of Kaplan–Meier survival analysis based on the GSE63157 dataset suggested that patients in the ANXA1 low expression group also had worse OS time, which was consistent with the training set (P=0.001). (C) Univariate Cox regression analysis for ANXA1 in the GSE17679 dataset (PD) Multivariate Cox regression analysis for ANXA1 in the GSE17679 dataset (PE) Univariate Cox regression analysis for ANXA1 in the GSE63157 dataset (PF) Multivariate Cox regression analysis for ANXA1 in the GSE63157 dataset (P

Figure 1. ANXA1-based survival curves, univariate and multivariate Cox regression analysis. (A) The results of Kaplan–Meier survival analysis based on the GSE17679 dataset suggested that patients in the ANXA1 low expression group had worse OS times than those in the low expression group (P<0.001). (B) The results of Kaplan–Meier survival analysis based on the GSE63157 dataset suggested that patients in the ANXA1 low expression group also had worse OS time, which was consistent with the training set (P=0.001). (C) Univariate Cox regression analysis for ANXA1 in the GSE17679 dataset (P<0.001). (D) Multivariate Cox regression analysis for ANXA1 in the GSE17679 dataset (P<0.001). (E) Univariate Cox regression analysis for ANXA1 in the GSE63157 dataset (P<0.001). (F) Multivariate Cox regression analysis for ANXA1 in the GSE63157 dataset (P<0.001).

Univariate and multivariate Cox regression analysis of ANXA1

The results of univariate analysis based on sex, age and tumor status suggested that ANXA1 was significantly associated with the survival and status of ES patients (P<0.001) and that the low HR value (HR<1) of ANXA1 indicated that it may be a low risk factor for prognosis in ES patients (Figure 1C). Multifactorial regression analysis indicated that ANXA1 was a low-risk factor for prognosis in ES patients (P<0.001) (Figure 1D). Univariate and multifactorial regression analyses from the validation group similarly demonstrated that ANXA1 could be an independent prognostic factor for predicting prognosis in ES patients (P<0.001) (Figure 1E, 1F).

The expression of ANXA1 in clinical subgroups

After dividing the clinical characteristics into specific subgroups, we analyzed the expression of ANXA1 in a clinical subgroup from the GSE17679 dataset. In the subgroup of age, the median expression of ANXA1 in the ≤20-year-old group was higher than that in the >20-year-old group (P=0.018) (Figure 2A). In the gender subgroup, there was no significant difference in ANXA1 expression levels between male and female patients (P>0.05) (Figure 2B). In the subgroup of tumor status, the median expression of ANXA1 in the metastasis ES group was significantly lower than that in the primary ES group (P=0.0072), but there was no significant difference between the primary group and recurrence group (P>0.05) (Figure 2C). Additionally, 22 ES patients from combined the GSE45544 and GSE73166 datasets were used to verify the age, gender and tumor status of the subgroup, and the results showed that the median expression of ANXA1 in the metastasis ES group was significantly lower than that in the primary ES group (P=0.00034) (Figure 2D).

ANXA1 expression in different subgroups based on the GSE17679, GSE45544 and GSE73166 datasets. (A) In different age groups, the expression levels of the ANXA1 gene were higher in the ≤ 20-year-old group than in the > 20-year-old group (P = 0.018). (B) In different gender groups, the expression levels of ANXA1 gene were not significantly different (P = 0.067). (C) In different type groups, the expression levels of the ANXA1 gene were higher in the primary group than in the metastasis group (P = 0.0072), and there was no significant difference between the remaining groups. (D) 22 ES samples from the GSE45544 and GSE73166 datasets validated that the expression levels of the ANXA1 gene were higher in the primary group than in the metastasis group.

Figure 2. ANXA1 expression in different subgroups based on the GSE17679, GSE45544 and GSE73166 datasets. (A) In different age groups, the expression levels of the ANXA1 gene were higher in the ≤ 20-year-old group than in the > 20-year-old group (P = 0.018). (B) In different gender groups, the expression levels of ANXA1 gene were not significantly different (P = 0.067). (C) In different type groups, the expression levels of the ANXA1 gene were higher in the primary group than in the metastasis group (P = 0.0072), and there was no significant difference between the remaining groups. (D) 22 ES samples from the GSE45544 and GSE73166 datasets validated that the expression levels of the ANXA1 gene were higher in the primary group than in the metastasis group.

Identification and enrichment analysis of DEGs

A total of 95 DEGs were identified, of which 81 were upregulated genes and 14 were downregulated genes. The visual volcano plot and heatmap showed that the high-expression genes and the low-expression genes could be clearly distinguished according to the median expression of ANXA1 (Figure 3A, 3B). Further analysis of the correlation between ANXA1 and DEGs showed that 81 genes were positively correlated with ANXA1, and 14 genes were negatively correlated with ANXA1 (Figure 3C). In GO enrichment analysis (Figure 4A, 4B), DEGs were mainly involved in biological processes (BPs) such as the formation of extracellular matrix organization (16 genes), formation of extracellular structure organization (16 genes), ossification (17 genes), formation of external encapsulating structure organization (16 genes), and formation of collagen fibril organization (8 genes). In cellular components (CCs), DEGs were mainly involved in collagen-containing extracellular matrix (36 genes), collagen trimer (10 genes), fibrillar collagen trimer (4 genes), banded collagen fibril (4 genes), and basement membrane (7 genes). In terms of molecular functions (MFs), DEGs were significantly involved in extracellular matrix structural constituent (26 genes), collagen binding (10 genes), integrin binding (11 genes), extracellular matrix structural constituent conferring tensile strength (6 genes), and extracellular matrix structural constituent conferring compression resistance (5 genes). Additionally, the enrichment results of the KEGG pathway analysis showed that DEGs were related to phagosomes (8 genes), Staphylococcus aureus infection (6 genes), and protein digestion and absorption (6 genes) (Figure 4C). The top 5 pathways are listed in Table 1.

Visualization of DEGs and correlation analysis of ANXA1. (A) Volcano plot of DEGs. Red dots indicate upregulated genes, green dots indicate downregulated genes, and black dots indicate genes with insignificant differences. (B) Heatmap of DEGs. Red indicates high expression, blue indicates low expression, and white indicates moderate expression. (C) Correlation coefficient heatmap of DEGs. Red represents positive correlation and green represents negative correlation.

Figure 3. Visualization of DEGs and correlation analysis of ANXA1. (A) Volcano plot of DEGs. Red dots indicate upregulated genes, green dots indicate downregulated genes, and black dots indicate genes with insignificant differences. (B) Heatmap of DEGs. Red indicates high expression, blue indicates low expression, and white indicates moderate expression. (C) Correlation coefficient heatmap of DEGs. Red represents positive correlation and green represents negative correlation.

GO functional annotation and KEGG pathway analysis of DEGs. (A) The top 5 biological processes (BP), cellular components (CC), and molecular functions (MF) of DEGs. (B) Enrichment circle diagram of DEGs. (C) KEGG pathway analysis of DEGs.

Figure 4. GO functional annotation and KEGG pathway analysis of DEGs. (A) The top 5 biological processes (BP), cellular components (CC), and molecular functions (MF) of DEGs. (B) Enrichment circle diagram of DEGs. (C) KEGG pathway analysis of DEGs.

Table 1. TOP 5 KEGG pathway of differentially expressed genes (DEGs).

IDDescriptionGene countP valueGenes
hsa03040Spliceosome232.08×10-9HNRNPM/SART1/DDX42/HNRNPA3/ISY1/SF3B2/SNRPD3/PUF60/SRSF4/PLRG1/HNRNPA1/CCDC12/SYF2/RNU4-1/WBP11/THOC3/RNU5A-1/LSM3/SNRPA1/SRSF5/RNU2-1/SNRPA/RNU4-2
hsa05014Amyotrophic lateral sclerosis325.63×10-6ATXN2/PSMC5/NUP214/NDUFS6/ATF4/PPP3R1/COX5B/HNRNPA3/PSMC1/COX4I1/ATP5MC1/NDUFC2/HNRNPA1/PSMD6/NDUFA12/NDUFA4/PFN1/NDUFA2/PSMA5/RAB1A/NDUFA5/GPX7/PSMA4/COX8A/PSMB3/ATP5F1A/PSMA2/NUP107/COX7A2/SIGMAR1/ACTR1B/COX7B
hsa05014Prion disease265.04×10-6PSMC5/NDUFS6/ATF4/PPP3R1/COX5B/PSMC1/COX4I1/ATP5MC1/NDUFC2/PSMD6/NDUFA12/NDUFA4/NDUFA2/PSMA5/NDUFA5/CREB3L2/PSMA4/COX8A/PSMB3/ATP5F1A/PSMA2/CSNK2A2/COX7A2/COX7B/CYBA/NCF2
hsa05014Oxidative phosphorylation172.25×10-6NDUFS6/COX5B/COX4I1/ATP5MC1/NDUFC2/NDUFA12/NDUFA4/NDUFA2/NDUFA5/COX8A/ATP5F1A/ATP6V0E1/COX7A2/ATP6V0B/ATP6V1F/COX7B/ATP6V1C1
hsa05012Parkinson disease242.91×10-5PSMC5/GNAI3/NDUFS6/ATF4/COX5B/PSMC1/COX4I1/ATP5MC1/NDUFC2/SLC39A7/PSMD6/NDUFA12/NDUFA4/NDUFA2/PSMA5/NDUFA5/PSMA4/COX8A/PSMB3/ATP5F1A/PSMA2/COX7A2/TXN/COX7B

To further identify the expression differences in the involved KEGG pathways, GSEA and GSVA enrichment analyses revealed that ANXA1 was involved in the cell adhesion molecules cams, cell cycle, cytokine-cytokine receptor interaction, DNA replication, ECM receptor interaction and focal adhesion pathways (Figure 5A, 5B).

Single-gene gene set enrichment analysis (ssGSEA) and gene set variation analysis (ssGSVA) pathway analysis of ANXA1. (A) The top 6 underlying related KEGG enrichment pathways of ANXA1 through single-gene GSEA. (B) Results of single-gene GSVA analysis involving KEGG pathways.

Figure 5. Single-gene gene set enrichment analysis (ssGSEA) and gene set variation analysis (ssGSVA) pathway analysis of ANXA1. (A) The top 6 underlying related KEGG enrichment pathways of ANXA1 through single-gene GSEA. (B) Results of single-gene GSVA analysis involving KEGG pathways.

PPI network construction and modular analysis

The DEGs were introduced into STRING to draw the PPI network plot. After Cytoscape processing, the module with the highest MCODE score (score=13.875) was constructed (Figure 6A), which included 71 nodes and 304 edges, as well as 69 upregulated genes and 1 downregulated gene. Next, a new module was constructed using ANXA1 as the node, which included COL1A2, MMP9, VIM, S100A11 and S100A4, all of which were upregulated genes (Figure 6B).

Protein interaction network diagram and survival curves of module genes. (A) The module with the highest MCODE score. (B) The module associated with AXNA1. Red indicates upregulated genes, green indicates downregulated genes, and hub gene is indicated in blue. The expressions of ANXA1-related module genes used to construct the prognosis of ES. (C) S100A11; (D) MMP9; (E) COL1A2; (F) VIM; (G) S100A4. P values

Figure 6. Protein interaction network diagram and survival curves of module genes. (A) The module with the highest MCODE score. (B) The module associated with AXNA1. Red indicates upregulated genes, green indicates downregulated genes, and hub gene is indicated in blue. The expressions of ANXA1-related module genes used to construct the prognosis of ES. (C) S100A11; (D) MMP9; (E) COL1A2; (F) VIM; (G) S100A4. P values <0.05 indicated that differences were significantly significant.

Survival analysis of DEGs between the two groups

The five modular genes were differentially expressed in RD-ES cells or A673 cells compared to MSCs, and we further established survival curves to explore the effects of these genes on the prognosis of ES. The results revealed that except for S100A11 (Figure 6C), the expression of COL1A2, MMP9, VIM and S100A4 significantly affected the overall survival time of patients with ES, in which the OS time of the high expression group of MMP9 and COL1A2 was significantly lower than that of the low expression group (Figure 6D, 6E), and the OS time of the low expression group of VIM and S100A4 was significantly lower than that of the high expression group (Figure 6F, 6G).

Immune infiltration analysis of ANXA1 and five module genes

To explore the relationship between the 23 immune cells and the hub genes, a correlation heatmap was constructed (Figure 7A) and it suggested that 5 target genes other than the VIM gene had a strong correlation with 23 immune cells. Specifically, ANXA1, COL1A2, MMP9, S100A4 and S100A11 high-expression groups were associated with high infiltration of activated dendritic cells, gamma delta T cells, immature B cells, immature dendritic cells, MDSCs, macrophages, natural killer cells, regulatory T cells, T follicular helper cells, and type 1 T helper cells (Figure 7B–7G).

Immune infiltration results of the 6 hub genes. (A) Heatmap of the correlation between immunocyte abundance and 6 hub genes. Box diagram of immunocyte abundance in the ANXA1 (B), MMP9 (C), COL1A2 (D), S100A4 (E), VIM (F) and S100A11 (G) high- and low-expression groups.

Figure 7. Immune infiltration results of the 6 hub genes. (A) Heatmap of the correlation between immunocyte abundance and 6 hub genes. Box diagram of immunocyte abundance in the ANXA1 (B), MMP9 (C), COL1A2 (D), S100A4 (E), VIM (F) and S100A11 (G) high- and low-expression groups.

Previous studies have shown that tumor tissues contain tumor cells and lymphocytes, and TILs can predict the survival status of tumor patients [17, 18]. To further investigate the correlation between the expression levels of ANXA1 and TILs, we explored the correlation between the expression levels of ANXA1 and the infiltrating abundance of TILs in sarcoma patients through the TIMER database. In sarcoma, the expression levels of ANXA1 were positively correlated with those of MMP9 (r=0.204, P=9.43e-04), COL1A2 (r=0.349, p=8.87e-09), S100A4 (r=0.574, P<0.001), VIM (r=0.48, P<0.001) and S100A11 (r=0.638, P<0.001) (Figure 8A).

ANXA1 expression is correlated with 5 module genes and the level of immune infiltration in sarcoma. (A) ANXA1 expression is correlated with 5 module genes in sarcoma. (B) ANXA1 expression is correlated with the level of immune infiltration in sarcoma. (C) Kaplan–Meier plots of immune infiltration and ANXA1 expression levels in sarcoma.

Figure 8. ANXA1 expression is correlated with 5 module genes and the level of immune infiltration in sarcoma. (A) ANXA1 expression is correlated with 5 module genes in sarcoma. (B) ANXA1 expression is correlated with the level of immune infiltration in sarcoma. (C) Kaplan–Meier plots of immune infiltration and ANXA1 expression levels in sarcoma.

The expression levels of ANXA1 were negatively correlated with tumor purity (r=-0.271, P=1.64e-05), but positively correlated with the levels of infiltrating B cells (r=0.04, P=5.4e-01), CD8+ T cells (r=0.175, P=6.64e-03), CD4+ T cells (r=0.256, P=6.15e-05), macrophages (r=0. 469, P=2.97e-14), neutrophils (r=0.348, P=2.57e-08) and dendritic cells (r=0.306, P=1.25e-06), with the strongest correlation with macrophages (Figure 8B). We then analyzed the correlation between the expression levels of ANXA1 and immune cells and the prognosis of sarcoma patients, and the Kaplan–Meier plots suggested that the expression levels of CD4+ T cells (P=0.003) and neutrophils (P=0.01) were significantly associated with the prognosis of sarcoma patients (Figure 8C).

Next, we explored the correlation between the expression levels of ANXA1 and immune cells through the TISIDB database. Figure 9A shows a heatmap of correlation between the expression levels of ANXA1 and TILs in different types of tumors. In 263 sarcoma patients, the expression levels of ANXA1 were significantly correlated with the infiltration abundance of 26 TILs, with the strongest correlation with tcm_CD8, tcm_CD4, macrophage, treg, tfh, MDSC, tgd, NKT, act_DC and monocyte (Figure 9B and Supplementary Figure 2).

Correlation of ANXA1 expression with immune cells and immunoinhibitors in sarcoma. (A) Heatmap of the correlation between ANXA1 expression and TILs in sarcoma: red represents a positive correlation, blue represents a negative correlation. (B) ANXA1 expression was positively correlated with infiltrating abundance of tcm

Figure 9. Correlation of ANXA1 expression with immune cells and immunoinhibitors in sarcoma. (A) Heatmap of the correlation between ANXA1 expression and TILs in sarcoma: red represents a positive correlation, blue represents a negative correlation. (B) ANXA1 expression was positively correlated with infiltrating abundance of tcm_CD8, tcm_CD4, macrophage, treg, tfh, MDSC, tgd, NKT, act_DC and monocyte. (C) Heatmap of the correlation between ANXA1 expression and immunoinhibitors in sarcoma: red represents positive correlation, blue represents negative correlation. (D) ANXA1 expression was positively correlated with infiltrating abundance of HAVCR2, CSF1R, IL-10, LGALS9, IL10RB, PDCD1LG2, TGFB1, TGFBR1, TIGIT and CD96.

Finally, we revealed the correlation between the expression levels of ANXA1 and immunoinhibitors in sarcoma patients by a heatmap of correlation (Figure 9C), and the results suggested that the expression levels of ANXA1 were associated with 20 immunoinhibitors, among which HAVCR2, CSF1R, IL-10, LGALS9, IL10RB, PDCD1LG2, TGFB1, TGFBR1, TIGIT and CD96 showed strong positive correlations (Figure 9D and Supplementary Figure 3). These results suggested that ANXA1 may play an important role in the immunotherapy of sarcoma.

The expression levels of ANXA1 and five module genes

The RT-qPCR results showed that six hub genes, including ANXA1, were significantly differentially expressed in either RD-ES cells or A673 cells (Figure 10A). The relative expression levels of ANXA1 were significantly higher in RD-ES cells than in MSCs, while the relative expression levels of MMP9, COL1A2, S100A4, VIM and S100A11 were significantly lower than in MSCs. The relative expression levels of ANXA1, MMP9 and COL1A2 were significantly higher in A673 cells, while the relative expression levels of S100A4, VIM and S100A11 were significantly lower than those in MSCs (Figure 10B). Subsequent Western blot analysis verified that the protein levels of the six hub genes were all significantly differentially expressed in either RD-ES cells or A673 cells (Figure 10C). The expression of ANXA1 increased in ED-ES and A673 cells, while the expression of MMP9, COL1A2, S100A4, VIM and S100A11 decreased in ED-ES cells. It is worth mentioning that the protein levels of ANXA1, MMP9 and COL1A2 were significantly increased in A673 cells, while the protein levels of S100A4, VIM and S100A11 were significantly decreased in A673 cells, which was consistent with the results of RT-qPCR (Figure 10D).

The expression levels of ANXA1-related module genes in ES cell lines and MSCs by RT-qPCR and Western blot and functional validation of ANXA1 in RD-ES cell using CCK-8 and EdU assays. The relative mRNA levels of ANXA1, MMP9, COL1A2, S100A4, VIM and S100A11 in RD-ES (A) or A673 (B) cells compared to MSCs. (C) The protein levels of ANXA1-related module genes in ES cell lines and MSCs by Western blot. (D) Quantitative analysis of protein expression levels of module genes. (E) The efficiency of si-ANXA1 was detected by RT-qPCR and Western blot analysis. (F) ES-RD cell viability at 24 h, 48 h and 72 h after downregulation of ANXA1 was measured using a CCK-8 assay. (G) The effect of ANXA1 on ES-RD cell proliferation was assessed by an EdU assay. P values are shown as follows: * PPPP

Figure 10. The expression levels of ANXA1-related module genes in ES cell lines and MSCs by RT-qPCR and Western blot and functional validation of ANXA1 in RD-ES cell using CCK-8 and EdU assays. The relative mRNA levels of ANXA1, MMP9, COL1A2, S100A4, VIM and S100A11 in RD-ES (A) or A673 (B) cells compared to MSCs. (C) The protein levels of ANXA1-related module genes in ES cell lines and MSCs by Western blot. (D) Quantitative analysis of protein expression levels of module genes. (E) The efficiency of si-ANXA1 was detected by RT-qPCR and Western blot analysis. (F) ES-RD cell viability at 24 h, 48 h and 72 h after downregulation of ANXA1 was measured using a CCK-8 assay. (G) The effect of ANXA1 on ES-RD cell proliferation was assessed by an EdU assay. P values are shown as follows: * P<0.05, ** P<0.01, *** P<0.001, and **** P<0.0001.

ANXA1 inhibits the proliferation of RD-ES cells in vitro

Since ANXA1 downregulation was associated with worse prognosis and metastasis in ES patients, we further investigated the effect of ANXA1 downregulation on RD-ES cells in vitro. Figure 10E showed the efficiency test of si-ANXA1. CCK-8 and EdU assays revealed that, compared with NC group, the proliferation ability of cells in the HMGB1 knockdown group was increased (Figure 10F, 10G). The above results suggested that downregulated ANXA1 may affect the prognosis and metastasis of ES patients by promoting cell proliferation.

Discussion

As a common malignancy of bone and soft tissue, the prognosis is not promising for the pediatric and adolescent population, even after standardized chemotherapy regimens [19, 20]. With the development of genomics and high-throughput technologies providing a novel direction for the prognosis and treatment of ES, it is crucial to explore key biomarkers for the prognosis of ES patients. In this study, we first screened the ANXA1 gene using independent prognostic analysis and clinical correlation analysis from 88 ES samples obtained from the GSE17679 dataset. Then, survival analysis and clinical correlation analysis based on ANXA1 determined that low ANXA1 expression was associated with a shorter survival time, and this result was successfully validated by 85 ES samples from the GSE63157 dataset. In the clinical subgroups, the expression levels of ANXA1 were significantly lower in the metastatic group than in the primary group. We further analyzed 96 DEGs with ANXA1 expression as the median value and GO and KEGG pathway analyses suggested that DEGs were mainly associated with the extracellular matrix. Finally, we found that the expression of five ANXA1-related modular genes was significantly different between ES cell lines and MSCs by in vitro RT-qPCR, and four genes were significantly associated with the prognosis of ES patients.

ANXA1 is a phospholipid-binding protein located on chromosome 9q12-q21.2 and is expressed in many tissues and cells [21, 22]. Its C-terminus is composed of four annexin repeats with calcium binding sites, while the N-terminus contains important phosphorylation regulatory sites, which are unique to ANXA1 [23, 24]. In addition to mediating inflammatory responses, ANXA1 is involved in the development of multiple tumors, metastasis and drug resistance and may be used as a potential biomarker for tumor diagnosis, treatment, and prognosis [25, 26]. However, the role of ANXA1 in tumors remains uncertain, both as an antitumor factor and as a promoter of tumor progression [27]. To date, there is no literature on the relationship between ANXA1 expression and the prognosis of ES. It shows high expression levels in breast [28], colorectal [29] and prostate cancers [30]. This is consistent with the results of the present study, in which we found that ANXA1 showed significantly high expression in either RD-ES or A673 cells and a significant correlation with the OS time of patients with ES. In contrast, ANXA1 shows low expression levels in thyroid cancer [31] and nasopharyngeal carcinoma [32]. A recent study showed that ANXA1 expression was upregulated in patients with small cell lung cancer bone metastases and promoted bone metastasis of small cell lung cancer cells in mice, suggesting that ANXA1 may be a potential biomarker for lung cancer bone metastasis [33]. Interestingly, low expression levels of ANXA1 in this study were significantly associated with ES metastasis, suggesting that the downregulated ANXA1 gene was associated with the proliferation, invasion, and migration of ES cells.

The tumor microenvironment (TME) plays a key role in tumor and developmental processes and therapeutic susceptibility, and the typical immune profile of ES is a large infiltration of myeloid cells in the TME [34, 35]. Myeloid-derived suppressor cells (MDSCs) are involved in a variety of immunosuppressive responses as immature monocytes and granulocytes, and Zhang et al. [36] identified a subpopulation of MDSCs with significant immunosuppressive effects. Patients with ES showed significant depletion of CD4+ and CD8+ T cells in the peripheral blood compared to the healthy population [37]. Consistent with these results, ANXA1 downregulation in this study was significantly associated with reduced infiltration of macrophages, MDSCs and T cells, suggesting that ANXA1 is involved in the regulation of immune responses in the TME of ES patients. Previous studies showed that ANXA1 was involved in the proliferation and invasion of a variety of tumor cells [38, 39], and we also found that downregulated ANXA1 promoted the proliferation of ES cells, which may lead to progression and metastasis in ES patients.

Additionally, the COL1A2 gene was shown to be upregulated in ES tissues [40], and we also obtained consistent results in RD-ES cells. Zhang et al. [41] found that COL1A2 was involved in the development of ES and was significantly associated with the survival time of ES patients. Our survival curves also suggest that patients with low COL1A2 expression have shorter survival times. The S100A11 gene was found to be significantly upregulated in clear cell sarcoma of soft tissue [42]. In contrast, we found higher levels of S100A11 expression in ES-RD and A673 cells. High expression of MMP9 is associated with human tumor invasion or metastasis, and knockdown of MMP inhibits the migration of ES cells [43]. In this study, we verified the relationship between the expression of six key genes and ES and found that five of them were significantly associated with the prognosis of ES. At present, there are no literature reports of MMP9, S100A11, VIM and S100A4 as prognostic markers for ES, and their prognostic relationships need to be further investigated.

In this study, to screen for the hub gene with prognostic significance in ES patients, DEGs were screened using survival analysis and Cox regression analysis based on the OS time and status of ES patients. To reflect the clinical application value of hub genes, patients were further grouped according to their age, gender, and tumor status to explore the hub gene significantly associated with clinical characteristics and prognosis. Finally, differences in the gene and protein expression levels of the hub gene were examined in different ES cell lines and BMSCs using independent datasets to validate its prognostic value, which may provide theoretical support for prognosis and molecularly targeted therapy in clinical ES patients.

Admittedly, our study has some limitations that need to be considered. Firstly, both datasets used for analysis were small and a larger chip set should be utilized to validate the prognostic role of ANXA1 in ES. Secondly, the GEO dataset of the validation cohort lacked detailed clinical information about the patients, making it difficult to validate the relationship between ANXA1 and the clinical characteristics of ES patients. Finally, some basic experiments are also required to further investigate the effect of ANXA1 on the phenotypes of ES cells such as migration.

Conclusions

Taken together, we screened ANXA1 as an independent prognostic gene for ES based on the GEO database using multifactorial Cox regression analysis. ES patients with low expression of this gene had a shorter survival time and were significantly associated with ES metastasis and immunosuppressive microenvironment. In addition, the expression levels of ANXA1 were higher in ES cell lines. This study may facilitate the development of new prognostic biomarkers for targeted treatment of ES, although further experimental validations are needed.

Materials and Methods

Collection and processing of GEO data

Data on ES patients and clinical characteristics were obtained through the Gene Expression Omnibus website (https://www.ncbi.nlm.nih.gov/geo/). Four microarray datasets, GSE17679, GSE63157, GSE45544 and GSE73166, were extracted and the platform files of both datasets were downloaded. 88 ES samples and 18 normal samples from the GSE17679 dataset were obtained from the GPL570 platform, 85 ES samples from the GSE63157 dataset were obtained from the GPL5175 platform, and a total of 22 ES samples from the GSE45544 and GSE73166 datasets were obtained from GPL6244. To extract the survival time and survival status of ES patients, non-ES samples and ES cell lines were excluded. The GSE17679 dataset was the experimental group and the GSE63157, GSE45544 and GSE73166 datasets were the validation group. The results of survival analysis were validated by the GSE63157 dataset, and the results of clinical correlation analysis were validated by the GSE45544 and GSE73166 datasets after eliminating the batch effect.

Screening of hub gene related to survival and prognosis

First, the data from the gene expression matrix were corrected using the “impute” package of R software (version 4.1.3) [44]. For data with large values of gene expression, log 2 processing is needed. Next, the survival time and survival status of patients were extracted, and the gene expression data and survival information were merged using the “Perl” program [45]. The merged data were subjected to Cox regression analysis to filter genes according to overall survival (OS) time and survival status with a P value < 0.001. After obtaining the genes, survival analysis was performed through the “survival” package to screen for survival-associated genes with a P value threshold of 0.001. Finally, the samples were divided into two subgroups based on clinical characteristics such as sex and age, and three subgroups based on tumor status. Further clinical correlation analysis of genes was performed to identify the hub gene with a significant number (SigNum) of differences.

Survival analysis of hub gene and validation

Based on the median expression of the hub gene, 88 ES samples were divided into two groups for survival analysis. Furthermore, survival analyses for patients with different types of ES were performed. The high expression group indicated that the expression of the hub gene was greater than the median, and the low expression group indicated that the expression was less than or equal to the median. Moreover, 85 ES samples from the GSE63157 dataset were used to verify the accuracy of the survival curve.

Univariate and multivariate Cox regression analysis of hub gene

The hub gene expression data, clinical characteristics, and survival data of 88 ES samples were integrated, and univariate and multivariate Cox regression analyses were conducted using the “survival” package to identify characteristics that independently guided ES prognosis and to validate the results with the GSE63157 dataset.

Clinical correlation analysis of hub gene

ES patients were divided into two subgroups by age using 20 years as the threshold, two subgroups by gender, and three subgroups by tumor status: primary, recurrent, and metastatic. The hub gene was analyzed for differences among the subgroups.

Identification and correlation analysis of DEGs

Based on the median expression level of the hub gene, ES samples were divided into high and low expression groups. DEGs were screened using the “limma” package, and the concentrated area of upregulated or downregulated genes was indicated by a volcano plot and heatmap. The screening criteria were an absolute value of |logFC| > 1 and an adjusted P value < 0.05. To explore the relationship between the hub gene and DEGs, a significant P value was obtained by correlation test. Further visualization was performed in the form of a heatmap and correlation coefficient heatmap using the “corrplot” package.

GO and KEGG enrichment analysis

Gene Oncology (GO) and Kyoto Encyclopedia of Genes and Genomes Pathway Enrichment (KEGG) enrichment analyses based on DEGs were conducted using the “ClusterProfile” and “enrichplot” packages [46, 47]. GO enrichment analysis focused on the cytological components (CC), molecular functions (MF), and biological processes (BP) of DEGs, while KEGG enrichment analysis focused on the pathways of DEGs. The filtering condition was set to a P value < 0.05.

Single-gene gene set enrichment analysis (GSEA) enrichment

To further analyze the related pathways and potential biological functions of the hub gene in osteomyelitis, we used the “enrichplot” and “clusterProfiler” packages to perform GSEA enrichment analysis for each signature gene, with two gene sets, “c5.go.symbols.gmt” and “c2.cp.kegg.symbols.gmt”, as the predefined sets. The top 6 pathways with significant enrichment were visualized and the screening threshold was set at a P value < 0.05.

Single-gene gene set variation analysis (GSVA) enrichment

GSVA is a method to estimate the variation in pathway activity in samples with an unsupervised way due to its stability and is often used in the data analysis of gene expression profiles [48]. The hub gene was analyzed by GSVA based on predefined sets “c5.go.symbols.gmt” and “c2.cp.kegg.symbols.gmt”. Firstly, the samples are scored and corrected. Then, the samples were divided into groups according to the expression of the target gene, and the difference in GSVA score between the high expression group and the low expression group was further analyzed. The screening conditions for significant differences were |t| >2, P value < 0.05.

Construction of the PPI network and module analysis

The protein-protein interaction network of the hub genes was constructed using STRING (http://string-db.org) [49]. Cluster analysis was further performed using the Molecular Complex Detection (MCODE) plugin in Cytoscape (version 3.9.1) [50] to identify key modules in the PPI network to construct subnetworks.

Immune infiltration analysis and survival analysis of module genes

We further used ssGSEA in the “gsva” package to combine the 23 immune gene datasets with “high-low discriminant analysis” and calculate the immune infiltration score for each ES sample [51]. The correlation between hub genes and immune cell infiltration was represented by a heatmap. Based on the median expression level of each gene, the samples were divided into low and high expression groups, and boxplots were used to observe whether there was a significant difference in the abundance of immune cells between the high and low expression groups. Based on the median expression values of the hub gene, survival analysis was divided into high and low expression groups.

TIMER database

The TIMER database serves as an online visualization tool to assess the correlation between immune cell infiltration and target genes in different cancers (https://cistrome.shinyapps.io/timer/) [52]. In the current study, we analyzed the correlations between the infiltration abundance of B cells, CD8+ T cells, CD4+ T cells, macrophage, neutrophil and dendritic cells and expression levels of ANXA1 in sarcoma by TIMER database.

TISIDB database

The TISIDB database is an integrated repository portal for tumor-immune system interactions (http://cis.hku.hk/TISIDB/) [53]. Correlations between the infiltration abundance of tumor infiltrating lymphocytes (TILs) and the expression levels of ANXA1 in sarcoma evaluated by the TISIDB database.

Cell culture

ES is a malignant bone tumor that invades bone or muscle tissue, and mesenchymal stem cells are its main precursor cells [54]. Therefore, we selected RD-ES and A673 cells (Procell Bio, Wuhan, China) as the experimental group and MSCs as the control group for the experiment. Cells were maintained in Dulbecco’s modified Eagle’s medium (Gibco, NY, USA) supplemented with 10% fetal bovine serum (Gibco, NY, USA), 100 U/mL penicillin and 100 mg/mL streptomycin. Cultures were incubated at 37° C and 5% CO2.

RNA interference

Specific small interfering RNAs (siRNAs) were obtained from Guangzhou RiboBio Co., Ltd. ANXA1 siRNA or negative control (NC) siRNA was transfected into cells. The siRNA-ANXA1 sequences were as follows: #1, 5′-CAUAAGGCCAUAAUGGUUAAATT-3′; #2, 5′-UUUAACCAUUAUGGCCUUAUGTT-3′; #3, 5′-GCAUUCUAUCAGAAGAUGUAUTT-3′.

Cell survival test (CCK-8 assay)

Cells were inoculated into 96-well plates and cultured overnight, followed by transfection of cells with si-ANXA1 and si-NC. Ten microliters of CCK-8 were added to each well, and all experimental procedures were performed according to the instructions of the CCK-8 kit (Sigma-Aldrich, St. Louis, MO, USA).

Cell proliferation test (EdU assay)

Cells were inoculated into 24-well plates, and EdU reagent was added to each well at a ratio of 1:1000. The cells were fixed and stained according to the instructions of the EdU Cell Proliferation Kit (Beyotimebio, Shanghai, China), followed by observation of the cells using a fluorescence microscope.

RNA extraction and real-time quantitative PCR

The TRIzol (Ambion LLC, Austin, TX, USA) was used to extract total RNA from cells, which was reverse transcribed into cDNA using a reverse transcription kit (Service Bio, Guangzhou, China). Real-time quantitative PCR was performed using Universal Blue SYBR Green qPCR Master Mix (Service Bio, Guangzhou, China). After brief centrifugation, reverse transcription was performed on a general PCR instrument under the following conditions: denaturation at 95° C for 20 s, annealing at 55° C for 20 s, and extension at 72° C for 30 s. The expression levels of GAPDH were used as an internal control. The primer sequences are listed in Table 2.

Table 2. The primer used for hub gene and module genes.

Targeted geneForward (5′-3′)Reverse (3′-5′)
ANXA1CTTTCTCTTGCTAAGGGTGATGGTGGTAAGGATGGTATTG
COL1A2CAAAGGAGAGAGCGGTAACAGAAGACCACGAGAACCAGGA
S100A4CCACCTTCCACAAGTACTCGGCTTCATCTGTCCTTTTCCC
VIMTGACCGCTTCGCCAACTATTCGGCTTCCCCTCTCTG
S100A11CCCTGATTGCTGTCTTCCGGGTCCTTCTGGTTCTTT
MMP9ATGAGCCTCTGGCAGCCCCTGGTCCGGACCAGGGGCTGCCAGAGGCTCAT
GAPDHCCCATCACCATCTTCCAGGCATCACGCCACAGTTTCCC

Protein extraction and Western blotting analysis

RIPA lysis buffer (Servicebio, Wuhan, China) supplemented with phosphatase inhibitor (Servicebio, Wuhan, China) was used to extract total cellular proteins from cultured cells on ice for 10 min. The above lysates were centrifuged for 15 min at 16,000 rpm to remove cellular debris, and the supernatant containing proteins was collected and separated into 80 μl aliquots. A BCA Protein Assay Kit (Beyotimebio, Shanghai, China) was used to measure the protein concentration in the supernatant, and the separated proteins were transferred to PVDF membranes (Millipore, MA, USA). Subsequently, the membranes were blocked with 5% (w/v) skimmed milk (Servicebio, Wuhan, China) in Tris-buffered saline with Tween 20 (TBST; Solarbio, Beijing, China) for 60 min at 37° C. Primary antibodies were then probed overnight at 4° C in 1% (w/v) skimmed milk in TBST. After three washes with TBST, membranes were incubated with the appropriate horseradish peroxidase (HRP)-linked secondary antibodies for 60 min. Antibody reactivity was detected using the ECL Chemiluminescent substrate (Millipore, MA, USA). Gray values were analyzed using ImageJ software. Blots were representative of 3 independent experiments, with quantified results expressed as the means ± standard deviations (SD). The primary antibodies were as follows: ANXA1 (YT0234; Immunoway, Suzhou, China, 1:1500), MMP9 (YT1892; Immunoway, Suzhou, China, 1:1500), COL1A2 (YM4409; Immunoway, Suzhou, China, 1:1500), S100A4 (YM1458; Immunoway, Suzhou, China, 1:1000), VIM (YT4879; Immunoway, Suzhou, China, 1:1500), S100A11 (10237-1-AP; Proteintech, Wuhan, China, 1:1500), and β-Actain (66009-1-Ig; Proteintech, Wuhan, China, 1:25000).

Statistical analysis

The external validation of the genes was repeated three times. Statistical analysis was performed using GraphPad Prism (version 8.0), and Student’s t-test was used for comparisons between two groups. Bioinformatics analysis was performed using R 4.1.3.

Availability of data and materials

The microarray data used to support the findings of this study can be downloaded from the GSE17679, GSE63157, GSE45544 and GSE73166 datasets (https://www.ncbi.nlm.nih.gov/geo). The processed data are available from the corresponding author upon request.

Author Contributions

X S and Y X conceived the study and wrote the manuscript. Y W and L T carried out the data collection and data analysis. X S and H N contributed to the data curation, methodology, and validation. All authors reviewed the results and approved the final version of the manuscript.

Conflicts of Interest

All the authors declare that they have no conflicts of interest.

Funding

This study was funded by National Natural Science Foundation of China (Grant No.81772367, 82072392); the Yunnan Traumatology and Orthopedics Clinical Medical Center (Grant No. ZX20191001); the Grants from Yunnan Orthopedics and Sports Rehabilitation Clinical Medicine Research Center (Grant No. 202102AA310068).

Editorial Note

&

This corresponding author has a verified history of publications using a personal email address for correspondence.

References

  • 1. Balamuth NJ, Womer RB. Ewing’s sarcoma. Lancet Oncol. 2010; 11:184–92. https://doi.org/10.1016/S1470-2045(09)70286-4 [PubMed]
  • 2. Esiashvili N, Goodman M, Marcus RB Jr. Changes in incidence and survival of Ewing sarcoma patients over the past 3 decades: Surveillance Epidemiology and End Results data. J Pediatr Hematol Oncol. 2008; 30:425–30. https://doi.org/10.1097/MPH.0b013e31816e22f3 [PubMed]
  • 3. Jawad MU, Cheung MC, Min ES, Schneiderbauer MM, Koniaris LG, Scully SP. Ewing sarcoma demonstrates racial disparities in incidence-related and sex-related differences in outcome: an analysis of 1631 cases from the SEER database, 1973-2005. Cancer. 2009; 115:3526–36. https://doi.org/10.1002/cncr.24388 [PubMed]
  • 4. Ordóñez JL, Osuna D, Herrero D, de Alava E, Madoz-Gúrpide J. Advances in Ewing’s sarcoma research: where are we now and what lies ahead? Cancer Res. 2009; 69:7140–50. https://doi.org/10.1158/0008-5472.CAN-08-4041 [PubMed]
  • 5. Van der Graaf WTA, Orbach D, Judson IR, Ferrari A. Soft tissue sarcomas in adolescents and young adults: a comparison with their paediatric and adult counterparts. Lancet Oncol. 2017; 18:e166–75. https://doi.org/10.1016/S1470-2045(17)30099-2 [PubMed]
  • 6. Kauer M, Ban J, Kofler R, Walker B, Davis S, Meltzer P, Kovar H. A molecular function map of Ewing’s sarcoma. PLoS One. 2009; 4:e5415. https://doi.org/10.1371/journal.pone.0005415 [PubMed]
  • 7. Pappo AS, Dirksen U. Rhabdomyosarcoma, Ewing Sarcoma, and Other Round Cell Sarcomas. J Clin Oncol. 2018; 36:168–79. https://doi.org/10.1200/JCO.2017.74.7402 [PubMed]
  • 8. Alvarez-SanNicolas J, Gracia-Alegria I, Trullols-Tarrago L, Peiro-Ibañez A, Lamas-Gomez C. Prognostic factors and survival in Ewing’s sarcoma treated by limb salvage surgery. Clin Transl Oncol. 2019; 21:1374–82. https://doi.org/10.1007/s12094-019-02067-1 [PubMed]
  • 9. Rosen G, Wollner N, Tan C, Wu SJ, Hajdu SI, Cham W, D’Angio GJ, Murphy ML. Proceedings: Disease-free survival in children with Ewing’s sarcoma treated with radiation therapy and adjuvant four-drug sequential chemotherapy. Cancer. 1974; 33:384–93. https://doi.org/10.1002/1097-0142(197402)33:2<384::aid-cncr2820330213>3.0.co;2-t [PubMed]
  • 10. Van Mater D, Wagner L. Management of recurrent Ewing sarcoma: challenges and approaches. Onco Targets Ther. 2019; 12:2279–88. https://doi.org/10.2147/OTT.S170585 [PubMed]
  • 11. Gaspar N, Hawkins DS, Dirksen U, Lewis IJ, Ferrari S, Le Deley MC, Kovar H, Grimer R, Whelan J, Claude L, Delattre O, Paulussen M, Picci P, et al. Ewing Sarcoma: Current Management and Future Approaches Through Collaboration. J Clin Oncol. 2015; 33:3036–46. https://doi.org/10.1200/JCO.2014.59.5256 [PubMed]
  • 12. Stahl M, Ranft A, Paulussen M, Bölling T, Vieth V, Bielack S, Görtitz I, Braun-Munzinger G, Hardes J, Jürgens H, Dirksen U. Risk of recurrence and survival after relapse in patients with Ewing sarcoma. Pediatr Blood Cancer. 2011; 57:549–53. https://doi.org/10.1002/pbc.23040 [PubMed]
  • 13. Doyle LA. Sarcoma classification: an update based on the 2013 World Health Organization Classification of Tumors of Soft Tissue and Bone. Cancer. 2014; 120:1763–74. https://doi.org/10.1002/cncr.28657 [PubMed]
  • 14. Grünewald TGP, Cidre-Aranaz F, Surdez D, Tomazou EM, de Álava E, Kovar H, Sorensen PH, Delattre O, Dirksen U. Ewing sarcoma. Nat Rev Dis Primers. 2018; 4:5. https://doi.org/10.1038/s41572-018-0003-x [PubMed]
  • 15. Li W, Dong S, Lin Y, Wu H, Chen M, Qin C, Li K, Zhang J, Tang ZR, Wang H, Huo K, Xie X, Hu Z, et al. A tool for predicting overall survival in patients with Ewing sarcoma: a multicenter retrospective study. BMC Cancer. 2022; 22:914. https://doi.org/10.1186/s12885-022-09796-7 [PubMed]
  • 16. Li W, Jin G, Wu H, Wu R, Xu C, Wang B, Liu Q, Hu Z, Wang H, Dong S, Tang ZR, Peng H, Zhao W, Yin C. Interpretable clinical visualization model for prediction of prognosis in osteosarcoma: a large cohort data study. Front Oncol. 2022; 12:945362. https://doi.org/10.3389/fonc.2022.945362 [PubMed]
  • 17. Azimi F, Scolyer RA, Rumcheva P, Moncrieff M, Murali R, McCarthy SW, Saw RP, Thompson JF. Tumor-infiltrating lymphocyte grade is an independent predictor of sentinel lymph node status and survival in patients with cutaneous melanoma. J Clin Oncol. 2012; 30:2678–83. https://doi.org/10.1200/JCO.2011.37.8539 [PubMed]
  • 18. Sinnamon AJ, Sharon CE, Song Y, Neuwirth MG, Elder DE, Xu X, Chu EY, Ming ME, Fraker DL, Gimotty PA, Karakousis GC. The prognostic significance of tumor-infiltrating lymphocytes for primary melanoma varies by sex. J Am Acad Dermatol. 2018; 79:245–51. https://doi.org/10.1016/j.jaad.2018.02.066 [PubMed]
  • 19. Chin M, Yokoyama R, Sumi M, Okita H, Kawai A, Hosono A, Koga Y, Sano H, Watanabe H, Ozaki T, Mugishima H, and Japan Ewing Sarcoma Study Group (JESS). Multimodal treatment including standard chemotherapy with vincristine, doxorubicin, cyclophosphamide, ifosfamide, and etoposide for the Ewing sarcoma family of tumors in Japan: Results of the Japan Ewing Sarcoma Study 04. Pediatr Blood Cancer. 2020; 67:e28194. https://doi.org/10.1002/pbc.28194 [PubMed]
  • 20. Gargallo P, Juan A, Yáñez Y, Dolz S, Segura V, Castel V, Cañete A. Precision medicine in Ewing sarcoma: a translational point of view. Clin Transl Oncol. 2020; 22:1440–54. https://doi.org/10.1007/s12094-020-02298-7 [PubMed]
  • 21. Guo C, Liu S, Sun MZ. Potential role of Anxa1 in cancer. Future Oncol. 2013; 9:1773–93. https://doi.org/10.2217/fon.13.114 [PubMed]
  • 22. Araújo TG, Mota STS, Ferreira HSV, Ribeiro MA, Goulart LR, Vecchi L. Annexin A1 as a Regulator of Immune Response in Cancer. Cells. 2021; 10:2245. https://doi.org/10.3390/cells10092245 [PubMed]
  • 23. Bai F, Zhang P, Fu Y, Chen H, Zhang M, Huang Q, Li D, Li B, Wu K. Targeting ANXA1 abrogates Treg-mediated immune suppression in triple-negative breast cancer. J Immunother Cancer. 2020; 8:e000169. https://doi.org/10.1136/jitc-2019-000169 [PubMed]
  • 24. Rosengarth A, Luecke H. A calcium-driven conformational switch of the N-terminal and core domains of annexin A1. J Mol Biol. 2003; 326:1317–25. https://doi.org/10.1016/s0022-2836(03)00027-5 [PubMed]
  • 25. Boudhraa Z, Bouchon B, Viallard C, D’Incan M, Degoul F. Annexin A1 localization and its relevance to cancer. Clin Sci (Lond). 2016; 130:205–20. https://doi.org/10.1042/CS20150415 [PubMed]
  • 26. Silva-Oliveira R, Pereira FF, Petronilho S, Martins AT, Lameirinhas A, Constâncio V, Caldas-Ribeiro I, Salta S, Lopes P, Antunes L, Castro F, de Sousa SP, Henrique R, Jerónimo C. Clinical Significance of ARID1A and ANXA1 in HER-2 Positive Breast Cancer. J Clin Med. 2020; 9:3911. https://doi.org/10.3390/jcm9123911 [PubMed]
  • 27. Foo SL, Yap G, Cui J, Lim LHK. Annexin-A1 - A Blessing or a Curse in Cancer? Trends Mol Med. 2019; 25:315–27. https://doi.org/10.1016/j.molmed.2019.02.004 [PubMed]
  • 28. Sobral-Leite M, Wesseling J, Smit VT, Nevanlinna H, van Miltenburg MH, Sanders J, Hofland I, Blows FM, Coulson P, Patrycja G, Schellens JH, Fagerholm R, Heikkilä P, et al, and kConFab/AOCS Investigators. Annexin A1 expression in a pooled breast cancer series: association with tumor subtypes and prognosis. BMC Med. 2015; 13:156. https://doi.org/10.1186/s12916-015-0392-6 [PubMed]
  • 29. Su N, Xu XY, Chen H, Gao WC, Ruan CP, Wang Q, Sun YP. Increased expression of annexin A1 is correlated with K-ras mutation in colorectal cancer. Tohoku J Exp Med. 2010; 222:243–50. https://doi.org/10.1620/tjem.222.243 [PubMed]
  • 30. Mota STS, Vecchi L, Alves DA, Cordeiro AO, Guimarães GS, Campos-Fernández E, Maia YC, Dornelas BC, Bezerra SM, de Andrade VP, Goulart LR, Araújo TG. Annexin A1 promotes the nuclear localization of the epidermal growth factor receptor in castration-resistant prostate cancer. Int J Biochem Cell Biol. 2020; 127:105838. https://doi.org/10.1016/j.biocel.2020.105838 [PubMed]
  • 31. Petrella A, Festa M, Ercolino SF, Zerilli M, Stassi G, Solito E, Parente L. Annexin-1 downregulation in thyroid cancer correlates to the degree of tumor differentiation. Cancer Biol Ther. 2006; 5:643–7. https://doi.org/10.4161/cbt.5.6.2700 [PubMed]
  • 32. Rodrigo JP, Garcia-Pedrero JM, Fernandez MP, Morgan RO, Suárez C, Herrero A. Annexin A1 expression in nasopharyngeal carcinoma correlates with squamous differentiation. Am J Rhinol. 2005; 19:483–7. [PubMed]
  • 33. Chen P, Min J, Wu H, Zhang H, Wang C, Tan G, Zhang F. Annexin A1 is a potential biomarker of bone metastasis in small cell lung cancer. Oncol Lett. 2021; 21:141. https://doi.org/10.3892/ol.2020.12402 [PubMed]
  • 34. Ren B, Cui M, Yang G, Wang H, Feng M, You L, Zhao Y. Tumor microenvironment participates in metastasis of pancreatic cancer. Mol Cancer. 2018; 17:108. https://doi.org/10.1186/s12943-018-0858-1 [PubMed]
  • 35. Zhang Z, Shi Y, Zhu Z, Fu J, Liu D, Liu X, Dang J, Tao H, Fan H. Characterization of myeloid signature genes for predicting prognosis and immune landscape in Ewing sarcoma. Cancer Sci. 2022. [Epub ahead of print]. https://doi.org/10.1111/cas.15688 [PubMed]
  • 36. Zhang H, Maric I, DiPrima MJ, Khan J, Orentas RJ, Kaplan RN, Mackall CL. Fibrocytes represent a novel MDSC subset circulating in patients with metastatic cancer. Blood. 2013; 122:1105–13. https://doi.org/10.1182/blood-2012-08-449413 [PubMed]
  • 37. Gabrilovich DI, Ostrand-Rosenberg S, Bronte V. Coordinated regulation of myeloid cells by tumours. Nat Rev Immunol. 2012; 12:253–68. https://doi.org/10.1038/nri3175 [PubMed]
  • 38. Feng J, Lu SS, Xiao T, Huang W, Yi H, Zhu W, Fan S, Feng XP, Li JY, Yu ZZ, Gao S, Nie GH, Tang YY, Xiao ZQ. ANXA1 Binds and Stabilizes EphA2 to Promote Nasopharyngeal Carcinoma Growth and Metastasis. Cancer Res. 2020; 80:4386–98. https://doi.org/10.1158/0008-5472.CAN-20-0560 [PubMed]
  • 39. Zhu JF, Huang W, Yi HM, Xiao T, Li JY, Feng J, Yi H, Lu SS, Li XH, Lu RH, He QY, Xiao ZQ. Annexin A1-suppressed autophagy promotes nasopharyngeal carcinoma cell invasion and metastasis by PI3K/AKT signaling activation. Cell Death Dis. 2018; 9:1154. https://doi.org/10.1038/s41419-018-1204-7 [PubMed]
  • 40. Zhang J, Huang J, Liu W, Ding L, Cheng D, Xiao H. Identification of Common Oncogenic Genes and Pathways Both in Osteosarcoma and Ewing’s Sarcoma Using Bioinformatics Analysis. J Immunol Res. 2022; 2022:3655908. https://doi.org/10.1155/2022/3655908 [PubMed]
  • 41. Tang M, Liu P, Wu X, Gong J, Weng J, Gao G, Liu Y, Gan L. COL3A1 and Its Related Molecules as Potential Biomarkers in the Development of Human Ewing’s Sarcoma. Biomed Res Int. 2021; 2021:7453500. https://doi.org/10.1155/2021/7453500 [PubMed]
  • 42. Schaefer KL, Brachwitz K, Wai DH, Braun Y, Diallo R, Korsching E, Eisenacher M, Voss R, Van Valen F, Baer C, Selle B, Spahn L, Liao SK, et al. Expression profiling of t(12;22) positive clear cell sarcoma of soft tissue cell lines reveals characteristic up-regulation of potential new marker genes including ERBB3. Cancer Res. 2004; 64:3395–405. https://doi.org/10.1158/0008-5472.CAN-03-0809 [PubMed]
  • 43. Sancéau J, Truchet S, Bauvois B. Matrix metalloproteinase-9 silencing by RNA interference triggers the migratory-adhesive switch in Ewing’s sarcoma cells. J Biol Chem. 2003; 278:36537–46. https://doi.org/10.1074/jbc.M304300200 [PubMed]
  • 44. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, Smyth GK. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015; 43:e47. https://doi.org/10.1093/nar/gkv007 [PubMed]
  • 45. Morris JA, Gayther SA, Jacobs IJ, Jones C. A suite of Perl modules for handling microarray data. Bioinformatics. 2008; 24:1102–3. https://doi.org/10.1093/bioinformatics/btn085 [PubMed]
  • 46. Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012; 16:284–7. https://doi.org/10.1089/omi.2011.0118 [PubMed]
  • 47. Gene Ontology Consortium. The Gene Ontology (GO) project in 2006. Nucleic Acids Res. 2006; 34:D322–6. https://doi.org/10.1093/nar/gkj021 [PubMed]
  • 48. Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics. 2013; 14:7. https://doi.org/10.1186/1471-2105-14-7 [PubMed]
  • 49. Szklarczyk D, Franceschini A, Kuhn M, Simonovic M, Roth A, Minguez P, Doerks T, Stark M, Muller J, Bork P, Jensen LJ, von Mering C. The STRING database in 2011: functional interaction networks of proteins, globally integrated and scored. Nucleic Acids Res. 2011; 39:D561–8. https://doi.org/10.1093/nar/gkq973 [PubMed]
  • 50. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, Amin N, Schwikowski B, Ideker T. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003; 13:2498–504. https://doi.org/10.1101/gr.1239303 [PubMed]
  • 51. Jin Y, Wang Z, He D, Zhu Y, Chen X, Cao K. Identification of novel subtypes based on ssGSEA in immune-related prognostic signature for tongue squamous cell carcinoma. Cancer Med. 2021; 10:8693–707. https://doi.org/10.1002/cam4.4341 [PubMed]
  • 52. Li T, Fan J, Wang B, Traugh N, Chen Q, Liu JS, Li B, Liu XS. TIMER: A Web Server for Comprehensive Analysis of Tumor-Infiltrating Immune Cells. Cancer Res. 2017; 77:e108–10. https://doi.org/10.1158/0008-5472.CAN-17-0307 [PubMed]
  • 53. Ru B, Wong CN, Tong Y, Zhong JY, Zhong SSW, Wu WC, Chu KC, Wong CY, Lau CY, Chen I, Chan NW, Zhang J. TISIDB: an integrated repository portal for tumor-immune system interactions. Bioinformatics. 2019; 35:4200–2. https://doi.org/10.1093/bioinformatics/btz210 [PubMed]
  • 54. Tirode F, Laud-Duval K, Prieur A, Delorme B, Charbord P, Delattre O. Mesenchymal stem cell features of Ewing tumors. Cancer Cell. 2007; 11:421–9. https://doi.org/10.1016/j.ccr.2007.02.027 [PubMed]
Download PDF
Corresponding Author

Yongqing Xu
docxuyongqing@sina.com
https://orcid.org/0000-0001-8724-1111

Keywords
bone tumor inflammation immunotherapy
Table of Contents