ATG16L1 is a Potential Prognostic Biomarker and Immune Signature for Osteosarcoma: A Study Based on Bulk RNA and Single-Cell RNA-Sequencing

Introduction

Osteosarcoma is a common bone tumor in children and adolescents, accounting for 15% of solid extracranial cancers in this age group.^1,2 Osteosarcoma typically originates from the long bones of the limbs near the metaphyseal growth plate and shows a high propensity for metastasizing to the lungs. Use of chemotherapy in combination with surgical resection has helped improve the 5-year survival rate of patients with osteosarcoma.³ Nevertheless, in the past few decades, there has been no major improvement in the overall survival rate of these patients (5-year overall survival rate of patients with lung metastasis: 20%).^4,5 Patients with metastasis and recurrence show poor response to chemotherapy and radiotherapy.

An increasing number of studies have indicated the potential value of autophagy in clinical treatment.^6–8 Use of chemotherapy activates autophagy in tumor cells, which accelerates development of resistance to drug cytotoxicity and renders cancer cells insensitive to treatment.⁹ Use of autophagy inhibitors in combination with anti-tumor drugs can have greater anti-tumor effects by Increasing the sensitivity of osteosarcoma cells to drugs.¹⁰

Autophagy is a dynamic homeostatic mechanism that helps maintain the stability of the intracellular environment by delivering unnecessary organelles to the lysosome for degradation.¹¹ It is a multistep regulatory event in response to cancer treatment. Autophagy can play a dual role in cancer cells. Autonomous cellular autophagy can inhibit tumorigenesis while it can promote tumorigenesis in established tumors, especially those in metastatic stages. Early deletion of autophagy-related genes (ATGs) seems to mainly inhibit the lymphatic lineage and lead to bone marrow proliferation. Autophagy is stimulated and required for peripheral T cell and effector survival and CD8⁺ memory T cell formation in naive T cells. CD8⁺ T cells play an important role in anti-tumor immunity, and autophagy is an autonomous regulator of CD8⁺ T cell metabolism.¹²

The development of high-throughput sequencing technology has allowed transcriptome analysis at the molecular level.¹³ Bulk RNA sequencing can help detect RNA transcripts in the tissue or mass of cells. Advances in scRNA-seq have facilitated in-depth characterization of the composition and heterogeneity of the tumor microenvironment (TME) at the single-cell level, promoting the research and development of new biomarkers and drugs. Recently, applications of scRNA-seq in the context of bone tumors and autophagy have been reported.^14–19

In the present study, we explained how ATGs affect osteosarcoma by mediating CD8⁺ T cells through the traditional bulk RNA sequencing and scRNA-seq.

Materials and Methods

Acquisition of Gene Expression Datasets

