Identication of a Potential Autophagy-Associated Therapeutic Target in Osteosarcoma using Bioinformatic Analyses

Background: Osteosarcoma, primarily resulting from mesenchymal cell differentiation, is the most common primary malignancy of the bone in children and adolescents. Metastases to other sites, such as the lung, often occur in the early stages and progress rapidly. Autophagy functions as a tumor-suppressing mechanism in the early stages of oncogenesis, however, in later stages, autophagy may promote tumor growth, spread, and increase treatment resistance. Methods: In the present study, we aimed to screen new key autophagy-related biomarkers that may serve as therapeutic targets for osteosarcoma. The expression prole of the GSE42352 dataset, including 103 OS cell lines and 15 MSC lines from Gene Expression Omnibus (GEO), was analyzed to identify the differentially expressed genes. Results: WGCNA was used to construct the gene co-expression network in osteosarcoma, identify co-expression modules for protein interactions (PPI), and screen the key genes. A total of 757 differentially expressed genes were identied by WGCNA analysis including one key autophagy-related module. Conclusions: Further, we performed the PPI analysis for module genes and identied a key gene. We also theoretically and functionally validated its expression using the validation dataset GSE14359, and in vitro experiments, respectively.


Introduction
Osteosarcoma (OS) is one of the most common malignancies of the bone in children and adolescents under the age of 20, accounting for approximately 5% of childhood tumors; it is associated with a high lethality rate (1)(2)(3). Before 1970, amputation was the standard treatment for osteosarcoma, but the veyear overall survival rate was still low, with most patients dying within one year of diagnosis (4). The current 5-year survival rate of patients ranges from 20-60% or 70%(5, 6).
Autophagy was discovered and its mechanism was proposed by Ashford and Porter in 1962(7).
Autophagy is a form of cell death where the cells engulf their cytoplasmic proteins or organelles and encapsulate them into vesicles where the contents are eventually degraded(8). It is a dynamic process that promotes organelle and protein turnover. It also generates metabolic precursor molecules through lysosome-dependent degradation of macromolecules, organelles, and other cellular components (9). It is regulated by a variety of autophagy genes, also known as type II programmed cell death. Some studies show that inhibition or activation of autophagy in tumor cells can promote apoptosis and signi cantly enhance chemotherapeutic e cacy (10). Autophagy promotes tumor growth by mediating levels of cytokines, mTOR, Beclin-1, miRNAs, etc. for the removal of damaged proteins and organelles. For example, during chemotherapy, autophagy is inhibited by miR-22 in osteosarcoma cells (11); Beclin 1 gene plasmid transfection into the MG63 human osteosarcoma cell line showed that with overexpression of Beclin 1, LC3B expression and autophagic activity were enhanced, indicating that Beclin 1 promotes the autophagic process in this tumor cell line (12); the classical autophagic signaling pathway functions through mTOR, downstream of PI3K/AKT, and is important for controlling translation, cell cycle progression, and negative regulation of autophagy (13). The PI3K/AKT/mTOR pathway has received great attention from researchers in recent years (14,15)as a major upstream regulator of the autophagic pathway. mTOR-containing protein complexes and other proteins in this pathway are known to be direct or indirect targets of many miRNAs. Thus, a deeper understanding of their underlying mechanisms of action would enhance the treatment strategy for osteosarcoma.
Autophagy can reduce the likelihood of apoptosis by maintaining a range of normal mitochondrial functions such as antioxidant role, anti-DNA damage, and adenosine triphosphate production which promote tumor metastasis(16). It similarly inhibits tumor metastasis by suppressing in ammation and inhibiting macrophage in ltration by parallelly reducing metabolic stress and preventing necrosis-induced death along with apoptosis (17). Thus, autophagy is a complex process involving multiple molecules. It is important to uncover and explore the key factors underlying autophagy for tumor treatment and patient prognosis. In this study, we identi ed autophagy-associated genes using a bioinformatic approach for osteosarcoma, which can serve as potential biomarkers for the diagnosis and treatment of osteosarcoma.

Differential gene expression
The"affy" and the "impute" R packages in R/Bioconductor software were used for GEO data processing, differential gene expression was analyzed with the Limma package (3.40.2) in R. Adjusted P-values were evaluated in GEO to correct for false-positive results. "adjusted P-value < 0.05 and |logFC| ≥ 1" were de ned as the cut-off threshold for differential mRNA expression; heat and volcano plots were plotted using the 'ggplot2' package in R.

Autophagy-associated consensus gene clustering
Consensus analysis was performed using the ConsensusClusterPlus package in R (v1.54.0). The maximum number of clusters up to 6 and 100 repetitions were used to extract 80% of the total samples using clusterAlg = "hc", innerLinkage='ward.D2'. The heat map for clustering was analyzed by pheatmap package in R (v1.0.12). The genes with expression variance above 0.1 were included in the heat map; if the number of input target genes was above 1000, the top 25% genes were extracted for display after sorting in descending order of variance. All of these analyses were performed using R v4.0.3.

