Integral Analyses of Competing Endogenous RNA Mechanisms and DNA Methylation Reveal Regulatory Mechanisms in Osteosarcoma

Background: Osteosarcoma (OS) is the most common primary malignant bone tumour in children and adolescents, with rapid growth, frequent metastasis, and a poor prognosis, but its pathogenesis has not been fully elucidated. Exploring the pathogenesis of OS is of great significance for improving diagnoses and finding new therapeutic targets. Methods: Differentially expressed circRNAs (DECs), miRNAs (DEMs), methylated DNA sites (DMSs), and mRNAs (DEGs) were identified between OS and control cell lines. GSEA of DEGs and functional enrichment analysis of methylated DEGs were carried out to further identify potential biological processes. Online tools were used to predict the miRNA binding sites of DECs and the mRNA binding sites of DEMs, and then construct a circRNA-miRNA-mRNA network. Next, an analysis of the interaction between methylated DEGs was performed with a protein-protein interaction (PPI) network, and hub gene identification and survival analysis were carried out. The expression pattern of circRNA-miRNA-mRNA was validated by real-time PCR. Results: GSEA and functional enrichment analysis indicated that DEGs and methylated DEGs are involved in important biological processes in cancer. Hsa_circ_0001753/has_miR_760/CD74 network was constructed and validated in cell lines. Low expression levels of CD74 are associated with poor overall survival times and show good diagnostic ability. Conclusion: Methylated DEGs may be involved in the development of OS, and the hsa_circ_0001753/has_miR_760/CD74 network may serve as a target for the early diagnosis of and targeted therapy for OS.