The gene expression datasets GSE21257 and GSE152048 were downloaded from the Gene Expression Omnibus database (https://www.ncbi.nlm.nih.gov/geo/). GSE21257 includes microarrays of pre-chemotherapy biopsy of patients with osteosarcoma that metastasized within 5 years (n=34) and in patients with osteosarcoma that did not metastasize within 5 years (n=19). The data in the profiles were obtained from the GPL10295 platform. Our study was approved by the Medical Ethics Committee of the First Affiliated Hospital of the Guangxi Medical University (No. 2021KY-E-230).

Consensus Clustering Analysis

The GSE21257 dataset was clustered using the “ConsensusClusterPlus” R package (version 1.50.0). The different gene expression patterns were subsequently confirmed through principal component analysis (PCA). Then, the “survival” R package (version 3.2.3) was applied to analyze the overall survival of patients.²⁰

CIBERSORT Analysis

The CIBERSORT algorithm is a general method to predict immune infiltration. The matrix (immune cell) was filtered based on CIBERSORT P<0.05 and then the relative expressions of 22 types of immune cells were identified in 53 osteosarcoma samples in GSE21257. “Pheatmap” and “vioplot” packages were used to visualize the corresponding results.

Identification of Differentially Expressed Genes (DEGs)

DEGs were identified using the “Limma” package (version 3.42.0)²¹ on the basis of the expression levels in the osteosarcoma patients. The selection criteria of DEGs were as follows: log2 fold change (FC) was greater than 2 or less than −2 and p value less than 0.05. The analysis results were visualized using “Pheatmap” package.

Protein-Protein Interaction (PPI) Network Analysis

STRING (https://string-db.org/) was conducted to predict interactions between functional proteins and establish a PPI network. The network was conducted with Cytoscape software (version 3.8.0).

Prognostic Signature Generation and Prediction

LASSO regression analysis was performed to set up a risk signature using “glmnet” package (version 4.1–3). Univariate and multivariate Cox regression analyses were used to assess the risk score by “survival” package (3.2–11). The dataset was grouped into two parts according to the median-risk score: high-risk group and low-risk group. The prognosis of the two groups were compared using Kaplan–Meier analysis and Log rank test. Then, 5-year receiver-operating characteristic (ROC) curves were established to assess the predictive accuracy of risk signatures.

Quantitative Reverse Transcription Polymerase Chain Reaction (qRT-PCR)

The human osteoblasts hFOB1.19 and human osteosarcoma cells U2, MG63, 143B were obtained from the FuHeng Cell Center (FuHeng Cell Center, Shanghai, China). The hFOB1.19 was considered as the control group. All cells were grown in Dulbecco’s modified Eagle’s medium (DMEM, Gibco, Carlsbad, CA, USA) in the presence of 10% fetal bovine serum (FBS, Gibco) and 1% streptomycin/ penicillin (Solarbio, Beijing, China) in a humidified atmosphere containing 5% CO2. The osteosarcoma cells were maintained at 37 °C, and hFOB1.19 cells were maintained at 33.5 °C. The RNA was extracted from the induced samples using Hipure Total RNA Mini kit (Magen, China). The extracted RNA was reverse transcribed into cDNA using a cDNA synthesis kit (Takara, Japan). qRT-PCR was performed to analyze the gene expression levels using SYBR Green (LightCycler^® 480 SYBR Green I Master, Roche, UK). Glyceraldehyde 3-phosphate dehydrogenase (GADPH, Abcam, USA) was used as the internal control. The experiments were performed in triplicate. The primer sequences are listed in Table 1.

Table 1 The Primers of ATG16L1 and GAPDH

Analysis Based on scRNA-Seq Transcriptome Datasets

The scRNA-seq transcriptome datasets of GSE152048 were obtained from the NCBI web site (https://www.ncbi.nlm.nih.gov/geo). Seven primary osteosarcoma samples were selected for analysis. Then the “Seurat” package (version 3.1.1) was used to perform the quality control of single-cell data, which was merged via the “Merge” function.^22,23 Low-quality cells were screened out in accordance with the criteria of gene number < 300 and > 5000 and those with a mitochondrial gene number > 10%. In addition, the batch effect was eliminated using “Harmony” package (version 1.0).²⁴ Then, the “FindClusters” function of Seurat package with a resolution parameter of 0.05 was used to achieve cell clusters. T cells were acquired using the “Subset” function of the “Seurat” package for clustering, and marker genes were applied to annotate cell clusters including CD8⁺ T cells. Then, CD8⁺ T cells were clustered into two group based on the absence or presence of ATG16L1 expression namely, ATG16L1⁺ CD8⁺ T cells and ATG16L1⁻CD8⁺ T cells, to explore the relation between ATG16L1 and cytotoxic effector molecules. Gene set variation analyses (GSVA) was used to estimate the variation of pathways using the “GSVA” package (version,1.34.0).²⁵

Statistical Analysis

All statistical analyses were carried out using R (version 3.6.3). Survival analysis was evaluated using Cox regression analysis. The risk ratio (hazard ratio [HR]) and 95% confidence interval (CI) were calculated to identify the genes associated with overall survival. Survival outcomes were computed using Kaplan–Meier method and between-group differences were assessed using the Log rank test.Two-sided P values < 0.05 were considered indicative of statistical significance.

Results

Consensus- Clustering Analysis

In the consensus- clustering analysis, according to the similarity of the gene expression, it was the optimal choice to divide the osteosarcoma samples into two clusters (Figure 1A–C). Kaplan–Meier survival analysis showed that the overall survival rate in cluster 1 was obviously higher than that in cluster 2 (Figure 1D). PCA confirmed the results (Figure 1E). This suggested that the data can be divided into two groups.

Figure 1 Classification of osteosarcoma patients. (A) Consensus CDF for k = 2–10. (B) Relative change in area under CDF curve (k = 2–10). (C) Osteosarcoma patients were separated into two clusters (k = 2). (D) Patients in cluster 1 group had better overall survival. (E) The clusters were confirmed by PCA. Pink represents cluster 1 and cyan represents cluster 2. P < 0.05 was considered indicative of significant difference. Abbreviations: CDF, Cumulative distributive function; PCA, principal component analysis.

Immune Cell Fraction Estimation

The immune cell fraction was estimated in all samples (Figure 2A–D). Immune cells, except eosinophils, were counted and exhibited in the two clusters. Macrophages M0, macrophages M2, and mast resting cells accounted for the largest components of immune cells. The fraction of T cells CD8⁺ was apparently higher in cluster 1.

Figure 2 Results of immunoinfiltration analysis showing the immune cells expressed in the osteosarcoma samples. (A) Histogram showing the expression rate of the immune cells in the 53 samples. (B) Heatmap showing the comparison of the immune cells between the two clusters. (C) Identification of differentially infiltrated tumor-infiltrating immune cells between the two clusters via CIBERSORT algorithm. Blue represents cluster 1 and red represents cluster 2. (D) Pearson correlation among the immune cells. P < 0.05 was considered indicative of significant difference.

Construction and Evaluation of DEGs

DEGs (Figure 3A) were tested to construct a prognostic prediction signature for osteosarcoma patients, and then the interaction network among the DEGs was analyzed (Figure 3B–C). The Pearson correlation among the DEGs was used (Figure 3D). Finally, the following 9 DEGs were obtained: PTPN2, PIK3R4, PARK2, LRRK2, GPR65, C9orf72, ATP6AP2, ATG16L1, and APP.

Figure 3 The interaction and correlation among the DEGs. (A) Heatmap of the DEGs between cluster 1 and cluster 2. (B) PPI network reduced the DEGs following STRING database. (C) Heatmap showing the DEGs between cluster 1 and cluster 2. (D) Pearson correlation among the DEGs. *P < 0.05, **P < 0.01, ***P < 0.001. P < 0.05 was considered indicative of significant difference. Abbreviation: DEGs, differentially expressed genes.

LASSO Analysis of the Hub Genes

Nine genes were selected for lasso regression analysis to clarify the internal relationship. Finally, a risk signature of five genes was found (Figure 4A–B). Kaplan–Meier survival analysis showed that the low-risk group was associated with significantly greater osteosarcoma survival rate than the high-risk group (Figure 4C; P= 3.079×10⁻²). The area under the curve (AUC) of this risk signature was 0.749 (Figure 4D), which indicated the good sensitivity and specificity of the risk signature. Heatmaps comprising clinical information showed that in the high-risk group, the expression levels of the ATG16L1, PTPN2, PIK3R4, and APP were higher and the expression level of C9orf72 was lower (Figure 4E). Moreover, high -risk group had a larger proportion of metastatic disease. Huvos grade, gender, age, and survival status showed diffused distribution in the two groups.

Figure 4 Prognosis model associated with genes. (A and B) Univariate Cox regression analysis identified the remaining 5 DEGs with LASSO coefficient profiles. (A) The partial likelihood deviance is shown against Log Lambda. (B) LASSO coefficients profiles of DEGs. (C) Kaplan-Meier survival analysis in low- or high-risk osteosarcoma patients. (D) Time-dependent ROC curve analysis indicating the sensitivity and specificity for predicting (AUC=0.749). (E) Heatmap summarizing the clinical information of low- or high-risk osteosarcoma patients. P < 0.05 was considered indicative of significant difference. Abbreviations: DEGs, differentially expressed genes; LASSO, least absolute shrinkage and selection operator; ROC, receiver-operating characteristic; AUC, area under the curve.

Risk Score and Other Clinicopathological Characteristics

Only patients with complete survival data and clinicopathological characteristics were included. Finally, “age,” “gender,” “grade,” and the risk score were incorporated into the univariate and multivariate cox regression analyses (Figure 5A–B). Grade and risk score showed a significant correlation with overall survival (HR: 0.482; CI: 0.270–0.861; P = 0.014 and HR: 2.255; CI: 1.231–4.130; P = 0.008, respectively).

Figure 5 Cox regression analyses of the risk score and clinicopathological parameters to select factors significantly associated with overall survival. (A) Forest plot of univariate Cox regression analysis in osteosarcoma patients. (B) Forest plot of multivariate Cox regression analysis in osteosarcoma patients.

Functional Analysis of ATG16L1

Survival analysis showed that only the ATG16L1 gene had significantly improved in the high-risk group (Figure 6A–E; p = 3.40×10⁻²). The AUC of the ATG16L1 gene’s ROC curve was 0.684 (Figure 6F).

Figure 6 Construction of DEGs and assessment of its prognostic value. (A–E) Kaplan-Meier survival curve showing the survival outcomes in the high-risk and low-risk groups of osteosarcoma patients based on the DEGs. (F) Time-dependent ROC curve analysis (AUC=0.684). (G) Results of qRT-PCR. The expressions of ATG16L1 in three osteosarcoma cell lines were significantly lower than that in hFOB 1.19. *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001. p < 0.05 was considered indicative of significant difference. Abbreviations: DEGs, differentially expressed genes; ROC, receiver-operating characteristic; AUC, area under the curve; qRT-PCR, quantitative reverse transcription polymerase chain reaction.

The immune cell fraction (Figure 7) was estimated between the group with expression of the ATG16L1 gene (ATG16L1⁺ group, red) and the group without expression of the ATG16L1 gene (ATG16L1⁻ group, green). CD8⁺ T cells, monocytes, and helper T cells were higher in the ATG16L1⁻ group, whereas macrophage M0 and memory resting CD4 T cells were higher in the ATG16L1⁺ group. CD8⁺ T cells exhibit potent cytotoxic activity against cancer cells.²⁶

Figure 7 Violin plot for comparison of the immune cell fraction difference between the group with expression of the ATG16L1 gene (red) and the group without expression of the ATG16L1 gene (green). P < 0.05 was considered indicative of significant difference.

qRT-PCR

To confirm the gene expression in osteosarcoma, we conducted in-vitro experiments. The results of qRT-PCR revealed that the expression of ATG16L1 gene was significantly lower in osteosarcoma cell lines (U2, MG63,143B) than in hFOB1. 19 (Figure 6G), which suggested that the expression of the gene was significantly lower in osteosarcoma cell lines than in normal control cells. This suggested that the gene may possibly not impact the osteosarcoma but rather the TME immune cells.

ATG16L1 Inhibited the Expression of GZMB and Promoted the Expression of ITGA1 in CD8+ T Cell

The scRNA-seq transcriptome datasets of 7 osteosarcoma cases were collected from the NCBI to investigate the effect of the ATG16L1 gene on CD8⁺ T cells. After normalization, PCA was used to achieve dimensionality reduction. Uniform Manifold Approximation and Projection (UMAP) analysis was performed to display the dimensional-reduced data in a 2D space, and then the cells were overlapped into 7 cell clusters (Figure 8A–B). With the existing cell types as reference, the cell clusters were annotated by marker genes (Figure 8C). Therefore, the cell types, marker genes, and proportions of these 7 cell clusters were as follows: chondroblastic cells (Runx2, Sox9, IBSP, Acan),^16,27 macrophage (CD163, LYZ),²⁸ osteoclasts (ACP5, CTSK, ATP6V0D2),^16,29,30 osteoblastic cells (ALPL),³¹ endothelial cells (EGFL7, PLVAP),^32,33 T cells (CD3D, CD3E, CD3G),³⁴ and cancer-associated fibroblasts (PDGFRB, ACTA2).^32,35 The number and proportion of these cells are shown in Figure 8B, in which the number and proportion of T cells were 3288 and 4.73%, respectively. Similarly, the data of T cells were extracted to perform UMAP analysis to obtain 3 cell clusters (Figure 8D), which were annotated by marker genes. The cell types and marker genes included CD8⁺ T cells (CD3E, CD3G, CD8A, CD8B),^36,37 Tregs (CTLA4, TIGIT, CD4)³⁸ and T cells (CD3D, LAG3)³⁹ (Figure 8E). Then, CD8⁺ T cells were divided into the ATG16L1⁺ and ATG16L1^– groups. The cytotoxic effector molecules of CD8⁺ T cells were screened, and the results revealed that the expression levels of granzyme B (GZMB) and integrin subunit alpha 1 (ITGA1) were significantly different between the groups (Figure 8F–G). In the ATG16L1⁺CD8⁺ T cell group, the expression of GZMB was lower, whereas the expression of ITGA1 was higher. GSVA revealed that the CD8 TCR Pathway, CD8 TCR Downstream Pathway, Effector VS exhausted CD8 T cell UP, PD1 high vs PD1 low CD8 T cell DN, Cancer immunotherapy by PD1 blockade, Cancer immunotherapy by CTLA4 blockade were upregulated, whereas the Naive vs EFF CD8 T cell UP, Tumor invasiveness 2D, and PD1 high vs PD1 low CD8 T cell UP were downregulated in the ATG16L1⁺CD8⁺ T cell group (Figure 8H).

Figure 8 Heterogeneity of CD8+ T cell populations in osteosarcoma. (A) UMAP plot showing the cells in the TME, with the different colors representing the different cell types. (B) Composition of the cells in the TME. (C) Marker genes of cells in the TME. (D) UMAP plot showing the CD8+ T cells, T cells and Tregs, with the different colors representing the different cell types. (E) Marker genes of the CD8+ T cells, T cells, and Tregs. (F) (G) The expressions of GZMB and ITGA1 between the ATG16L1+CD8+ T and ATG16L1−CD8+ T cells groups. (H) GSVA showing the function of the ATG16L1+CD8+ T cells group and ATG16L1−CD8+ T cells group. *P < 0.05. P < 0.05 was considered indicative of significant difference. Abbreviations: UMAP, uniform manifold approximation and projection; GSVA, gene set variation analysis; TME, the tumor microenvironment.

Discussion

The primary goal of our research was to provide a signature and therapeutic target to help improve the prognostic accuracy of osteosarcoma patients. Currently, metastasis and relapse of osteosarcoma are unavoidable risks. Lung metastasis and recurrence are the main challenge in the treatment of osteosarcoma. Most patients with osteosarcoma have locally-advanced disease or lung metastasis at the time of diagnosis.⁴⁰ Advances in understanding of the genetic determinants of osteosarcoma can facilitate risk stratification of patients, which potentially can outperform the current pathological staging system.⁴¹ Osteosarcoma cells are formed in the TME, which contains a complex mixture of other cells (MSCs, fibroblasts, osteoblasts, and myeloid immune cells such as CD8⁺ T cells).⁴² This special TME environment is suitable for the development and metastasis of osteosarcoma.

At present, the common methods used for analysis of the TME include bulk and single-cell RNA sequencing. In bulk RNA sequencing, we analyzed the GSE152048 dataset, which contained microarrays of pre-chemotherapy biopsies of osteosarcoma patients. Bioinformatics analysis was used to identify the hub genes. Fifteen hub genes were identified by limma analysis, and combined with PPI network, only nine genes left (Figures 3 and 4); the ATG16L1 gene was significantly upregulated in the high-risk group. Survival analysis (Figure 5A–E) showed that only the ATG16L1 gene was highly expressed in the high-risk group. The immune cell fraction was estimated. CD8⁺ T cells, helper T cells, and monocytes were lower in the ATG16L1⁺ group (Figure 7). The results of qRT-PCR showed significantly lower expressions of ATG16L1 in osteosarcoma cells compared with that in the control group (Figure 5F). These results suggest that ATG16L1 gene may influence osteosarcoma by acting on immune cells in the TME.

In the past few years, the ability to perform single-cell expression profiles to thousands of cells in a single experiment has been demonstrated.⁴³ ScRNA-seq offers advantages in characterizing tumor heterogeneity and functional analysis of immune cells.^44,45 In the present paper, we discuss the role of ATG16L1 gene at the single-cell level. Analysis of the scRNA-seq transcriptome datasets showed that the expression levels of GZMB in the ATG16L1⁺CD8+ T cell group were significantly lower than those in the ATG16L1⁻CD8⁺ T cell group, and the expression levels of ITGA1 were significantly higher in the ATG16L1⁺CD8⁺ T cell group than in the ATG16L1⁻CD8⁺ T cell group (Figure 8A–H). These findings indicate that ATG16L1 can inhibit the expression of GZMB and promote the expression of ITGA1.

It has been reported that GZMB is usually expressed in immune cells, such as CD8⁺ T cells and NK cells. GZMB is one of the most potent cytotoxic molecules, which mainly has an anti-tumor effect in the TME.^46,47 CD8⁺ T cells with high expression of GZMB were shown to be associated with prolonged survival of osteosarcoma patients.⁴⁸ ITGA1 was shown to be is a malignant biomarker in pancreatic cancer.⁴⁹ Some studies have shown that ITGA1 acts as an oncogene in the context of gastric cancer, colorectal cancer and lung cancer.^50–52 Therefore, ATG16L1 can inhibit the expression of GZMB and promote the expression of ITGA1, subsequently reducing the lethality of CD8⁺ T cells to tumor cells.

ATG16L1 is a component of a large protein complex necessary for autophagy. Autophagy is widespread phenomenon in cells. Loss of autophagy can enhance rejection of tumors mediated by CD8⁺ T cells. CD8⁺ T cells play an essential role in anti-tumor immunity. Suppression of autophagy induces glycolytic phenotype in T cells leading to a decrease in S-adenosyl methionine. As a result, autophagy-deficient T cells transcribe and recombine immune response genes to an effector memory state. ATG16L1 plays a crucial role at various stages of autophagy.⁵³ Mice deficient in the essential autophagy gene ATG16L1 were shown to display a dramatic impairment in the growth of syngeneic tumors. Autophagy represents a new direction in the treatment of osteosarcoma. Regulators or inhibitors of autophagy may provide a potential target for cancer treatment.²⁶ Furthermore, ATG16L1 can prevent the activity of CD8⁺T cells through mediating by dendritic cells.⁵⁴ Thus, the ATG16L1 gene can play an important role in the development and metastasis of osteosarcoma by mediating CD8⁺T cells.

Immunotherapy has become one of the most effective treatments for many tumors. Currently, there is a lack of specific immunosuppressive drugs for osteosarcoma. Immunotherapy (such as immune checkpoint inhibitor therapy) is related to T cell response. CD8⁺ T cells mediated by ATG16L1 may affect immune function and promote the survival of osteosarcoma.

Some limitations of our study should be acknowledged. Firstly, the bulk RNA sequencing and scRNA-seq datasets have relatively small samples; Secondly current research methods can only analyze 22 types of immune cells, but cannot determine the characteristics of uncommon immune cell populations. Third, scRNA-seq reduces the sensitivity of each cell. It may not be able to recognize subtle transcriptional differences between cells, which may affect gene expression.

Conclusion

Our results demonstrate that ATG16L1 can affect the prognosis of osteosarcoma by affecting CD8+ T cells, especially in patients with metastases. ATG16L1 is a potential prognostic biomarker and may be a novel therapeutic target for osteosarcoma.