Exploration and Validation of a Novel Inflammatory Response-Associated Gene Signature to Predict Osteosarcoma Prognosis and Immune Infiltration

Introduction

Osteosarcoma was the predominant malignant primary bone tumor that mainly occurred in adolescents. It was the second leading cause of tumor-associated death in teenagers.^1,2 In order to overcome the lethality of osteosarcoma, the radical treatment which included multi-drug neoadjuvant chemotherapy and aggressive surgical resection was applied in clinical. Since then, the 5-year survival rates of osteosarcoma patients were significantly improved.³ However, numerous patients still developed metastasis before or after treatment and their survival rates dramatically decreased to less than 20%.⁴ To improve the prognosis of these metastatic patients, novel therapeutic strategies such as immune checkpoint and molecular targeted inhibitors had been experimentally applied in clinical. Unfortunately, little progress had been achieved since the mechanism of metastasis was still elusive.

Inflammatory response, always accompanied with immune alterations in the tumor microenvironment, was involved in the development of tumor. In 1863, Virchow first found that “lymphoreticular infiltrate” was existed in the origin of tumor and this chronic inflammatory response might promote the progression of tumor.⁵ After that, the relationships between inflammation and cancer were widely investigated. To some extent, cancer was considered as the wound that did not heal and chronic inflammatory response would cause gene mutation and epigenetic alterations, which promoted the malignancy of cancer cells.^6,7 Numerous studies also revealed that clinical blood inflammatory parameters were associated with the prognosis of cancers. In osteosarcoma, several elevated blood inflammatory-associated markers such as serum C-reactive protein (CRP) and erythrocyte sedimentation rate (ESR) were negatively related to patients’ prognosis.^8,9

In addition, plenty of inflammatory response-associated genes were identified by microarray and next-generation sequencing technology in recent years.¹⁰ Several studies revealed that these genes could be combined as the signature to predict the prognosis of patients with malignant tumors and might be applied in clinical treatment.^11,12 However, few studies focused on the relationship between osteosarcoma and inflammatory response.

In the present study, we identified inflammatory response-associated differentially expressed messenger RNA (DEmRNA) related to osteosarcoma metastasis in Gene Expression Omnibus (GEO) database. Univariate and multivariate cox regression analyses were conducted to screen out the prognosis-related genes and establish the prognostic signature. Next, patients were divided into low- and high-risk groups according to the median value of risk scores obtained previously. The Kaplan–Meier (KM) survival analysis was performed and the relationships between individual genes and clinical features were explored. In addition, the nomogram was established and the immune infiltrative features of signature were investigated. All the above results were validated in Therapeutically Applicable Research to Generate Effective Treatments (Target) database. The DEmRNAs in signature were also validated in different cell lines by quantitative real-time PCR (qRT-PCR) and in osteosarcoma patients’ samples by immunohistochemistry (IHC). The flow chart of this study is shown in Figure 1. We believed that our research constructed a reliable gene signature to reveal the relationship between inflammatory response and osteosarcoma patients’ prognosis. We hoped that this signature could act as the novel prognostic biomarker and would be applied in clinical practice in the future.

Figure 1 The flow chart of the present study. (A) The osteosarcoma inflammatory response-associated DEmRNAs were identified in GSE21257 and MSigDB databases. (B) The inflammatory response-associated gene signature was established through univariate and multivariate cox regression analyses. (C) The application of inflammatory response-associated gene signature in clinic and validation in Target database. (D) Identification and functional analysis of DEmRNAs related to inflammatory response-associated risk. (E) The association between immune infiltration and inflammatory response. Abbreviations: DEmRNA, differentially expressed messenger RNA; MSigDB, Molecular Signature Database.

Materials and Methods

Screening of Datasets

