Profiles of immune cell infiltration and immune-related genes in the tumor microenvironment of osteosarcoma

1Osteosarcoma (OS), a common primary bone malignancy, tends to metastasize. Between 1973 and 2012, the overall incidence rate of OS was 4.5 per million in the United States. According to statistics reported in the United States from 2007 to 2013, the five-year relative survival rates of OS patients were 69.8% and 65.5% in the ages from birth to 14 years and 15 to 19 years, respectively. Currently, new OS cases are administered neoadjuvant chemotherapy and surgery to remove the primary and overt metastatic tumors, with postoperative adjuvant chemotherapy; this has resulted in increased overall survival in OS. However, drug resistance has worsened patient prognosis. Therefore, it is important to develop additional efficient therapeutics to improve survival in OS.

2As an emerging treatment, immunotherapy has shown promising results for some cancers, including hepatocellular carcinoma and breast cancer. The tumor microenvironment (TME), a mixture that consists of mesenchymal cells, tumor-infiltrating immune cells (TIICs), endothelial cells, extracellular matrix molecules and inflammatory mediators, provides all metabolites and factors for controlling proliferation, dissemination, dormancy, and drug resistance in OS cells. It was suggested that the TME plays a critical role in OS development. In the TME, TIICs constitute the major type of non-tumor components reported to be valuable for prognostic assessment in OS. Thus, improving immunotherapy efficacy in OS by systematically assessing the TME’s immune properties and determining TIIC distribution and functions is of prime importance.