Weighted Gene Co-expression Network Analysis (WGCNA)
For the GSE42352 dataset, the WGCNA package in R was used to calculate the Pearson correlation coe cient among the differentially expressed genes. An appropriate soft threshold β was selected to construct the network consistent with the criteria of a scale-free network. A one-step method was used to construct the gene network, transform the adjacency matrix into a topological overlap matrix (TOM), and generate a hierarchical clustering tree of genes. Gene and module signi cances were calculated. These indicated the signi cance of genes with clinical information. The module identity (MM) of each gene was calculated and it indicated the importance of the individual gene in the module. Cut-offs were set at |MM|> 0. 8 and |GS|> 0. 2 for gene identi cation.

Functional annotation and pathway enrichment analysis
Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed for all genes in the red module using Database for Annotation, Visualization, and Integrated Discovery (DAVID) (v6.8).

Protein-protein interaction network (PPI) and hub genes
The STRING database (https://string-db.org/) was used to analyze known and predicted protein-protein interactions. All genes in the red module were analyzed and PPI networks were constructed using STRING. Further, MCODE and Cytohubba in Cytoscape (v3.8.2) software were executed to screen the hub genes in PPI networks.

Single gene enrichment analysis (GSEA)
The samples were divided into two groups based on the median gene expression as high-and lowexpression groups. Gene set enrichment analysis(GSEA) was performed to investigate the functions correlated with different risk groups by GSEA 4.1.0, and the software was downloaded from the website for GSEA (http://www.gseamsigdb.org/gsea/downloads). The association of gene expression with tumors was examined using GSEA. The threshold cut-offs were set at |NES|≥1, FDR<0.25 and P<0.05.

Cell proliferation assay
Cell proliferation was assessed using the Enhanced CCK-8 kit (Beyotime Biotechnology, C0042). Cells were inoculated at a density of 2×10 3 cells/well into 96-well plates for 24h, 48h, 72h, and 96h. 10µL of enhanced CCK-8 solution was added to each well and incubated for 4h at room temperature(20℃-25℃). The absorbance at 450nm was measured using an enzyme marker.

Trans-well invasion assay
Cells were inoculated at a density of 1×10 5 cells/well into the upper chamber of trans-well pre-wrapped with Matrigel. The lower chamber contained the medium supplemented with 10% FBS The cells were incubated for 24h at room temperature; xed using paraformaldehyde for 15min; stained with crystal violet for 10min, cells were visualized and counted from six randomly selected elds under a microscope (200 × magni cation).

Cell scratching assay
Cells were inoculated in 24-well plates at a density of 1×10 5 cells/well. After growing to con uence, the cells were scratched using a sterile small gun tip, washed, and incubated for 24 h. The wound healing was observed under the microscope.

Statistical analysis
All experiments in this study were repeated at least three times. Statistical analyses were performed using GraphPad Prism 8.0 software. Data were expressed as mean ± standard deviation (x ± s), and a t-test was used for the comparison of two independent data groups. The differences were considered statistically signi cant at P<0.05.

DEGs in GSE42352
A total of 103 cancer-and 15-normal cell samples were obtained from the GSE42352 microarray. Heat maps were plotted for the expression of genes in each sample ( Figure 1A); a total of 439 up-regulated and 318 down-regulated genes were obtained as presented using the volcano map ( Figure 1B). The differentially expressed genes were taken for further functional enrichment analysis, and the top 20 annotated functions and enriched pathways are shown in Figure 1C.

Consensus clustering based on autophagy-associated genes
To analyze the relationship between the expression of 222 autophagy-related genes and osteosarcoma subtypes, we performed a consensus clustering analysis for 118 samples from GSE42352. On increasing the clustering variable (k) from 2 to 6, we found that at k = 4, the intra-group correlation was the highest and the inter-group correlation was the lowest. This indicated that the 118 samples were well classi ed into the four clusters (Figure 2A-C). The heat map in Figure 2D shows the gene expression pro les, clinical characteristics, and the ConsensusClusterPlus consistency clustering at k = 4.

Identi cation of gene co-expression modules
To attain the prerequisites of scale-free network distribution, the values of the adjacency matrix weight parameter β needed to be examined. By setting the range of network construction parameters selection, the scale-free distribution topology matrix was calculated. In this study, a soft threshold value of 3 ( Figure  3A) was selected to ensure that the scale-free network satis ed the hierarchical clustering, and a total of 36 modules were identi ed ( Figure 3B). The clustering correlation between modules is shown in Figure  3C. The red module was con rmed to be signi cant by analyzing the correlation of each module with autophagic subtypes and types ( Figure 3D-E). Subsequently, KEGG and GO functional enrichment analyses were performed for 1991 genes in the red module, and the signi cant pathways enriched for the top 20 genes are shown ( Figure 3F-G).
(A: soft threshold screening results of constructing scale-free network; B: identi cation results of TOMbased gene system clustering tree, different colors in the gure represent different multi-gene modules. c: clustering relationship between WGCNA modules and modules; D: heat map of correlation between related traits and module features; E: correlation between red modules and Gene signi cance; F: KEGG enrichment (analysis of red module genes; G: GO enrichment analysis of red module genes;)

Construction of PPI network
The DEGs in the red module were analyzed online using the STRING database to obtain the protein interaction network. We visualized the gene information and network construction as shown in Figure 4A. The MCODE plug-in of Cytoscape was used to calculate the characteristics of each node in the network graph. We selected the largest sub-network and visualized it as shown in Figure 4B. MCODE1 functional enrichment is shown in Table 1. Core genes were analyzed with the CytoHubba plug-in; the top 10 genes were identi ed as core genes using the Closeness algorithm ( Figure 4C). These included TLR4, TLR8, ITGAX, ICAM1, CD86, CXCR4, ITGAM, TLR2, TNF, and PTPRC. The overlap between the two analyses showed that CXCR4 was the common core gene ( Figure 4D).
(A: PPI plots of genes in the red module; B: MCODE plug-in ltering out the highest-scoring sub-networks; C: Cytohubba plug-in screening the top ten hubba gene networks of the proximity algorithm; D: intersection genes of MCODE1 genes and proximity genes;) 3.5 CXCR4 expression in different datasets CXCR4 expression was analyzed in different datasets. In the GSE42352 dataset, the expression was signi cantly high in osteosarcoma cells as compared to the normal cells ( Figure 5A). CXCR4 expression differences were statistically signi cant in all four consensus clustering subtypes of the GSE42352 dataset ( Figure 5B); CXCR4 expression was upregulated in the GSE14359 dataset, in the osteosarcoma samples as compared to normal samples ( Figure 5C).
(A: CXCR4 expression in the GSE42352 dataset; B: CXCR4 expression in different consensus clustering subgroups; C: CXCR4 expression in the GSE14359 dataset;)

Functional annotation of single gene CXCR4
The ve KEGG and HALLMARK pathways that were most signi cantly associated with high CXCR4 expression are shown in Figure 6. Among these, CXCR4 was mainly enriched in arachidonic acid metabolism, natural killer cell-mediated cytotoxicity, systemic lupus erythematosus, toll-like receptor signaling pathway, and B-cell receptor signaling pathway ( Figure 6A). Additionally, CXCR4 high expression was positively associated with complement pathway, apoptosis, coagulation, IL-6/JAK/STAT3 signaling system, and heme metabolism ( Figure 6B).

In vitro CXCR4 assay
In the present study, we evaluated CXCR4 expression in cell lines by quantitative real-time PCR,the results showed that CXCR4 siRNA led to a marked reduction of CXCR4 mRNA expression level( Figure 7A). CCK-8 assay was used to detect the effect of CXCR4 on cell proliferation. The results indicated that the proliferation of U2OS cells transfected with si-CXCR4 was signi cantly inhibited ( Figure 7B). Trans-well invasion assay and cell scratch assay were used to analyze the effect of CXCR4 on cell invasion and migration. The results showed that the invasion ability of U2OS cells transfected with si-CXCR4 was signi cantly reduced and the scratch spacing was signi cantly increased as compared to the untransfected U2OS cells and si-NC group. Thus, the invasion and migration of the cells were diminished ( Figure 7C-7D). We used western blotting to examine the expression of autophagic pathway-related proteins in U2OS cells. si-CXCR4 transfection signi cantly decreased the levels of PI3K, AKT, and mTOR phosphorylation and elevated the levels of Caspase-3 protein expression in U2OS cells ( Figure 7E). These results suggested that CXCR4 silencing induced autophagic cell death through the inhibition of the PI3K/AKT/mTOR signaling pathway but did not inhibit the growth in osteosarcoma.

Conclusion
In recent years, the number of studies based on public databases such as GEO and TARGET has increased. These are useful to mine potential biomarkers for osteosarcoma. For example, Dai et al. (20)screened candidate genes to predict the chemotherapy resistance response of osteosarcoma by miRNA-mRNA interaction network. Jiang et al. (21) identi ed autophagy and immune-related gene markers TRIM68, PIKFYVE, and DYNLL2, which could predict the prognosis of osteosarcoma. In the present study, using the GSE42352 dataset and by differential WGCNA analysis, we screened autophagyrelated modular genes and construct a PPI network. A key gene, CXCR4, was identi ed by MCODE and Closeness algorithms. 222 autophagy-related genes were obtained from the online database. We performed concordance analysis for samples in GSE42352, and four subgroups were obtained; CXCR4 was identi ed in all four subgroups. Using the validation dataset GSE14359, we showed that CXCR4 was indeed highly expressed in osteosarcoma as compared to normal samples. Niu et al. (22) have also identi ed CXCR4 as a possible therapeutic agent to inhibit the progression of osteosarcoma by bioinformatic analysis.
Chemokine receptor-4 (CXCR4) is a speci c receptor for chemokine stromal cell-derived factor-1 (CXCL12) (23). There is increasing evidence that CXCR4 plays a crucial role in osteosarcoma progression and metastasis (24,25). In addition, CXCR4 is associated with poor survival in patients with osteosarcoma and is considered to be an important clinical prognostic indicator(26). In the present study, the module selected for WGCNA analysis showed the highest signi cant correlation with autophagy. Autophagy, a highly evolutionarily conserved cellular homeostatic dynamic recycling system, is used for cytoplasmic recycling, degradation, and reuse (27). Autophagy selectively degrades long-lived proteins, removes dysfunctional organelles under basal conditions, and promotes cell survival(28). It has been recognized as a cytoprotective process that promotes chemoresistance in osteosarcoma. Several studies have focused on autophagy inhibition in osteosarcoma chemo-sensitization (29). Coly et al. (30) found that, in HEK293 and U87 cells, activation of CXCR4 leads to a decrease in the number of autophagosomes, suggesting that CXCR4 exerts its anti-autophagic effect by regulating calcium-activating enzymes and preventing pre-autophagosomal vesicle formation. In the present study, our in vitro assays also showed that CXCR4 could affect the progression of osteosarcoma by regulating the PI3K-Akt-mTOR signaling pathway. The main regulator of autophagy is mTORC1; mTORC2 also contributes to the autophagic control through AKT pathway regulation (31). Under normal conditions, mTORC1 directly inactivates the ULK1/2 protein complexes. These complexes are critical and essential regulators of cellular pathways that control the initiation of gene translation and ribosome biogenesis and exhibit important monitoring roles for cellular metabolism, lipolysis, and autophagy (32,33). In general, there is a negative interaction between autophagy and apoptosis, as autophagy prevents apoptosis induction and apoptosis-associated caspase activation inhibits the autophagic process. Some studies also show that the negative correlation between CXCR4 and autophagy is dependent on PI3K/AKT/mTOR signaling pathway (34).
In the present study, by single gene enrichment analysis, we showed that CXCR4 is mainly enriched in arachidonic acid metabolism, natural killer cell-mediated cytotoxicity, and IL-6/JAK/STAT3 signaling pathway in osteosarcoma. Previous studies (35) have shown that CXCR4 is itself involved in NK cellmediated cytotoxicity in chronic granulocytic leukemia. IL-6/JAK/STAT3 is a stress-related in ammatory signaling pathway, with rapid response. Dysregulated cytokine signaling contributes to cancer development (36). In addition, the STAT family of transcription factors can be regulated by other kinases, such as mTOR, MAPK which phosphorylate Ser727 of STAT3 for sustained activation (37). STAT3 is potentially pro-oncogenic (a proto-oncogene) and is consistently expressed across several cancers. In the pathogenesis of atherosclerosis, CXCR4 is involved in chronic in ammation of the arterial wall and is characterized by a chemokine-mediated inward ow of leukocytes(38). Chronic in ammation and local in ltration of CXCR4-expressing immune cells can promote esophageal carcinogenesis (39). Activating downstream cascades involving JAK/STAT, PI3K/Akt, MAPK, JNK pathway considered important targets for development of new therapeutic strategies, and the activation can be triggered by the binding between CXCR4 and its ligand CXCL12. These critical points control stemness, chemotaxis and cell survival, proliferation, migration (40). CXCL12 (41,42) is a homeostatic chemokine produced by mesenchymal stem cells,which binds to chemokine receptors to regulate hematopoietic stem cell (HSC) transport and secondary lymphoid tissue structure. CXCR4 overexpression has been found in many tumors (43). There are evidences that CXCR4/CXCL12 axis activates multiple downstream pathways and plays an important role in tumor progression (41,44).Besides, previous study have found that the overwhelming majority of metastases from osteosarcoma are to the lung, a site that expresses high levels of the ligand CXCL12 (45).
Taken together, our results systematically demonstrated high CXCR4 expression in osteosarcoma. We speculated that it could induce autophagic cell death by regulating the PI3K/AKT/mTOR signaling pathway. Indeed, functional enrichment analysis showed that CXCR4 in in ammatory, immune-related contexts could provide new strategies for gene-targeted therapy for osteosarcoma.

Declarations
Ethics approval and consent to participate: none Consent for publication: none Availability of data and materials: none Competing interests: none Funding: none