The datasets of osteosarcoma patients were searched from Target (https://ocg.cancer.gov/programs/target) and GEO (http://www.ncbi.nlm.nih.gov/geo) databases. The inclusion criteria for dataset were shown as follows: (1) the patients were diagnosed as osteosarcoma by pathology; (2) the datasets should contain complete prognostic information; (3) the number of patients in datasets should be more than 50. Finally, one GEO dataset (GSE21257, platform GPL10295, Illumina human-6 v2.0 expression beadchip Illumina, Inc., San Diego, CA, United States) was screened out as the training group and another dataset in Target database was selected as the validation group. The clinical characteristics of patients in two datasets are shown in Supplementary Table 1. The inflammatory response-associated genes were obtained from the hallmark gene sets in the Molecular Signature Database (MSigDB, https://www.gsea-msigdb.org/gsea/msigdb/).

Identification of Inflammatory Response-Associated DEmRNAs

The mRNA expression file and patients’ clinical information were downloaded from each dataset. Patients in GSE21257 were separated into metastatic and non-metastatic groups and the DEmRNAs among these groups were screened out through Linear Models for Microarray Analysis (limma) package in R software. The cutoff values of screening were |fold change (FC)| >1.5 and P < 0.05. The characteristic distribution of DEmRNAs was shown by “pheatmap” package in R software. Inflammatory response-associated DEmRNAs were identified through overlapping DEmRNAs and inflammatory response-associated genes.

Construction and Validation of the Inflammatory Response-Associated Prognostic Signature

Inflammatory response-associated DEmRNAs related to osteosarcoma prognosis were identified by univariate cox regression analysis. The multivariate cox regression analysis was utilized to construct the inflammatory response-associated gene signature. The least numbers of DEmRNAs that represented the signature were screened through akaike information criterion (AIC).¹³ The risk score of signature was calculated by the following formulation:

Risk score = ∑i Coefficient mRNAi * Expression mRNAi

Then the patients in GSE21257 were divided into high- and low-risk groups according to the median value of risk score. The relationships between risk score and patients’ survival were conducted by KM survival curve in R software. The reliability of the signature was evaluated by time-related receiver operating characteristic (ROC) curves.

After that, the univariate and multivariate cox regression analyses were applied to evaluate the prognostic value of this gene signature. The correlations between individual mRNAs in signature and clinical features were explored through the follow-up information. The KM survival analysis of mRNAs in signature was also conducted. The above features of signature were validated in Target database and the KM survival analysis of individual mRNAs was examined in R2 database (https://hgserver1.amc.nl/cgi-bin/r2/main.cgi) once more. P < 0.05 was considered as statistically significant.

Establishment of the Nomogram

Recently, the nomogram that contained patients’ clinical information was constructed to stratify the different risk patients in various studies.¹⁴ The nomogram of the osteosarcoma patients in the present study was established through “rms” package in R software. The predictive value of clinical factors such as gender, age and metastasis were calculated by multivariate cox regression analysis. The reliability of the nomogram was evaluated by Harrell’s concordance index (C-index). The nomogram was also established and validated in Target database.

Functional Analyses of DEmRNAs Related to Inflammatory Response-Associated Prognostic Signature

The DEmRNAs related to this gene signature were screened out through high- and low-risk groups in Target and GSE21257 datasets. The overlapped DEmRNAs among these two datasets were considered as the risk score-related DEmRNAs. The heatmaps of DEmRNAs in two datasets were drawn by the “pheatmap” package in R software. The functional analyses of these DEmRNAs were performed by “clusterProfiler” and “enrichplot” packages in R software. P- and q-values less than 0.05 were defined as statistically significant.

Assessment of Tumor Immune Infiltration Status Related to the Inflammatory Response

The infiltrative levels of immune and stromal cells in osteosarcoma samples were evaluated through immune and stromal scores by ESTIMATE software. The relationships between these scores and inflammatory response-associated risks were conducted by Spearman correlate analysis.¹⁵ Next, the compositions of 22 immune cells in osteosarcoma were explored through the CIBERSORT deconvolution algorithm.^16,17 In brief, the mRNA signature matrix of 22 immune cells was downloaded from the CIBERSORT platform (https://cibersortx.stanford.edu). Then, it was matched with the mRNAs from GSE21257 and Target datasets to generate the infiltrative proportions of immune cells in osteosarcoma.

After that, the immune infiltrative differences between low- and high-risk patients were analyzed by single-sample gene set enrichment analysis (ssGSEA). The intersected differences with the same alternative trends in GSE21257 and Target datasets were regarded as the inflammatory response-associated immune changes in osteosarcoma. The immune infiltrative features between lowly and highly expressed genes in signature were also investigated.

Cell Culture and Reagents

Human osteosarcoma cell lines 143B, MNNG/HOS, SAOS2, U2OS, MG63 and human osteoblast cell line hFOB1.19 were obtained from ATCC (Manassas, VA, USA). Osteosarcoma cell line WELL5 was a kind gift from Dr. Wu Zhang of Shanghai Institute of Hematology.¹⁸ All osteosarcoma cell lines were cultured in high glucose DMEM (Gibco, 11965092) supplemented with 10% fetal bovine serum (FBS, Gibco, 10099–141) and 1% penicillin-streptomycin (Gibco, 15070063) at 37°C in 5% CO₂. The hFOB1.19 cell was cultured in DMEM/F12 medium (Gibco, 11320033) containing 10% FBS and 0.3mg/mL G418 (Gibco, 10131035) at 34°C in 5% CO₂.

Quantitative Real-Time PCR Analysis

The mRNA expressions of MYC and CLEC5A in osteosarcoma and osteoblast cell lines were quantified by quantitative real-time PCR (qRT-PCR) analyses. In brief, whole cellular RNA was extracted by Trizol reagent (Invitrogen, 15596018) and the cDNA templates were obtained by reverse transcription. The amplification was conducted through Roche LightCycler System. The whole procedure was performed according to the protocol of TB Green Premix Ex Taq II kit (Takara RR820A, Japan). The primer sequences of genes were: MYC: 5ʹ- CGACGAGACCTTCATCAAAAAC-3ʹ (forward) and 5ʹ-CTTCTCTGAGACGAGCTTGG-3ʹ (reverse); CLEC5A: 5ʹ-GAAAAGGATCCACATTGGCAAT-3ʹ (forward) and 5ʹ-GTCGCACAGTTGAAATTCTGA-3ʹ (reverse); GAPDH: 5ʹ-GCACCGTCAAGGCTGAGAAC-3ʹ (forward) and 5ʹ-TGGTGAAGACGCCAGTGGA-3ʹ (reverse). The expression of genes was calculated by formula 2^−ΔΔCT. The expressive difference between osteosarcoma and osteoblast cell lines or high- and low-metastatic osteosarcoma cell lines were calculated by Student’s t test. P value less than 0.05 was considered as statistically significant.

Immunohistochemistry Analysis

Ten patients who were diagnosed as osteosarcoma by pathology in Ruijin Hospital were selected for IHC assay. Among these patients, five developed metastasis in 3 years and the others lived without metastasis. Each osteosarcoma sample was fixed with 10% formalin, embedded in paraffin and sliced as 4 μm thickness. The sections were deparaffinized and rehydrated by xylene and ethanol. After heat-induced antigen retrieval, endogenous peroxidase activity blockage and serum sealing, the samples were incubated with rabbit anti-CLEC5A (Biosis, bs-2663R, 1:200) and anti-MYC (Affinity, AF0358, 1:100) antibodies at 4°C overnight. Then, slices were covered with horseradish peroxidase-coupled goat anti-rabbit secondary antibody (CST, 7074, 1:200) at room temperature for 50 minutes and stained with diaminobenzidine. The nuclei were counterstained by hematoxylin. At last, the samples were dehydrated and mounted. All slices were visualized by the microscope (Nikon DS-U3).

Results

Identification of Inflammatory Response-Associated DEmRNAs Related to Osteosarcoma Metastasis

The expressions of mRNAs in GSE21257 that included 34 metastatic and 19 non-metastatic osteosarcoma patients were downloaded from GEO database. Two hundred and forty-seven up-regulated and 217 down-regulated DEmRNAs related to osteosarcoma metastasis were identified (Figure 2A). Two hundred inflammatory response-associated genes were obtained from hallmark gene sets in MSigDB and then overlapped with previously obtained DEmRNAs (Figure 2B). In the end, 24 inflammatory response-associated DEmRNAs were screened out (Table 1).

Table 1 The Inflammatory Response-Associated DEmRNAs in GSE21257

Figure 2 Identification of inflammatory response-associated DEmRNAs related to osteosarcoma metastasis. (A) DEmRNAs related to osteosarcoma metastasis in GSE21257 database. Red dots indicated up-regulated DEmRNAs (P < 0.05, FC > 1.5) and green dots indicated down-regulated DEmRNAs (P > 0.05, FC < −1.5). (B) Twenty-four inflammatory response-associated DEmRNAs were screened out in GSE21257 and MSigDB databases. Abbreviations: DEmRNA, differentially expressed messenger RNA; FC, fold change; MSigDB, Molecular Signature Database.

Construction and Validation of the Inflammatory Response-Associated Prognostic Signature in GEO and Target Databases

Based on the 24 DEmRNAs acquired previously, univariate and multivariate cox regression analyses were utilized to establish the inflammatory response-associated prognostic signature. In univariate cox regression analysis, three genes were associated with osteosarcoma prognosis (Figure 3A). Then, two of these three DEmRNAs were screened out to construct the inflammatory response-associated prognostic signature according to the multivariate cox regression result (Table 2). The risk score formula was established through mRNAs’ regression coefficients multiplied by expression levels:

Table 2 The Multivariate Cox Regression Analysis of GSE21257

Figure 3 Construction and validation of the inflammatory response-associated prognostic signature in GSE21257 and Target datasets. (A) Three inflammatory response-associated prognostic genes were identified in univariate cox regression analysis. (B and E) Upper: The distribution and median value of risk scores in GSE21257 (B) and Target (E) datasets. Bottom: The scatter plot of patients’ risk scores and survival status in GSE21257 (B) and Target (E) datasets. (C and F) Kaplan-Meier survival plot of patients in high- and low-risk groups in GSE21257 (C) and Target (F) datasets. (D and G) Time-dependent ROC curve at 1, 3, and 5 years of inflammatory-associated prognostic signature in GSE21257 (D) and Target (G) datasets. Abbreviations: Target, Therapeutically Applicable Research to Generate Effective Treatments; ROC, receiver operating characteristic; AUC, area under curve.

Risk score = (0.41*MYC) + (−0.96*CLEC5A)

Based on the median value of risk score, 26 patients were assigned to high-risk group and the rest 27 were regarded as low-risk. At the same time, patients in high-risk group showed worse mortality rates and shorter survival times than low-risk patients (Figure 3B). The KM survival curve also indicated poorer survival rates in high-risk patients (Figure 3C). Finally, the area under the time-dependent ROC curves (AUC) were 0.65 at 1-year, 0.80 at 3- and 5-year, which demonstrated the robustness of this inflammatory response-associated prognostic signature (Figure 3D).

The reliability of the inflammatory response-associated gene signature was validated in Target database. Ninety-five patients with integrated follow-up information were included in the test. Based on the median value of risk scores, 48 patients were separated into low-risk group and the rest were considered as high risk (Figure 3E). The risk plot (Figure 3E) and KM survival analysis (Figure 3F) showed that low-risk patients had better prognoses. In addition, the time-dependent ROC curve revealed the signature was robust (Figure 3G).

Exploration of the Relationship Between Clinical Features and Inflammatory Response-Associated Prognostic Signature

The prognostic value of gene signature was investigated through univariate and multivariate cox regression analyses. Figure 4A and B showed that the risk score was the independent prognostic factor in osteosarcoma. The same result was validated in Target database (Supplementary Figure 1A and 1B).

Figure 4 Relationships between inflammatory response-associated gene signature and clinical parameters in GSE21257. (A and B). Univariate and multivariate cox regression analyses revealed the risk score was the independent prognostic factor in GSE21257. (C–F). Boxplots of the relationship between CLEC5A expression and clinical characteristics. (G–J). The correlation between MYC expression and clinical characteristics. Abbreviation: CLEC5A, C-type lectin domain family 5 member A.

After that, the correlations between clinical features and genes in signature were investigated. In GSE21257, CLEC5A was associated with patients’ metastasis status as well as inflammatory response-associated risk and MYC was related to gender, metastasis status and inflammatory response-associated risk (Figure 4C–J). In Target database, CLEC5A and MYC were only correlated to risk scores (Supplementary Figure 1C–J). The KM survival analyses indicated the expression of CLEC5A positively related to the survival rates of osteosarcoma patients in GSE21257, Target and online R2 databases (Supplementary Figure 2A, C, E). Instead, the higher expression of MYC related to the worse prognosis of osteosarcoma patients (Supplementary Figure 2B, D, F).

Establishment of the Nomogram

The nomogram was a powerful tool that integrated clinical variables to evaluate the prognosis of patients. In this research, the nomogram was constructed in GSE21257 and Target databases to predict the 3- and 5-year survival rates of osteosarcoma patients (Figure 5A and Supplementary Figure 3A). The 3-year and 5-year calibration curves were close to the standard curve in two datasets (Figure 5B, C and Supplementary Figure 3B, C). The C-index was 0.99 in GSE21257 and 0.77 in Target dataset, which indicated the high accuracy of the predictive model.

Figure 5 Establishment of the nomogram in GSE21257. (A) The nomogram predicted the 3- and 5-year survival risk in osteosarcoma patients. (B) The calibration curve of the 3-year survival. (C) The calibration curve of the 5-year survival.

Functional Analysis of DEmRNAs Related to Inflammatory Response-Associated Prognostic Signature

The DEmRNAs among high- and low-risk patients were explored in GSE21257 and Target databases. A total of 736 DEmRNAs were screened out in GSE21257 and 1619 DEmRNAs were identified in Target database (Figure 6A and B). Then, 216 DEmRNAs overlapped in two databases were selected for further analysis (Figure 6C).

Figure 6 Identification and functional analyses of DEmRNAs related to inflammatory response-associated risks. (A) The heatmap of DEmRNAs associated with risk score in GSE21257. (B) The heatmap of risk score-associated DEmRNAs in Target dataset. (C) The overlapped DEmRNAs in GSE21257 and Target datasets. (D and E) The top 10 GO functional and top 30 KEGG pathway analyses of DEmRNAs. Abbreviations: Target, Therapeutically Applicable Research to Generate Effective Treatments; DEmRNA, differentially expressed messenger RNA; GO, gene ontology; KEGG, Kyoto encyclopedia of genes and genomes.

The biological functions of these DEmRNAs were investigated through GO functional and KEGG pathway analyses (Figure 6D and E). A number of inflammatory response-associated biological processes such as neutrophil activation, regulation of leukocyte activation and antigen processing and presentation were significantly enriched. Furthermore, some immune-related functions (ie, regulation of immune effector process, Th17 cell differentiation) were also enriched, which implied the underlying relationships between inflammatory response and immune alterations in osteosarcoma.

Evaluation of the Correlation Between Inflammatory Response and Immune Infiltration

The relationship between tumor microenvironment and inflammatory response-associated risks was evaluated through ESTIMATE software. The results showed that the risk score was negatively related to the immune and stromal scores in GSE21257 and Target datasets (Figure 7A, B and Supplementary Figure 4A, B). Then the immune infiltration of 22 immune cells was investigated by CIBERSORT algorithm. The statistically significant abundance ratios of immune cells in osteosarcoma were shown in Figure 7C and Supplementary Figure 4C. Then the specific relationships between immune functions and inflammatory response-associated risks were evaluated through ssGSEA analysis. The outcome revealed that the vast majority of immune cells (14) and functions (11) were down-regulated in high-risk group in GSE21257 (Figure 7D and E). At the same time, eight different immune cells and twelve immune functions were dysregulated in Target dataset (Supplementary Figure 4D, E). The overlapped items in both datasets that included eight immune cells and eleven immune functions were considered as the alterations of inflammatory response-associated immune infiltration (Figure 7F). In addition, the expression of MYC was negatively related to 10 immune cells and 10 immune functions in GSE21257, while CLEC5A showed the positive correlations (Supplementary Figure 5A–D). The same dysregulated trends were observed in Target dataset (Supplementary Figure 5E–H).

Figure 7 The relationships between inflammatory response and immune infiltration in GSE21257. (A) The correlation between risk score and immune score. (B). The relationship between risk score and stromal score. (C) The infiltrative proportion of 22 immune cells in osteosarcoma with statistically significant. (D) The correlations between risk scores and different immune cells. (E) The association between risk scores and different immune functions. (F) The overlapped immune cells and functions in GSE21257 and Target datasets. *P < 0.05; **P < 0.01; ***P < 0.001. Abbreviations: Target, Therapeutically Applicable Research to Generate Effective Treatments; NS, not significant; aDCs, activated dendritic cells; DCs, dendritic cells; iDCs, immature dendritic cells; NK, natural killer; aDCs, activated dendritic cells; pDCs, plasmacytoid dendritic cells; Tfh, T-follicular helper cells; Th, T helper; TIL, tumor infiltrating lymphocytes; Treg, T regulatory cells; APC, antigen-presenting cells; CCR, cytokine-cytokine receptor interaction; HLA, human leukocyte antigen; MHC, major histocompatibility complex; IFN, interferon.

Validation of the Inflammatory Response-Associated Genes in Cell Lines and Osteosarcoma Samples

The expressive levels of inflammatory response-associated genes in human osteosarcoma and osteoclast cell lines were detected through qRT-PCR analyses. The results showed that MYC was up-regulated and CLEC5A was down-regulated in all osteosarcoma cell lines (Figure 8A and B). Then, osteosarcoma cell lines were divided into high- (143B, WELL5, MNNG/HOS) and low-metastatic (U20S, SAOS2, MG63) groups according to their metastasis ability. The MYC was up-regulated in high-metastatic osteosarcoma cell lines and CLEC5A showed the opposite trend but with no statistically significant (Figure 8C and D). At the same time, the expressions of MYC and CLEC5A in metastasis and non-metastasis patients were validated by IHC. The result revealed that MYC was down-regulated in non-metastasis patients (Figure 8E) while CLEC5A showed the opposite trends (Figure 8F).

Figure 8 The expression of MYC and CLEC5A in cell lines and osteosarcoma patients’ samples. (A and B) The expression of MYC (A) and CLEC5A (B) in osteosarcoma and osteoblast cell lines. (C and D) The expression of MYC (C) and CLEC5A (D) in high- and low-metastatic cell lines. (E and F) The expression of MYC (E) and CLEC5A (F) in metastasis and non-metastasis patients, scale bar: 20μm. CLEC5A: C-type lectin domain family 5 member A; *P < 0.05; **P < 0.01; ****P < 0.0001.

Discussion

Osteosarcoma, one of the most common bone malignant tumors, often developed lung metastasis and the 5-year survival rates of these patients were extremely low.² In clinical studies, the diagnosis of osteosarcoma metastasis mainly relied on the radiology tests such as X-ray, computerized tomography (CT) and positron emission tomography-computed tomography (PET-CT). Unfortunately, the sensitivity and specificity of these methods were not high enough and the patients with metastatic possibilities always lost the best chance to receive early treatment. Thus, the novel reliable evaluation of early-stage metastasis was urgently needed.

In 1889, Paget first described the mechanism of tumor metastasis which was widely known as the “seed and soil” theory.¹⁹ After that, various theories related to metastasis were discovered, such as epithelial-mesenchymal transition (EMT), angiogenesis, gene mutation and so on.²⁰ Among these hallmarks, the inflammatory response was extensively investigated. Numerous inflammatory biomarkers, such as neutrophils, lymphocytes, platelets and so on, had been implied as the prognostic predictor for different cancers.^21–23 In osteosarcoma, some researchers revealed that lymphocyte-to-monocyte ratio (LMR) was positively related to overall and event-free survival rates of patients, while the neutrophil-lymphocyte ratio (NLR) and platelet-to-lymphocyte ratio (PLR) showed the opposite trends.^24–26 Several studies also implied that nonsteroidal anti-inflammatory drugs (NSAIDs) could inhibit the progression of osteosarcoma.^27,28 In addition, the relationship between inflammatory response-associated genes and osteosarcoma had also been widely studied. For instance, Jiang et al found that the inflammatory gene polymorphisms were closely related to the occurrence and progression of osteosarcoma.²⁹ Recent researches also indicated that several inflammatory response-associated genes could combine as the prognostic signature in different cancers.^11,12 Unfortunately, none of this signature had been established in osteosarcoma.

In the present study, we first analyzed the relationship between osteosarcoma and inflammatory response-associated genes. We screened out 464 DEmRNAs related to osteosarcoma metastasis in GSE21257 dataset. By overlapping with MSigDB database, 24 of these DEmRNAs were identified as the inflammatory response-associated genes. Then, two of these genes were selected to construct the inflammatory response-associated gene signature through univariate and multivariate cox regression analyses. Next, patients in GSE21257 were separated into high- and low-risk groups according to the median value of risk scores. The KM survival curve indicated that the risk score was negatively correlated with patients’ survival rates and the AUC values of ROC curves revealed that this signature was reliable. Most importantly, the same results were obtained in Target database. Besides that, we verified that the risk score was the independent prognostic factor for osteosarcoma patients.

The nomogram was a convincible tool that incorporated clinical features to estimate the patients’ prognosis. It could stratify the patients with unfavorable prognoses and guide clinical treatments. In this research, the nomogram was established in both datasets and the calibration curves, as well as C indexes, indicating the robustness of this nomogram. To sum up, we believed that this inflammatory response-associated gene signature was a reliable prognostic predictor for osteosarcoma and it could be easier to apply in clinics because only two genes were involved in it.

In the last decades, numerous studies demonstrated that the tumor infiltrative immune cells and immune-related genes played vital roles in the development and progression of tumors.³⁰ In osteosarcoma, tumor cells could also recruit different immune infiltrating cells and establish the specific local immune microenvironment, which finally promoted tumor growth, drug resistance and metastasis.³¹ For instance, Cao et al found that the expression of BMPR2 was negatively related to the infiltration of CD8+ T cells, monocytes, and M2 macrophages, which resulted in the poor prognosis of patients.³² In addition, Deng et al revealed that after neoadjuvant chemotherapy, more immune cells such as CD3+ T cells, CD8+ T cells and Ki67+/CD8+ T cells were infiltrated in tumor microenvironment and this might benefit for immunotherapy.³³ However, no studies focused the relationship between inflammatory response and immune infiltration in osteosarcoma.

In this study, we screened out the overlapped DEmRNAs related to inflammatory response-associated risks in GSE21257 and Target datasets. The functional analyses of these genes indicated that they were involved in the dysregulated immune functions such as regulation of immune effector process, Th17 cell differentiation and so on. According to the above results, we speculated that the inflammatory response might transform the tumor microenvironment, especially for the immune system. We then used Estimate algorithm to evaluate the correlation between inflammatory response-associated risks and stromal or immune infiltration. Both stromal and immune scores were negatively related to risk scores, which indicated higher purity of osteosarcoma cells existed in high-risk groups. These results implied that high-risk osteosarcoma was poorly immunogenic infiltrated, which was widely known as immunological cold tumors and thus was not sensitive to immunogenic therapy.³⁴ To further investigate the specific dysregulation of immune functions, ssGSEA analysis was conducted. The results showed that high-risk patients related to decreased infiltrative immune cells and functions. These conclusions were validated in various studies. For example, Muraro et al found that osteosarcoma cells could inhibit the maturation and function of dendritic cells (DCs) and these effects would be revised by adding inflammatory-associated cytokines (rhIL-12) or compounds (Indometacin).³⁵ The follicular helper T (Tfh) cells were related to the prognosis of patients in several tumors.³⁶ In osteosarcoma, Gao et al revealed that impaired functions of Tfhs lead to lower IL-21 secretion in CD4+ T cells, which might result in the immunosuppressive microenvironment.³⁷ In a word, the inflammatory response was closely related to the osteosarcoma immune infiltrative alteration and the down-regulated immune functions might contribute the tumor development.

What is more, the expression of individual mRNAs in the signature was investigated. The results showed that both of them related to the survival rates of patients in GSE21257, Target and R2 online datasets. For in vitro validation, the qRT-PCR analyses revealed that the expression of MYC was up-regulated in all osteosarcoma cell lines while CLEC5A showed the opposite trend. Compared to low-metastatic cell lines, MYC was over-expressed in high-metastatic cell lines but CLEC5A exhibited the contrary tendency. For the IHC study, the same expressive trends of two genes were observed between metastasis and non-metastasis osteosarcoma patients’ samples. Besides that, the ssGSEA analysis revealed that in both datasets, CLEC5A was positively related to the immune cells’ infiltration and immune functions while MYC showed the opposite trends. These results indicated that the mRNAs in our inflammatory response-associated signature played vital roles in the dysregulation of immune system, which finally promoted tumor progression.

MYC, one of the most common oncogene, was mutated in more than 10% of osteosarcoma patients.³⁸ High expression of MYC promoted proliferation, migration and spheroid formation of osteosarcoma cells and it was also upregulated in metastatic samples.^39,40 In addition, Kortlever et al revealed that MYC promoted the inflammatory, angiogenic and immunosuppressive tumor microenvironment by elevating IL-23 and CCL9, which accelerated pulmonary tumor progression. C-type lectin domain family 5 member A (CLEC5A) was one of the C-type lectin receptors that involved in inflammatory diseases.⁴¹ It regulated the release of pro-inflammatory cytokines from macrophages in infection and was involved in the process of collagen-induced arthritis.⁴² It played vital roles in the innate immune system by regulating the differentiation and activation of macrophages as well as neutrophils.⁴³ Lu et al found that the expression of CLEC5A in osteosarcoma was higher than adjacent normal tissues and it could prevent the metastasis of osteosarcoma.⁴⁴ However, the specific mechanisms of CLEC5A in osteosarcoma progression were still unknown and further research should be carried out.

For all we know, this research first revealed the relationship between inflammatory-associated gene signature and osteosarcoma prognosis. However, several defects still existed. First, the sample size of datasets was relatively small. This might owe to the extremely low incidence of osteosarcoma and the large cohort researches in osteosarcoma were limited. Second, the correlations between immune infiltration and gene signature might be variable because osteosarcoma was high heterogeneity. Most importantly, our research was a retrospective study and the result should be verified in other prospective studies in the future.

Conclusion

In summary, a novel inflammatory response-associated gene signature that related to osteosarcoma prognosis was constructed and validated in our study. Besides that, we provided a new perspective on the correlation between immune infiltration and inflammatory response. Meanwhile, this signature could be regarded as the prognostic biomarker for osteosarcoma and acted as the therapeutic target in clinical studies. Further experimental and clinical studies should be conducted to illuminate the precise mechanisms of inflammatory response in osteosarcoma.