INTRODUCTION
Osteosarcoma (OS) is the most common malignant bone tumour and occurs primarily in the metaphysis of long bones in children and adolescents (Berner et al., 2015). The annual incidence of OS ranges from 1 to 4 in 1 million and is slightly higher in men than in women (Mirabello et al., 2009). Although the overall incidence of OS is low, its prognosis is poor. With the combined application of limb salvage surgery and neoadjuvant chemotherapy, the 5year survival rate for patients without metastases is approximately 65-70% (Siegel et al., 20182018). Approximately 15-20% of patients have evidence of metastasis at the time of diagnosis, mainly to the lung, and the prognosis for patients with metastatic disease remains poor (Meazza and Scanagatta, 2016). Many OS patients are insensitive to chemotherapy drugs or have developed resistance, and the recurrence and widespread metastasis of tumours also complicate the clinical treatment of osteosarcoma. To date, although molecular biological research related to osteosarcoma has provided a theoretical basis for its pathogenesis, the molecular mechanisms of the malignant biological behaviour of osteosarcoma is still unclear.
With the development of a new generation of sequencing technology, noncoding RNA (ncRNA) has gradually become a hot topic of research regarding tumour pathogenesis. circRNAs are ncRNAs that play a role in biological regulation, primarily through genetic regulation (Li et al., 2018). circRNA is a singlestranded endogenous noncoding RNA produced by reverse splicing; it has a covalent closed-loop structure that cannot be easily degraded by RNA nucleic acid exonuclease (Chen and Yang, 2015). circRNAs have abundant miRNA binding sites and can absorb miRNAs as sponges (Tay et al., 2014). circRNAs can compete with miRNAs to affect the stability of target RNAs or their translation, thus regulating gene expression at the posttranscriptional level (Zhong et al., 2018). At the same time, circRNAs are widely involved in the processes of cell proliferation, apoptosis, invasion and migration, playing an important regulatory role in human diseases (Wang et al., 2020a;Yu et al., 2020). Based on competing endogenous RNA (ceRNA) mechanisms, a variety of abnormally expressed circRNAs play an important regulatory role in cancer, including breast cancer (Sang et al., 2019;Yang et al., 2019), lung cancer (Cheng et al., 2019;Chen et al., 2020), liver cancer He et al., 2020) and gastric cancer Luo et al., 2020). As a ceRNA, circ-0001785 mediates miR-1200 to regulate the pathogenesis of OS by acting as a molecular sponge, thereby upregulating HOXB2 and mediating the malignant cell behaviour of OS cells . Another study revealed that the high expression of circ-03955 promotes EMT in osteosarcoma by acting as a miR-3662 sponge-mediated MTDH expression (Wang et al., 2020b). However, many biological functions of circRNAs in osteosarcoma are still unknown.
Many studies in recent years have shown that epigenetic changes play an important role in regulating the occurrence and development of tumours. DNA methylation is one of the most common epigenetic events involving the predominantly reversible addition of methyl groups to cytosines without altering the genomic DNA sequence in the context of CpG dinucleotides (Ding et al., 2019). Studies have shown that the promoter regions of tumour suppressor genes are hypermethylated and inhibited, while oncogenes are hypomethylated and abnormally active (Lin and Wang, 2014). Abnormal DNA methylation has been proven to regulate tumorigenesis and the progression of various cancers (Skvortsova et al., 2019;Wang et al., 2020c). Many studies have reported findings related to DNA methylation in OS. According to reports, overexpression of SENP3 can increase DNA methylation of E-Cad and then increase the proliferation, migration and invasion of osteosarcoma cells . The high expression level of WNT6 in primary osteosarcoma is mainly due to the low DNA methylation level of WNT6, while the DNA methylation of WNT6 shows a negative correlation with the prognosis of children with osteosarcoma (Li et al., 2017). Although the DNA methylation pattern in OS has attracted attention as an important biomarker and therapeutic target, the mechanism by which ceRNA and DNA methylation regulate OS is still not understood.
To better understand the molecular mechanisms of OS, we used gene microarray and bioinformatics analysis to study the potential mechanism of ceRNA and DNA methylation in OS. The flow chart summarizing this work is shown in Figure 1. In this study, OS-related datasets were selected from the gene microarray of the GEO database, and differentially expressed circRNAs (DECs), miRNAs (DEMs), methylated DNA sites (DMSs) and mRNAs (DEGs) were identified between OS and control cell lines. Functional enrichment analysis of DEGs and GSEA of methylated DEGs were carried out to further identify potential biological processes. On the basis of the ceRNA mechanism, online tools were used to predict the miRNA binding sites of DECs and the mRNA binding sites of DEMs and then construct a circRNA-miRNA-mRNA network. Next, a protein-protein interaction (PPI) network was established to analyse interactions between methylated DEGs, and then hub gene identification and survival analysis were carried out. Finally, the ceRNA mechanism of the circRNA-miRNA-mRNA network was verified in cell lines and tissues. The aims of our research were to identify and verify the ceRNA mechanism and methylated DEGs in OS. Our findings provide new insights into the potential carcinogenic effects of abnormal DNA methylation and identify potential biomarkers and therapeutic targets of OS. Illumina human-6 v2.0 expression BeadChip]) were downloaded from the GEO database. The GSE96964 dataset includes 7 human osteosarcoma cell lines (OS group) and 1 human osteoblast (control group); GSE70367 includes 5 human osteosarcoma cell lines (OS group) and 1 human mesenchymal stem cell line (control group); GSE36002 includes 19 human osteosarcoma cell lines (OS group) and 2 human osteoblasts (control group); GSE42351 includes 19 human osteosarcoma cell lines (OS group) and GSE33382 includes 3 human osteoblasts (control group).

Intragroup data Repeatability Test
Principal component analysis (PCA) is a general method of sample clustering analysis, which is often used for sample clustering based on various variable information, such as gene expression, resequencing, and diversity analysis. The intragroup data repeatability verification of the dataset was tested by PCA. Pearson's correlation test was used to verify the repeatability of the data in each group. The statistical analysis and visual heatmap were drawn via the R programming language and depict the correlations among intragroup data.

Microarray Analysis
The limma package (Ritchie et al., 2015) in R was used to identify differentially expressed circRNAs (DECs), miRNAs (DEMs), methylated DNA sites (DMSs) and mRNAs (DEGs) between OS and control cell lines. |Log2FC| ≥ 0.263, P value ＜ 0.05 and q value ＜ 1 were considered statistically significant differences. For the selected DECs, DEMs, DMSs and DEGs, we generated volcano plots and heatmaps using R software. Frontiers in Cell and Developmental Biology | www.frontiersin.org December 2021 | Volume 9 | Article 763347 3

