Establishment of immune prognostic signature and analysis of prospective molecular mechanisms in childhood osteosarcoma patients

Supplemental Digital Content is available in the text


Introduction
Osteosarcoma (OS) is the most common solid malignant disease in childhood (an annual incidence of 5.6 cases per million children). OS was rarely diagnosed before age of 5, and the incidence rate reached the peak with age until 10 to 14 years old. The second peak of the age distribution of OS occurred over 50 years old. [1][2][3][4] At present, the treatment of OS is mainly based on resection of primary lesions and combination of multiple The authors report no conflicts of interest.
Supplemental Digital Content is available for this article.
The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request. chemotherapies. The overall survival rate of OS was significantly improved after systematic treatment. [5,6] However, OS is prone to metastasis in the early stage, and about 30% of patients die of lung metastasis even within 5 years after diagnosis. [7,8] In addition, drug resistance worsens the prognosis of OS. [9] Therefore, it is essential to identify early diagnostic markers and therapeutic targets for OS. Based on the encouraging results of immunotherapy and gene therapy in OS, [10,11] the current study aimed to identify new immune targets associated with prognosis.
"Cancer Immunotherapy" was defined as Breakthrough of the Year 2013. [12] The interactions between the immune system and cancer include immune surveillance, immune cell infiltration, and tumor cytolysis. Immunotherapy counteracts the immune escape of tumors by targeting the tumor microenvironment and reactivates the patient's immune system to achieve the goal of recognizing and eliminating tumor cells. [13,14] In pediatric tumors, immunotherapy exhibits less toxicity than chemotherapy and radiation. The aim of the current study was to identify prognostically relevant immune targets for childhood OS.
A large amount of evidence indicates that C-C motif chemokine ligand genes (CCLs) play an important role in the development, progression, and metastasis of cancers, such as gastric cancer, [15] colorectal cancer, [16] breast cancer, [17] and oral squamous cell carcinoma. [18] C-C motif chemokine receptors (CCRs) are the key mediator of inflammation and immune response, and have significant association with multiple cancers. [19,20] Nevertheless, the relationship between CCLs, CCRs and the prognosis of childhood OS is still uncertain. Therefore, the aim of the current study was to identify immunerelated genes in CCLs and CCRs, and to explore the relationship between these genes and prognosis of childhood OS.