3An algorithm has been developed to predict the levels of TIICs using gene expression data from the cancer genome atlas (TCGA) (https://portal.gdc.cancer.gov/), and immune score could be calculated for predicting immune cell infiltration, by analyzing a specific gene expression signature of TIICs. Recently, several studies have applied this algorithm to glioblastoma multiforme and clear cell renal cell carcinoma, showing the feasibility of such big-data based algorithms, although the immune scores of OS cases from the TCGA database have not been investigated in detail. Moreover, Cell type Identification By Estimating Relative Subsets Of RNA Transcripts (CIBERSORT), can use the deconvolution technique to assess the levels of 22 TIICs in large amounts of heterogeneous samples. CIBERSORT has been successfully applied for identifying TIIC landscapes and their associations with prognosis in colorectal, gastric and breast cancer.

4To increase immunotherapy efficacy, determining immune-associated prognostic biomarkers is especially pivotal. Here, we calculated immune score of OS cohorts in the TCGA database by taking advantage of the algorithm, known as ESTIMATE, retrieved immune-associated differentially expressed genes (DEGs) in OS, and built a predictive risk model to estimate patient outcome. Importantly, we also evaluated the associations of the immune-related risk score with the levels of TIICs and immune pathways.

5We first determined immune score of the normalized matrix data of 85 OS samples with complete clinical data by applying the ESTIMATE algorithm. Subsequently, the OS cases were assigned to the high and low immune score groups respectively, according to the median value of immune scores. Kaplan-Meier curves revealed that the high immune score was significantly associated with improved outcome (P=0.002) (Figure 1). The five-year survival rates in cases with high and low immune score were 82.1% and 48.5%, respectively.

6Of all OS samples, 38 and 43 with low and high immune score, respectively, were eligible based on CIBERSORT P<0.05. The two most common TIICs in OS tissues were macrophages and T lymphocytes, which accounted for more than 80% of all TIICs. Specifically, the proportions of naive B cells (Z=-3.014, P=0.003) and M0 macrophages (Z=-3.095, P=0.002) were significantly lower in high immune score tissues compared with the low immune score group, while the proportions of M1 macrophages (Z=-3.047, P=0.002), M2 macrophages (Z=-3.785, P<0.001) and resting dendritic cells (Z=-2.251, P=0.024) were significantly higher (Figure 2). Furthermore, the ratio of M1 macrophages to total polarized macrophages (M1 and M2) showed no significant difference between high and low immune score tissues (Z=-1.427, P=0.154). Correlations among the 22 TIICs ranged from weak to moderate. Obviously, M0 macrophages showed highly negative correlations with M1 and M2 macrophages (Figure 3).

7Firstly, we compared 42 low and 43 high immune score OS samples of the normalized matrix data. Compared with the low immune score group, there were 607 upregulated and 459 downregulated DEGs in the high immune score group (Figure 4A). Subsequently, we identified immune-related DEGs. Compared with the low immune score group, there were 177 upregulated DEGs and 14 downregulated immune-related DEGs in high immune score specimens (Figure 4B).

8Univariable Cox regression analysis revealed 34 immune-related genes which were significantly associated with improved outcome (P<0.05) (Table 1). To assess multicollinearity among different covariates in the model, we excluded variables with variance inflation factor (VIF) >5 (Table 2). A total of 15 genes were excluded, while 19 were included in multivariate Cox regression analysis. Finally, a minimum Akaike information criterion (AIC) value of 210.64 was estimated by the R software to develop the optimal multivariate Cox regression model. The predictive model was then built with three genes, including peroxisome proliferator activated receptor gamma (PPARG), immunoglobulin heavy constant gamma 3 (IGHG3), and pyruvate dehydrogenase kinase 1 (PDK1). Based on relative coefficients in multivariable Cox regression analysis, the following formula was obtained for the risk score: (-0.7728 * PPARG expression level) + (-0.3620 * IGHG3 expression level) + (0.4210 * PDK1 expression level). The median value of risk scores was used as a cutoff to divide samples into two groups (Figure 5A). As shown in Figure 5B, the number of deaths was significantly higher while overall survival was shorter in high-risk cases compared with the low-risk group. Furthermore, in comparison with the low-risk group, high-risk cases had lower expression levels of IGHG3 and PPARG, and higher PDK1 amounts (Figure 5C). High risk score was significantly associated with poor outcome (P<0.01), revealing this score as a good predictive tool (Figure 6). The five-year survival rates of high and low risk score cases were 45.4% and 83.8%, respectively. Moreover, the areas under the curve (AUC) for the risk model in predicting 1, 3 and 5-year survival were 0.634, 0.781, and 0.809, respectively (Figure 7A–7C).

9As shown in Table 3, the risk score was positively correlated with the proportions of memory B cells (r=0.252, P<0.05), M0 macrophages (r=0.305, P<0.01) and resting dendritic cells (r=0.246, P<0.05), and negatively correlated with those of gamma delta T cells (r=-0.245, P<0.05) and M2 macrophages (r=-0.244, P<0.05).

10Two immune gene sets, including M19817 (immune response) and M13664 (immune system process), were retrieved from the Molecular Signatures Database v4.0 (http://software.broadinstitute.org/gsea/downloads.jsp). As shown in Figure 8A, 8B by gene set enrichment analysis (GSEA), both immune response and immune system process gene sets were significantly enriched in the low-risk group (P<0.001).

11As is calculated by the ESTIMATE algorithm, it is well admitted that elevated immune score is significantly correlated with poor prognosis in clear cell renal cell cancer patients. This aroused our interest in exploring a potential association of immune score with survival in OS patients. To the best of our knowledge, this is the first study building an immune-related risk model to predict outcome in patients with OS by mining the TCGA database. In the present study, we firstly evaluated approximate proportions of TIICs in OS TME, calculating immune score by applying the ESTIMATE algorithm. Importantly, a high correlation was found between the immune score and overall survival in OS patients. It has been demonstrated that TIICs are significantly relevant to the progression and prognosis of OS. In order to explore specific differences in the proportions of TIICs, OS cases were assigned to the high and low immune score groups. Then, the types of TIICs were assessed in both groups of OS tissues with CIBERSOTR. Subsequently, immune-related DEGs were screened between high and low immune score OS tissues, and an optimal immune-related risk model was built by univariate and multivariate Cox regression analyses. In this model, high risk score, calculated by the expression levels of three immune-related DEGs, was associated with poor outcome. Moreover, of these three genes, PDK1 overlapped with the gene signatures found on the CIBERSORT platform, which implied that immune-related risk score and the proportions of TIICs may be somehow associated. Fortunately, correlations of the risk score with the proportions of 5 TIICs were determined in this study. Finally, we applied GSEA to assess the associations of immune pathways with the determined risk score.

12It has been reported that tumor-associated macrophages (TAMs) and T-lymphocytes are the main components of the immune environment in OS, in agreement with the above results. We found that OS cases with elevated immune cell infiltration in the microenvironment had better prognosis. Compared with low immune score cases, patients with high immune score showed markedly decreased levels of M0 macrophages and significantly increased amounts of M1 and M2 macrophages, especially M2 macrophages. Furthermore, the risk score was negatively correlated with the proportions of M1 and M2 macrophages, and positively correlated with the proportion of M0 macrophages, suggesting that the polarization level of M0 to M1 or M2 macrophages may be associated with improved outcome in OS patients. In preclinical models of OS, M2-TAMs are associated with increased tumor growth, metastatic dissemination and vascularization. Excitingly, contrary to findings reported for other solid tumors, such as gastric cancer, lung adenocarcinoma, and colon cancer, studies by Anne Gomez-Brouchet et al. indicated that the presence of CD163-positive M2-polarized macrophages is essential for inhibiting OS progression, which represents an important discovery. However, Buddingh et al. described TAMs in OS as a heterogeneous cell population with both M2 pro-tumor and M1 anti-tumor characteristics. Interestingly, Cristiana Guiducci et al. reported the plasticity of TAMs, with CpG combined with anti-interleukin-10 receptor antibodies readily switching them from M2 to M1. Recently, switching TAMs from M2 to M1 has been suggested for developing novel treatments. It has been demonstrated that all-trans retinoic acid suppresses pulmonary metastasis of OS cells by inhibiting M2-like TAMs, which may lead to clinical application in metastatic OS. Based on the above findings, we further analyzed the balance between M1 and M2 macrophages in OS tissues. The ratio of M1 macrophages to total polarized macrophages (M1 and M2) was only slightly elevated in the high immune score group (6.040%) compared with the low immune score group (4.741%), but this difference was not statistically significant. Therefore, we speculated that small changes in the balance of polarized macrophages may be an important factor affecting the prognosis of OS patients.

13In this study, the proportion of resting dendritic cells in high immune score tissues was significantly higher than that of the low score group. Moreover, the risk score was positively correlated with the proportion of resting dendritic cells, implying that the activation level of dendritic cells may be associated with improved outcome in OS patients. Masanori Kawano et al. reported that combining agonist anti-glucocorticoid-induced tumor necrosis factor receptor (GITR) antibodies with tumor lysate-pulsed dendritic cells reduces the amounts of immunosuppressive cytokines in OS tissues as well as serum. Furthermore, it has been confirmed that pulsing of dendritic cells with LM8 cell lysate, derived from OS, efficiently enhances CD4+ and CD8+ T cell proliferation and decreases serum interleukin-4. While assessing synergistic effects with chemotherapy, it was found that combining doxorubicin, which induces immunogenic cell death, with resting dendritic cells boosts systemic immune reactions, leading to OS inhibition in mouse models. Currently, the dendritic cell-based vaccine, a form of active specific immunotherapy, is considered to confer a possible overall survival advantage in children with cancer, similar to findings in adults. The most promising clinical effects were observed in cases with limited disease or complete response, in whom the complete response state could be maintained upon dendritic cell-vaccination, preventing the tumor from recurring. Conversely, dendritic cell-vaccination shows reduced effects in cases with progressive disease or elevated residual tumor load, most likely for the extremely high immunosuppressive burden of malignant cells, with insufficient time to produce appropriate antitumor immune reactions.

14We screened differentially expressed immune-associated DEGs, and performed univariable and multivariable Cox analyses to generate a risk model for predicting the prognosis of OS patients. Three DEGs were used to construct the model. PPARG and IGHG3 are two protective immune-related DEGs, while PDK1 is a risk immune-related DEG. The AUC values for the risk model in predicting 1, 3 and 5-year survival were 0.634, 0.781, and 0.809, which indicates a good capability for predicting survival in OS patients of the three-gene combination. Cases were assigned to the high- and low-risk groups based on the median risk score. More than 40% of high-risk cases died within three years of diagnosis, while less than 14% died in the low-risk group. Therefore, extremely frequent follow-up and more aggressive treatments should be applied in the high-risk group. Finally, GSEA further confirmed the close connection of the risk signature with immune pathways. As shown above, immune response and immune system process gene sets were significantly enriched in the low-risk group, which suggests that immunosuppression may exist in high-risk OS patients, and is associated with poor outcome.

15PPARG is a ligand-activated transcription factor, belonging to the nuclear hormone receptor family. Accumulating evidence confirms that activation of PPARG could confer inhibitory effects on tumors. A meta-analysis identified an association of PPARG c.1347 C > T polymorphism with elevated risk of developing malignancies such as glioblastoma and esophageal cancer. In clear cell renal cell cancer, PPARG suppresses cell migration and proliferation and induces apoptosis by inhibiting SIX homeobox 2. Sabatino et al. reported that ring finger domains 1 regulates PPARG negatively and is associated with higher clonogenic, proliferative and migratory potential in colorectal cancer. These findings suggest that PPARG might also represent a tumor suppressor gene in OS. In addition, it has been demonstrated that PPARG is a direct target of miR-27a by luciferase reporter assays. Analysis of differentially expressed miRNAs and their target genes in OS samples showed that miR-324-5p targets PPARG, which needs to be verified in future experiments. Furthermore, PPARG expression is independently associated with prolonged survival in colorectal cancer, as demonstrated in two separate prospective cohorts, corroborating the findings of the present study in OS patients.

16IGHG3 is a member of the immunoglobulin G family. Several studies indicated that IGHG3 is overexpressed in multiple cancer types, such as prostate, breast, and lung cancers, which can differentiate tumor from normal tissues. Hsu et al. suggested that IGHG3 expression is tightly associated with improved outcome in breast cancer, which is consistent with our findings in OS. The relationship between IGHG3 expression and OS progression deserves further attention.

17PDK1 is a hypoxia-inducible factor-1α target antagonizing pyruvate dehydrogenase (PDH), a pivotal rate-limiting enzyme of the tricarboxylic acid cycle. In hypoxia, pyruvate transformation into acetyl-CoA is suppressed due to PDK1-dependent inhibition of PDH, reducing the amounts of glucose-derived pyruvate entering the tricarboxylic acid cycle. It was reported that PDK1 downregulation in metastatic breast cancer greatly alters the tumor cell capability to utilize glucose as an energy source for the mitochondria under hypoxic or limited glucose conditions. Moreover, for coping mechanisms against stress during metastasis, PDK1 mediates the adaptability of breast cancer cells metastasizing to the liver. Meanwhile, Liu et al. confirmed that downregulation of PDK1 could inhibit migration and metastasis in human breast cancer cells. In addition, overexpression of PDK1 has been reported in multiple myeloma, acute myeloid leukemia, breast cancer and OS. Li et al. demonstrated that overexpression of PDK1 promotes the proliferation of OS cells. More importantly, PDK1 is a direct target of miR-379, which functions as a tumor-inhibiting miRNA by targeting PDK1 in OS. Two novel PDK1 inhibitors were shown to concentration-dependently reduce the phosphorylation of the pyruvate dehydrogenase complex in MG-63 OS cells, whose proliferation was inhibited as a result.

18There were several limitations in this study. First, the number of OS tissue samples in the TCGA cohort was relatively small, which could lead to some bias. Secondly, the landscape differences of TIICs and immune-related DEGs between tumor and normal samples were not analyzed in the TCGA cohort because noncancerous samples were not included, and sampling normal bone tissue is subject to restriction in clinic to some extent. Thus, the present findings could only be applied to predict the prognosis of OS patients with definite diagnosis. Thirdly, because the clinical information in the database does not include the tumor stages of OS tissues, we could not perform subgroup analysis based on tumor stage. However, gene expression in OS tissues of different tumor stages may be different, and further research is warranted to address this issue. Finally, the present study performed no external validation based on other available databases; therefore, the current conclusion requires validation in future experiments.

19In summary, according to the ESTIMATE algorithm-based immune score that was significantly correlated with improved outcome, 22 TIICs in OS TME were assessed for their levels. Then, a list of immune-related DEGs was extracted, and three such genes (PPARG, IGHG3, and PDK1) were included in a predictive risk model, which could assist clinicians in assessing the prognosis of OS patients and selecting appropriate targets for immunotherapy.

20The gene expression quantification data of 88 OS samples were of the HTSeq-FPKM type, and downloaded from TCGA (version 18.0) on August 17, 2019. Updated clinical data related to these OS samples, such as age, gender, race, overall survival time and vital status, were also downloaded from the TARGET database (https://ocg.cancer.gov/programs/target) on August 17, 2019. In order to screen out the matrix data of mRNAs with gene properties, the gene expression profiles were compared with the human genome annotation GTF file, which was downloaded from the GenCode platform (https://www.gencodegenes.org/). Subsequently, the matrix data of the gene expression values was organized by the Perl software (version 5.24.3) (https://www.perl.org/).

21The matrix data of gene expression amounts were normalized with the limma package of the R software (version 3.5.2). Then, immune score was calculated by applying the ESTIMATE algorithm to the matrix data. Furthermore, the OS cases were assigned to high- and low immune score groups based on the median value of immune scores, to identify a possible association of immune score with overall survival.

22TIICs in OS samples from the TCGA cohort were assessed by applying the CIBERSORT deconvolution algorithm. The gene expression signature matrix of 22 TIICs was obtained from the CIBERSORT platform (https://cibersortx.stanford.edu/). The matrix data of gene expression levels were compared with those of the signature matrix of 22 TIICs from the CIBERSORT platform to generate a proportion matrix for the 22 TIICs in OS tissues of the high and low immune score groups using support vector regression. By Monte Carlo sampling, the algorithm derives a P-value for the deconvolution of each sample, offering a measure of confidence for the obtained data. The results of the inferred proportions of TIICs assessed by CIBERSORT were considered to be accurate at a threshold of P<0.05. Therefore, only samples with a CIBERSORT P<0.05 were deemed qualified for further analysis. Moreover, the number of permutations of the default signature matrix was set to 100.

23To further examine the DEGs between low and high immune score OS samples of the TCGA cohort, the normalized matrix data of gene expression levels were analyzed with the limma package of the R software. Log fold change >1 or <-1 and adjusted P<0.05 were set as cut-offs for filtering DEGs. In addition, immune-associated genes were retrieved from the ImmPort platform (https://www.immport.org/), and used to identify immune-related DEGs for constructing predictive risk model.

24Only samples with complete clinical data were included in survival analysis, and the logrank test was performed for comparing Kaplan-Meier curves between groups. Differences and correlations among TIICs were analyzed with the vioplot (https://cran.r-project.org/web/packages/vioplot/index.html) and corrplot (https://cran.r-project.org/web/packages/corrplot/index.html) packages of the R software. The differential proportions of the 22 TIICs in the TCGA cohort were evaluated by the Wilcoxon rank-sum test. A heat map was generated using the pheatmap package of the R software (https://cran.r-project.org/web/packages/pheatmap/index.html). The Cox proportional-hazards model was used for analyzing associations of the levels of immune-related DEGs with overall survival. Collinearity diagnostics was performed with the SPSS software (version 24.0) (SPSS, USA). The multicollinearity of each variable was estimated by calculating the VIF. A variable with VIF>5 was considered to show high collinearity, and would be excluded from multivariable Cox regression analysis. The optimal multivariable Cox regression model was selected according to the lowest AIC. Coefficients for each covariate were determined by the multivariable Cox regression model, and a total risk score was calculated. The specificity and sensitivity of survival prediction according to the determined risk score were obtained by time-dependent receiver operating characteristic (ROC) curves, with AUC values quantified with the survivalROC package (https://cran.r-project.org/web/packages/survivalROC/index.html). Next, the association of the immune-related risk score with TIIC levels was assessed by spearman rank correlation using the R software. GSEA (http://www.broadinstitute.org/gsea/index.jsp) was carried out to evaluate associations of immune pathways with the immune-related risk score using the GSEA software (version 4.0.1). P< 0.05 indicated statistical significance.