Gene Set Enrichment Analysis
According to the differentially expressed genes, GSEA (Subramanian et al., 2005) was used for gene collection enrichment analysis. Genome enrichment analysis is a computational method used to determine whether a given gene set differs significantly between different groups. The genes in these clusters are related in some way, such as having the same biological function, being located close to each other on the chromosome, or being self-defined according to certain criteria. Therefore, gene set enrichment analysis can compensate for the deficiencies of single-gene analysis.

GO and KEGG Analysis of Methylated DEGs
Gene Ontology (GO) analysis is widely used for annotating genes and determining biological characteristics, including three key biological aspects: cellular components (CC), molecular function (MF), and biological process (BP) (Ashburner et al., 2000). The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database is one of the most commonly applied databases for systematically analysing gene functions (Kanehisa and Goto, 2000). The Database for Annotation, Visualization and Integrated Discovery (DAVID, https://david.ncifcrf.gov/) is a gene functional classification tool designed to provide authors with a comprehensive set of functional annotation tools to understand the biological significance behind a large number of genes, and which can be used to perform GO and KEGG analysis. A P value < 0.05 was considered statistically significant.

PPI Network Construction and Hub Gene Identification
For genes in the ceRNA network, the Search Tool for the Retrieval of Interacting Genes (STRING) (Mering et al., 2003) (https://string-db. org/) online tool was used to analyse protein-protein interactions (PPIs). The required confidence (combined score) > 0.15 was selected as the threshold for protein-protein interactions. Based on the obtained PPI relationship pairs, Cytoscape was used to analyse the topology of the PPI relationship network. From the obtained biological networks, we observed that most biological networks conformed to the attributes of a scale-free network. Therefore, the primary node involved in the protein interaction relationship in the PPI network, namely, the hub gene (He and Zhang, 2006), could be obtained by using connectivity degree analysis in network statistics. In this study, the nodes of the network were analysed, and the hub genes were identified with degree > 2.

Survival Analysis of the Hub Genes
The study attempted to validate the hub genes from the TCGA datasets of OS. The expression level of each hub gene was extracted for further analysis from all included data. Survival analysis and ROC curve analysis were performed by using the Survival Roc package according to the OS expression values of key genes in TCGA datasets. Survival analysis was carried out through Kaplan-Meier analysis, and the log-rank test was performed to evaluate the statistical significance of differences. The significance cut-off was fixed at a p value < 0.05.

RNA Isolation and RT-PCR
Total RNA was isolated from cultured cells using TRIzol reagent (Invitrogen, Carlsbad, CA, United States) according to the manufacturer's instructions. The extracted RNA was quantified by spectrophotometry using NanoDrop equipment (Thermo Fisher Scientific, Waltham, MA, United States). Reverse transcription of mRNA and miRNA was conducted using random primers and tail primers in the TaKaRa system (Dalian, China), respectively. Real-time PCR was performed using TB Green PCR Master Mix (Takara) with an ABI StepOne Real-Time PCR System (Applied Biosystems 7500, Foster City, CA, United States). The expression levels of circRNA and mRNA were normalized according to GAPDH expression, and miRNA was normalized according to FIGURE 3 | Heatmaps and volcano plots of DECs, DEMs, DMSs and DEGs in OS. In the hierarchical cluster heatmap, red and blue indicate high and low expression, respectively, with columns representing samples and rows representing, respective RNAs or methylated sites. Volcano plots were used to distinguish the differentially expressed RNAs or methylated sites. The vertical lines correspond to a difference of 0.263-fold, and the horizontal line represents a P value of 0.05. (A,B) Twenty-one upregulated and 150 downregulated circRNAs were identified from differentially expressed circRNAs; (C,D) 15 upregulated and 66 downregulated miRNAs were identified from differentially expressed miRNAs; (E,F) 3,025 methylated DNA sites were identified from differentially methylated DNA sites; (G,H) 1,184 upregulated and 1,513 downregulated mRNAs were identified from differentially expressed mRNAs.
Frontiers in Cell and Developmental Biology | www.frontiersin.org December 2021 | Volume 9 | Article 763347 5 U6 expression by using an optimized comparative Ct (2−ΔΔCt) value method. Real-time PCRs were performed in triplicate.

Tissue Microarray and Immunohistochemistry
Paraffin-embedded TMAs used in the present study were purchased from Bioaitech (Bioaitech, Xian, China). Osteosarcoma and normal bone tissue TMAs contain 70 cancer samples and 20 normal bone tissue samples. IHC was utilized to determine the relative expression of CD74 in 90 samples by means of the manufacturer's instructions. Primary antibodies against CD74 (1:500, Cell Signaling Technology, United States) were used. HRP-labeled streptavidin solution was added into the samples for 15 min after both primary and secondary antibody incubations. The immunocomplex was visualized with DAB, and the nucleus were counterstained with haematoxylin. Pictures were taken with a biochip scanner (Pannoramic MIDI, 3D HISTECH, Hungary). The expressions of CD74 were quantified and analyzed. The staining was scored according to the previously described 4-point system (score 0-3) as follows: score 3, dark staining that is easily visible and present in >50% of cells; score 2, focal areas of dark staining (<50% of cells) or moderate staining of >50% of cells; score 1, focal moderate staining in <50% of cells or pale staining in any proportion of cells not easily observable at low power; and score 0, none of the aforementioned (Deng et al., 2019). A high level of expression was defined as a score of 2-3 and low level of expression was defined as a score of 0-1, as described previously. Bisulfite Sequencing PCR Genomic DNA from osteoblasts and osteosarcoma cells was isolated using the Ezup Column Animal Genomic DNA Purification Kit (Sangon Biotech, B518251) according to the manufacturer's instructions. By treatment of genomic DNA with sodium bisulfite, all unmethylated cytosine (C) is converted to uracil (U) while methylated cytosine remains unchanged. After the sulfite treatment of genomic DNA, BSP primers were designed to amplify the target fragments, at which time uracil (U) was all converted into thymine (T). Finally, the PCR products were sequenced to determine whether methylation of CpG sites occurred. The CpG island of 2000 bp sequence before exon 1 of CD74 gene was predicted. The primers used to amplify the bisulfitetreated DNA were designed using MethPrimer software. The BSP primers were as follows: F: 5′-TGTTGTTTAGTTAGTGTAGG-3'; R: 5′-AAAAAAACCAAACACRATAACTCA-3'. The PCR products were purified and cloned, and 10 clones from each group of samples were selected for sequencing, and sequencing was done by Sangon Biotech.

Statistical Analysis
Statistical analyses were conducted by GraphPad Prism 7.0 (GraphPad Software Inc., CA, United States). Data are expressed as the mean ± standard deviation (SD) and were assessed using unpaired Student's t-tests.

Validation of the Datasets
To further validate the repeatability of intragroup data, PCA and Pearson's correlation test were used. Based on the PCA, the intragroup data repeatability for the circRNA, miRNA, methylated DNA and mRNA datasets was acceptable. In all four datasets, the distance between the intragroup samples in the control and OS groups was close, and the distance between the control and case groups was considerable (Figure 2A-D). Based on Pearson's correlation test, it was found that in the circRNA, miRNA, methylated DNA and mRNA datasets, there was a strong correlation within the control group samples, and there was also a strong correlation within the OS group samples. The heatmap ( Figure 2E-H) and scatter diagram ( Figure 2I-L) depict Pearson's correlation.

DECs, DEMs, DMSs, DMGs, and DEGs
After expression profile processing, there were a total of 4,660 circRNAs, 1,303 miRNAs, 26,780 methylated DNA sites and 24998 mRNAs in the expression profiles used for identifying differential expression. According to the cut-off criteria of | Log2FC| ≥ 0.263, p value ＜ 0.05 and q value ＜ 1, a total of 171 differentially expressed circRNAs (DECs, 21 up-and 150

GSEA
To identify the biological functional mechanism in OS, the expression matrix constructed from DEG analysis was used for GSEA based on the C2 database. GSEA revealed that the most significantly enriched gene set in OS was the REACTOME DNA REPLICATION, and other significant gene sets included REACTOME MITOTIC METAPHASE AND ANAPHASE, REACTOME REGULATION OF MITOTIC CELL CYCLE, REACTOME REGULATION OF TLR BY ENDOGENOUS LIGAND, REACTOME ARACHIDONIC ACID METABOLISM and BIOCARTA EICOSANOID PATHWAY (Figure 4). These enriched pathways revealed that DEGs were closely associated with the malignant characteristics of OS, especially apoptosis and the cell cycle.
To further understand the gene function, we performed 134 methylated DEG enrichment functions determined by the DAVID online analysis tool, and the GO analysis of DEGs was mainly classified into three functional groups: BP, CC, and MF ( Figure 5B). In the BP group, the GO term analysis indicated that the methylated DEGs were primarily enriched in functions such as negative regulation of cartilage development,  Frontiers in Cell and Developmental Biology | www.frontiersin.org December 2021 | Volume 9 | Article 763347 8 regulation of the peroxisome proliferator activated receptor signalling pathway, regulation of cartilage development, transforming growth factor beta receptor complex assembly, negative regulation of the DNA damage response and signal transduction by the p53 class mediator (Figures 5C,D). The top 15 enriched GO-BP terms of the methylated DEGs are listed in Table 1. In the CC group, the results of GO analysis were enriched in extracellular exosomes, lysosomal lumen, extracellular vesicles, integral components of the luminal side of the endoplasmic reticulum membrane, and the luminal side of membranes ( Figures 5E,F). The top 15 enriched GO-CC terms of the methylated DEGs are listed in Table 2. In the MF group, the enriched GO terms were primarily involved in extracellular matrix structural constituents conferring compression resistance, chemokine receptor binding, chemokine activity, four-way junction DNA binding and E-box binding ( Figures   5G,H). The top 15 enriched GO-MF terms of the methylated DEGs are listed in Table 3.
KEGG pathway enrichment analysis showed that the top pathways related to the methylated DEGs were the chemokine signalling pathway, viral protein interaction with cytokine and cytokine receptors, the JAK-STAT signalling pathway, viral myocarditis, and the adipocytokine signalling pathway ( Figures 5I,J). The top 15 enriched KEGG pathway terms of the methylated DEGs are listed in Table 4. ceRNA Network of circRNA-miRNA-mRNA To better understand the interaction among the DECs, DEMs, and methylated DEGs of OS, the ceRNA network of circRNA-miRNA-mRNA was constructed. First, the miRNA binding sites of DECs were predicted by the starBase database, and a total of 7 upregulated DEC-downregulated DEM pairs and 198    downregulated DEC-upregulated DEM pairs were obtained after integration with the identified DEMs through a ceRNA mechanism. Second, the mRNA binding sites of DEMs were predicted by the miRanda, miRDB, miTarBase, TargetScan and miRMap databases, and a total of 152 upregulated DEMhypermethylated DEG pairs and 284 downregulated DEM-hypomethylated DEG pairs were obtained after integration with the identified methylated DEGs through a ceRNA mechanism ( Figure 6A). Finally, a circRNA-miRNA-mRNA regulatory network involving 1 downregulated circRNA, 1 upregulated miRNA and 20 downregulated mRNAs was visualized using Cytoscape software 3.7 ( Figure 6B). The results from the ceRNA network showed that hsa_circ_0001753 was derived from exons 12-18 of the CREB3L2 gene at chr7: 137,597,736-137,600,758 and eventually formed a mature sequence of 264 nt. The CSCD database showed that hsa_circ_0001753 contained the structure of miRNA binding sites (MREs) and RNA binding protein (RBP), suggesting that hsa_circ_0001753 is a key circRNA that potentially regulates OS development through spongy miRNA. Then, we explored the potential targeting relationship between hsa_circ_0001753 and hsa_miR_760 from the online database CSCD ( Figure 6C).

PPI Network Construction and Hub Gene Identification
We constructed the 20 mRNAs of ceRNA network expression proteins using the STRING online tool to construct PPI networks. When the mRNAs with required confidence (combined score) greater than 0.15 were submitted to STRING, a total of 12 PPI interactions were obtained, and the results are shown in Figure 7. The hub genes in the networks with a degree of connectivity greater than 2 were identified. A total of 7 hub genes were screened from the PPI network, which included CD74, FRZB, FSTL1, GDF5, MYO1C, PRELP and SEMA3F. The details of the hub genes are listed in Table 5.

Survival Analysis of the Hub Genes
We attempted to investigate Kaplan-Meier curves of the 7 hub genes, and the results revealed that the overall survival time of OS patients with low expression of CD74 was lower than that of patients with high expression (p < 0.05; Figure 8A). However, the effects of FRZB, FSTL1, GDF5, MYO1C, PRELP, and SEMA3F expression on OS were not statistically significant (p > 0.05; Figures 8B-G). To determine the accurate threshold for predicting OS by CD74, an ROC curve was constructed. The ROC curve of CD74 is shown in Figure 8H, and the area under the curve (AUC) was 1 year 0.674 and 3 years 0.673, which represented good diagnostic ability for OS.

Validation of hsa_circ_0001753/ has_miR_760/CD74 in Cells and Tissues
We constructed the hsa_circ_0001753/has_miR_760/CD74 network based on the ceRNA mechanism and the survival analysis of CD74. To test the reliability and stability of the bioinformatics analysis results, the expression patterns of this network were validated in osteoblast and OS cell lines by RT-PCR. Consistent with the bioinformatics analysis results, hsa_circ_0001753 was significantly downregulated, hsa_miR_760 was significantly upregulated and CD74 was significantly downregulated in the osteosarcoma cell lines  compared to the osteoblast cell line (Figures 9A-C). This expression trend revealed that this regulatory network of hsa_circ_0001753/has_miR_760/CD74 conforms to the ceRNA mechanism. Subsequently, we utilized human tissue microarrays (TMAs) containing 70 OS samples and 20 normal bone tissue samples tissues to determine the expression of CD74 by IHC. CD74 is relatively low expressed in osteosarcoma tissue as compared to normal bone tissue ( Figures 9D,E); of the 70 osteosarcoma tissue samples, 73% had low expression of CD74 ( Figure 9F).

The DNA Methylation Pattern of CD74
The promoter of the CD74 gene is generally considered as a 2,000 bp sequence upstream of the exon 1 site. The CpG island of the CD74 promoter was predicted using MethPrimer online software. We found that there was 1 CpG island located from 271 to 383 with 9 CpG sites in the promoter region of the CD74 gene, as shown in Figure 10A. Figure 10B shows that the BSP primer targeted the CpG island in the promoter region of the CD74 gene, and the CG dinucleotides were marked in red. DNA methylation status of the CD74 promoter were analyzed by BSP. The methylation status of the promoter region of the CD74 gene is shown in Figure 10C, and each row indicates the sequence of an individual clone. Black circles represent methylated CpG sites and white circles represent unmethylated CpG sites. The methylation rate of CD74 promoter DNA in osteoblasts was 96.7%, compared to 82.2% for U2OS, 91.1% for MG63 and 98.9% for HOS in the osteosarcoma group.

DISCUSSION
OS is the most common malignant bone tumour and seriously impacts the quality of life and life expectancy of patients. For osteosarcoma patients, early diagnosis and treatment of the disease is of great significance. Therefore, there is important clinical value and great market prospects for further exploring the pathogenesis of OS, finding possible therapeutic targets and realizing early diagnosis, targeted therapy and individualized treatments. Recently, microarray and high-throughput sequencing technology have been widely used in gene research on different diseases, bringing new hope for early detection, diagnosis and targeted therapy of tumours.
In our study, based on microarray datasets and using bioinformatics analysis, we identified 171 DECs, 81 DEMs, 3025 DMSs and 2,697 DEGs in the OS group compared to the control group. Using PCA, we found that the OS group was clearly separated from the control group, indicating that the screened DECs, DEMs, DMSs and DEGs were specific and could be used to identify relevant genes implicated in OS development.
To further investigate the biological pathways, we performed GSEA using the expression matrix of DEGs. The GSEA findings showed that DEGs are mainly involved in DNA replication, mitosis and cell cycle regulation signalling pathways. Genome instability is a sign of cancer, and DNA replication is the cellular process that most easily leads to cancer (Gaillard et al., 2015). Abnormal mitosis can promote the unlimited growth of defective cells, thus becoming the main method of tumour development (Mc Gee, 2015). Cell cycle disorder is the basis of abnormal cell proliferation, which is an important feature of cancer, and loss of cell cycle checkpoint control will promote genetic instability (Williams and Stoeber, 2012). Thus, the present study demonstrated that DEGs participate in vital cancer-related biological processes in OS. GO analysis showed that the methylated DEGs were mainly enriched in terms closely related to osteosarcoma, such as cartilage development, transforming growth factor beta, signal transduction by the p53 class mediator, extracellular matrix structural constituents, extracellular exosomes and chemokine receptors. It has been reported that in addition to a large number of osteoblasts, there are also many chondrocytes formed by cartilage, indicating the important role of cartilage development in the progression of OS (Lu et al., 2014). Transforming growth factor beta can promote angiogenesis, bone remodelling, and cell migration and inhibit immune surveillance, which is beneficial to the spread and metastasis of primary tumours and plays a role in promoting tumours in OS (Verrecchia and Rédini, 2018).
p53 is a nuclear protein related to malignant transformation in several tumour model systems, and p53 gene mutations have also been found in human osteosarcoma cell lines (Masuda et al., 1987). Structural changes in the extracellular matrix (ECM) are necessary for cell migration in the process of tissue remodelling, and its components play a key role in angiogenesis and are the basis of tumour invasion and metastasis (Roomi et al., 2006). It has been confirmed in the last 10 years that exosomes are The IHC staining scores were calculated in 70 OS tissues. IHC staining score < 2 was defined as low expression, while a score ≥ 2 was regards as high expression. Data shown are mean ± SD (n 3). (*p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001). Frontiers in Cell and Developmental Biology | www.frontiersin.org December 2021 | Volume 9 | Article 763347 15 extracellular vesicles released by cells that participate in cell signal crosstalk by transferring cellular components (RNA, proteins and lipids) across cells, and can induce changes in cell behaviour, playing an important role in the progression and metastasis of OS (Chicón-Bosch and Tirado, 2020). Chemokines and their receptors have been proven to be involved in the occurrence and development of malignant tumours. The functions of specific chemokines and their receptors are closely related to the biological environment and microenvironment of their expression, and they are considered active participants in osteosarcoma biology (von Luettichau et al., 2008). Another very interesting and meaningful research result is that KEGG pathway analysis showed that the methylated DEGs were also mainly enriched in the chemokine signalling pathway. Therefore, we speculate that methylated DEGs may be involved in the development of OS by affecting extracellular matrix components, cartilage development, transforming growth factors, exosomes, chemokine receptors, and signal transduction pathways of p53 mediators and may be a potential target for diagnosis and treatment.
Based on the ceRNA mechanism and methylated DEGs in OS, we constructed a circRNA-miRNA-mRNA regulatory network and identified the hub genes. The overall survival of CD74 with low expression levels in the 7 hub genes was shorter than that with high expression levels, while the area under the curve (AUC) of CD74 represented a good diagnostic capability for OS. We also validated the expression pattern of hsa_circ_0001753/ has_miR_760/CD74 in osteoblasts and OS cell lines, thus confirming that the ceRNA regulatory molecular network hsa_circ_0001753/has_miR_760/CD74 potentially plays important roles in the pathogenesis of OS. Compared with normal bone TMAs, the expression of CD74 in OS TMAs is lower. CD74 was hypomethylated in osteoblast and hypermethylated in osteosarcoma cells as compared to DNA methylation patterns in 19 osteosarcoma cell lines and 2 normal osteoblasts in the GSE36002 datasets. These results fully demonstrated that the expression level of hypermethylationmodified CD74 gene was decreased. Unfortunately, analysis of DNA methylation patterns in osteoblast and osteosarcoma cells using BSP did not reveal significantly different trend results. We suggested that this might be due to the fact that the type and number of osteoblastic and osteosarcoma cell lines we used were different from the GSE36002 datasets, so that no clear trend of difference was observed in our verification results. CD74, also known as the MHC class II-associated invariant chain, is a transmembrane glycoprotein that plays a vital role as a chaperone of MHC class II proteins during antigen presentation. CD74 is a cell surface receptor for macrophage migration inhibitory factor and is associated with tumour progression and metastasis. Its expression can be used as a prognostic factor for many cancers, and its high relative expression can be used as a marker of tumour progression (Gil-Yarom et al., 2017). Studies have shown that the overexpression of CD74 in thyroid cancer is associated with advanced tumour staging and can be used as a therapeutic target. Treatment of thyroid cancer cells with anti-CD74 antibodies inhibits cell growth, colony formation, cell migration and invasion, and secretion of vascular endothelial growth factors (Cheng et al., 2015). Another study revealed the mechanism of CD74 in the proliferation, invasion and angiogenesis of bladder cancer; indicating that it can be used as a potential therapeutic target for bladder cancer (Gai et al., 2018). Of interest, CD74 is upregulated in various forms of cancer and is associated with proliferative and metastatic potential but is downregulated in several cancers, including OS. High expression of CD74 in brain metastatic tumour cells can cause the processing of functional HLA class II cells, and is associated with a better prognosis (Zeiner et al., 2018). Another study confirmed that CD74 is a useful tumour cell prognostic protein marker closely related to recurrence-free survival and overall survival in stage III melanoma (Ekmekcioglu et al., 2016). We speculate that CD74 may also participate in the occurrence and development of OS by regulating malignant biological behaviours such as cell proliferation, migration, and invasion.
Before our study, there were also some related studies on the circRNA of OS. A study by Liu et al. (2017) on microarray expression profile and functional analysis of circRNA in osteosarcoma provided a good basis for our study. In their study, we first reported the comprehensive expression profile of circRNA in osteosarcoma and constructed a ceRNAs prediction network. However, on the basis of this study (GSE96964), we have performed a more in-depth analysis by combining miRNA (GSE70367), mRNA (GSE42351 and GSE33382) and DNA methylation data (GSE36002) from other osteosarcoma samples and constructed a new regulatory network of ceRNA. The large number of samples in multiple databases may provide more reliable research results. In addition, we verified the reliability of the regulatory network on cell lines, and conducted immunohistochemical verification on a large number of sample tissues. The hsa_circ_0001753/ has_miR_760/CD74 network constructed based on competing endogenous rna mechanisms and dna methodology may serve as a target for the early diagnosis and specific treatment of OS, and the findings can provide novel insights into the pathogenesis of OS.
Although rigorous bioinformatics analysis was adopted in this research and some results were achieved, it still has some limitations. Due to the small amount of cell sample used to detect the methylation pattern of CD74, and OS tumours are highly heterogeneous, our results need to be further verified with clinical samples to obtain more accurate results. To identify the potential molecular mechanism of OS, additional in vivo/in vitro experimental studies are necessary.

CONCLUSION
In summary, this study identified DECs, DEMs, DMSs and DEGs related to OS. Methylated DEGs were also mainly concentrated in pathways related to osteosarcoma. On the basis of the ceRNA mechanism and methylated DEGs, the hsa_circ_0001753/ has_miR_760/CD74 network was constructed and validated in cell lines and tissues. Low expression of CD74 is associated with poor overall survival time and presents a good diagnostic ability.
Frontiers in Cell and Developmental Biology | www.frontiersin.org December 2021 | Volume 9 | Article 763347 Methylated DEGs may be involved in the occurrence and development of OS, and hsa_circ_0001753/has_miR_760/ CD74 may be a target for early diagnosis of and specific treatment for OS. Using integrated bioinformatics analysis of integral competing endogenous RNA mechanisms and DNA methylation may provide valuable information for exploring the potential pathogenesis of OS.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

AUTHOR CONTRIBUTIONS
TW, HL, and WH conceived and designed the study. TW, TL,