Gene expression datasets and immune-related genes
The gene expression quantification data of 88 childhood OS samples were of the HTSeq-FPKM type, and downloaded from the University of California, Santa Cruz Xena (UCSC Xena; http://xena.ucsc.edu/) on December 20, 2019. The corresponding clinical data were downloaded from UCSC Xena at the same time. Immune-related genes were downloaded from the Immunology Database and Analysis Portal (ImmPort) database (https://www.immport.org/). This study was based on data from an open database. Ethics and patient consent are not applicable.

Bioinformatics analysis of immune-related genes in CCLs and CCRs
To learn more about the functionality of immune-related genes in CCLs and CCRs, we use the Database for Annotation, Visualization and Integrated Discovery (DAVID; https://david. ncifcrf.gov/; version 6.8) [21] to perform gene ontology [22] and Kyoto Encyclopedia of Genes and Genomes (KEGG) [23] pathway on these genes. And then, a protein-protein interaction network was constructed using the Search Tool for the Retrieval of Interacting Genes (STRING; https://string-db.org/; version 11) [24] to analyze the correlation between these genes. Finally, a gene-gene interaction network was constructed using Gene-MANIA (http://genemania.org/). [25] 2.3. Survival and correlation analysis According to gene expression level, patients were divided into high-group and low-expression group, and then different genes were analyzed for survival (Kaplan-Meier analysis with log-rank test). According to the results, the correlation between genes with significant significance for childhood OS prognosis was analyzed with the Pearson correlation coefficient in R (http://cran.rproject.org/; version 3.6.0) using corrplot package. Then the prognosis-related genes were analyzed for combined effect survival (Kaplan-Meier analysis with log-rank test), and a nomogram was constructed. The nomogram was constructed by rms package in R according to the expression of prognosisrelated genes.

Prognostic signature construction
According to the survival analysis results of the above genes, the genes significantly related to the prognosis of childhood OS were combined to construct a prognostic model based on gene expression level. The formula for the risk score was as follows: risk score = (expression value of gene A) Â b A + (expression value of gene B) Â b B + . . . (expression value of gene n) Â b n; b meant the regression coefficient. [26] In addition, the survival receiver operating characteristic (ROC) package of R software was used to generate time-dependent ROC curves to test the predictive performance of prognostic features.

Correlation analysis between significant genes expression, risk score, and tumor-infiltrating immune cells
CIBERSORT is a deconvolution algorithm based on gene expression. We downloaded gene expression feature matrix of 22 immune cells from CIBERSORT platform(https://cibersort. stanford.edu/), and then evaluate tumor-infiltrating immune cells of childhood OS samples with CIBERSORT deconvolution algorithm. Finally, the correlation between significant gene expression, risk score of the model, and tumor-infiltrating immune cells was analyzed.

Gene set enrichment analysis (GSEA)
To further explore the biological pathways of enrichment of these significant survival significance genes, we uploaded the relevant data required for these genes to GSEA (version 4.0.3), [27] and then used the databases of c2 and c5 in the Molecular Signature Database (MSigDB) (https://www.gsea-msigdb.org/) [28] for enrichment analysis. The enrichment gene sets in the GSEA that attained a false discovery rate (FDR) of <0.25 and P < .05 were considered statistical significance.

Data collection and collation
Through careful screening, we obtained 85 OS samples with complete clinical data. Secondly, 34 genes were identified as immune-related genes in CCLs and CCRs. Details of these genes including ID, name, synonyms, chromosome, and category (see Supplementary Table S1, http://links.lww.com/MD/F228, which illustrates the details of immune-related genes in CCLs and CCRs).

Bioinformatics analysis of immune-related genes in CCLs and CCRs
In Table 1 and Figure 1, the results of gene ontology and Kyoto Encyclopedia of Genes and Genomes pathway enrichment showed that these genes were enriched in cellular response to tumor necrosis factor, immune response, inflammatory response, positive regulation of T cell migration, mitogen-activated protein kinase (MAPK) cascade, positive regulation of I-kappaB kinase (IKK)/ nuclear factor-kappaB (NF-kB) signaling, Toll-like receptor (TLR) signaling pathway and NF-kB signaling pathway (see Supplementary Fig. S1, http://links.lww.com/MD/F222, which shows the map of TLR signaling pathway and see Supplementary Fig. S2, http://links.lww.com/MD/F223, which shows the map of NF-kB signaling pathway). In Figure 2, both STRING and GeneMANIA results showed that these immunerelated genes in CCLs and CCRs exhibit close associations, such as share protein domains, experimentally determined, coexpressed and gene fusions, etc.

Survival and correlation analysis
Survival analysis showed that CCL5, CCL8, CCR4, and CCR5 were assigned significant statistical significance. High expression of CCL5(Log-rank P = .020,), CCL8 (Log-rank P = .049,), CCR4 (Log-rank P = .016,) and CCR5(Log-rank P = .026,) have better prognosis than low expression. Based on the above 4 genes with significant statistical significance ( Table 2 and Fig. 3A-D), Pearson correlation analysis shows that CCL5, CCL8, and CCR4 have a high correlation, while CCR5 has a relatively low correlation with the other 3 genes. The size, number, and color of the circles all represent the correlation between individual genes (Fig. 3E). The combined effect survival analysis of significant genes showed that group C (corresponding genes were at high expression level) had a better prognosis than group A (corresponding genes were at low expression level) and group B (excluding group A and group C, ie, the expression levels of each gene were not identical). All results are shown in Figure 4. In addition, the results of nomogram also support CCL5, CCL8, CCR4, and CCR5, which were beneficial to the prognosis of childhood OS. The high expression of all genes was at a relatively low point (Fig. 5).

Construction and validation of prognostic signature based on significant genes
In order to more intuitively show the relationship between the expression of significant genes and prognostic risk, and to     . Since the regression coefficients are all negative, in order to make the results more intuitive, the constant is used to add to the end of the risk score formula. The results of Kaplan-Meier curves of high-and low-risk groups showed that the high-risk group was significantly associated with poor prognosis of childhood OS, and the difference was statistically significant (Log-rank P = .024, Fig. 6B). The of the time-dependent ROC curve was 0.54, 0.645, and 0.655 for 1-year, 3-year, and 5-year survival, respectively (Fig. 6C).
3.5. Identified potential relationships between significant genes expression, risk score, and tumor-infiltrating immune cells The results of correlation analysis between the expression of significant genes and tumor-infiltrating immune cells showed that the expression of CCL5, CCL8, CCR4, and CCR5 were all correlated with the level of tumor immune cell infiltration. The expression of CCL5 was significantly correlated with the infiltration levels of macrophages M0, macrophages M1, The results showed that the risk score was positively correlated with the expression of macrophages M0 (P < .001). In addition, the risk score was negatively correlated with the expression of macrophages M1 (P < .001), neutrophils (P = .015), CD8 + T cells (P < .001), and Tregs (P = .001). The results were shown in Figure 7.

GSEA
The analysis of c2 and c5 gene sets showed that CCL5 was mainly enriched in TLR signaling pathway, T to NK cells, Parvin-beta, MAPK family signaling cascade, T cell immune response, and NF-kB pathway (Fig. 8A-D

Discussion
Although immunotherapy has made tremendous progress in adult tumors, less success has been achieved in pediatric tumors. This may be associated with a significant reduction in mutational load in childhood cancers leading to a reduction in the number of neoantigens in immunotherapy. However, the relative low toxicity of immunotherapy makes the immunotherapy of pediatric cancer exhibit infinite charm. [14,29] Fortunately, bone biology is inextricably linked to the immune system, which has led to significant implications for immunotherapy of childhood OS. [30] Therefore, many researchers have carried out a large number of studies on OS immunotherapy, such as denosumab targeting RANK ligand, [31] adding decitabine to upregulate the expression of cancer antigens, [32] T-VEC and anti-PD-1 antibody pembrolizumab combination test, [33,34] etc. However, in the past 3 decades, progress in improving the prognosis of OS patients has been limited by the failure to identify new active agents or to optimize the use of existing drugs. [32] Therefore, it is necessary to identify potential immune prognostic targets for childhood OS.
The purpose of the present study was to identify the relationship between immune-related genes in CCLs, CCRs, and childhood OS patients' prognosis. Survival analysis of immune-related genes in CCLs and CCRs showed that high expression of CCL5, CCL8, CCR4, and CCR5 all represented a good prognostic signal for childhood OS. The biological function of CCL5 in tumor is not clear. On the one hand, the production of CCL5 is related to the induction of appropriate immune response to tumor, but on the other hand, CCL5 was related to the progression and metastasis of tumor. [35,36] It was reported that the activation of CCL5 will raise tumor-associated macrophages into tumor microenvironment. Tumor-associated macrophages interact with tumor cells to secrete multiple factors and promote the release of matrix metalloproteinase, thus activating NF-kB signaling pathway to induce epithelial-mesenchymal transition. [37] In addition, CCL5 blocking the selective depletion of Treg cells in tumor microenvironment can promote the anti-tumor immune response of pancreatic ductal adenocarcinoma patients. [38] But we also note that in esophageal squamous cell carcinoma, CCL5 induces CD8 + T cells to infiltrate tumor cells and is associated with better survival. [39] We also found that CCL8 expression leads to the migration of monocytes and T lymphocytes and plays an antitumor role. And CCL8 can directly inhibit the proliferation of tumor cells or transplant tumor cells into tissues. [40] Di Stasi et al found that T lymphocytes coexpressing CCR4 and targeting CD30 chimeric antigen receptor enhance homing and antitumor activity in Hodgkin's tumor model. [41] In addition, superagonists of CCR5 regulate upstream and downstream events  of antitumor responses by participating in adaptive and innate immunity. [20,42] These findings all effectively corroborate our current findings on CCL5, CCL8, CCR4, and CCR5, however, the current findings have limitations because they were obtained from a single cohort and therefore require an external validation cohort to prove the current results. Tumor-infiltrating immune cells are important for studying the interaction between tumor and immunity, so we explored the correlation between risk score and immune cell infiltration. The results showed that patients with high-risk scores had high levels of macrophages Mo infiltration, whereas patients with low-risk scores had high levels of macrophages M1, neutrophils, CD8 + T cells, and Tregs infiltration. Similar to our results, Yang et al found higher levels of macrophages M0 are a high-risk immune subtype for cancer recurrence in the digestive system. [43] Xiong et al found macrophages M1 has antitumor effect in colorectal cancer. [44] Granot et al found that tumor-entrained neutrophils inhibited lung dissemination before transplantation. [45] T cells are an important part of immune cells, which have antitumor effect, especially CD8 + and CD4 + T cells. [46] Previous studies have found that the presence of Tregs is associated with a positive prognosis in head-and-neck and gastric cancer. [47,48]  Interestingly, some subtypes of tumor-infiltrating immune cells (such as macrophages M2, B cells, and CD4 + T cells) infiltrate in different tumors, and the difference is given statistical significance, but we did not find any association in the current study, even the subtype infiltration results and our current study results have the opposite situation, which may be related to the difference between different tumors and the functional heterogeneity of immune cell subtypes in tumor progression. Comprehensive analysis of functional annotations and GSEA results showed that CCL5, CCL8, and CCR5 were significantly enriched in immunity, inflammation, TLR signaling pathway, MAPK family signaling cascade, and NF-kB pathway. The inseparable relationship between immunity, inflammation, and cancer has been widely reported. [49,50] Previous studies have shown that TLRs are negative regulators of cancer and activation of TLRs can lead to activation of MAPKs and NF-kB. [51] However, it has also been found that some TLRs are found in many tumors and produce the effect of tumor cell proliferation and resistance to chemotherapy, but depend on the TLR and tumor type. [52] The role of the NF-kB pathway in cancer is a double-edged sword. On the one hand, the activation of NFkappa B can cause immune defense and target transformed cells; on the other hand, NF-kB is constitutively activated in many types of cancer and can play a variety of cancer-promoting roles. [53] Schulze-Osthoff et al found that complex crosstalk effects occur between the NF-kB pathway and members of the MAPK family. [54] In addition, the results of GSEA showed that CCR4 with low expression was enriched at the initiation and elongation of eukaryotic translation. These pathways are also closely related to cancer. [55,56] We speculate that CCL5, CCL8, and CCR5 may inhibit the progression of childhood OS through the complex regulation of TLR signaling pathway, MAPK family signaling cascade, and NF-kB pathway. And CCR4 may influence the prognosis of childhood OS by regulating eukaryotic translation. However, these findings need further experimental verification.
However, there are some limitations at present. First, our analysis sample size is relatively small, which may lead to falsenegative results. Secondly, the current study is retrospective and needs further prospective study. Finally, the data of the current study come from a single cohort, and the combination of data from multiple cohorts will make the results more convincing.
Although the current study has some limitations, our study provides a new idea for immunotherapy of childhood OS. We identified immune-related genes in CCLs and CCRs that influence childhood OS prognosis, and secondly, we identified potential molecular mechanisms that influence childhood OS prognosis through functional annotation and GSEA. The selected genes and molecular mechanisms can be prioritized for experimental studies to demonstrate their clinical value.

Conclusion
CCL5, CCL8, CCR4, and CCR5 are potential biomarkers for the prognosis of childhood OS and have great clinical prospects. In addition, CCL5, CCL8, and CCR5 may affect the prognosis of childhood OS through complex regulation among TLR signaling pathway, MAPK family signaling cascade, and NF-kB pathway, while CCR4 may affect the prognosis of childhood OS by affecting eukaryotic translation.