Integral Analyses of Competing Endogenous RNA Mechanisms and DNA Methylation Reveal Regulatory Mechanisms in Osteosarcoma
Tingrui Wu1† 
Bo Wei1† Hao Lin1 Boan Zhou1 
Tao Lin1 
Qianzheng Liu1 Hongxun Sang2 
Huan Liu3,4* 
Wenhua Huang1,5,6*- 1Orthopaedic Center, Affiliated Hospital of Guangdong Medical University, Zhanjiang, China
- 2Department of Orthopaedics, Shenzhen Hospital, Southern Medical University, Shenzhen, China
- 3Department of Orthopaedics, Affiliated Traditional Chinese Medicine Hospital, Southwest Medical University, Luzhou, China
- 4National Traditional Chinese Medicine Clinical Research Base, The Affiliated Traditional Chinese Medicine Hospital of Southwest Medical University, Luzhou, China
- 5Guangdong Provincial Key Laboratory of Medical Biomechanics, National Key Discipline of Human Anatomy, Guangdong Engineering Research Center for Translation of Medical 3D Printing Application, School of Basic Medical Sciences, Southern Medical University, Guangzhou, China
- 6Guangdong Medical Innovation Platform for Translation of 3D Printing Application, The Third Affiliated Hospital of Southern Medical University, Guangzhou, China
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 5-year 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 single-stranded 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 post-transcriptional 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; Li et al., 2020; 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 (Yu et al., 2019; He et al., 2020) and gastric cancer (Zhang et al., 2019; 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 (Li et al., 2019). 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 (Yang et al., 2020). 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.
FIGURE 1. Study flow chart.
Materials and Methods
Microarray Datasets
The Gene Expression Omnibus (GEO; http://www.ncbi.nlm.nih.gov/geo) is a public platform for storing genetic data. One circRNA expression profiling dataset (GSE96964 [platform: GPL19978 Agilent-69978 Arraystar Human CircRNA microarray V1]), one miRNA expression profiling dataset (GSE70367 [platform: GPL16384 Affymetrix Multispecies miRNA-3 Array]), one DNA methylation expression profiling dataset (GSE36002 [platform: GPL8490 Illumina HumanMethylation27 BeadChip]), and two gene expression profiling datasets (GSE42351 and GSE33382 [platform: 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.
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.
CeRNA Network Construction
In this study, we used the starBase database (Li et al., 2014) (http://starbase.sysu.edu.cn/) to predict the miRNA binding sites of DECs. Based on DEMs, we further screened target mRNAs. Then, the interactions between miRNAs and mRNAs were predicted by the miRanda (John et al., 2004), miRDB (Wong and Wang, 2015), miTarBase (Hsu et al., 2011), TargetScan (Agarwal et al., 2015) and miRMap (Vejnar et al., 2013) databases. The Cancer-Specific CircRNA Database (http://gb.whu.edu.cn/CSCD/) was used to understand the structural model of circRNAs, including miRNA binding sites (MREs), RNA binding proteins (RBPs), and open reading frames (ORFs). According to the circRNA-miRNA and mRNA-miRNA interactions, Cytoscape software 3.7 (Shannon et al., 2003) (https://cytoscape.org/) was used to visualize the ceRNA network of circRNA-miRNA-mRNA.
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.
Cell Culture
Human osteoblasts (hFOB1.19) and OS cell lines (U2OS, MG63, HOS) were cultured according to a previously reported protocol (Lin et al., 2020).
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 U6 expression by using an optimized comparative Ct (2−ΔΔCt) value method. Real-time PCRs were performed in triplicate.
The primers used were as follows: hsa_circ_0001753 (forward, 5′-AGGACTCGTTCCGTCTGTCACTC-3′ and reverse, 5′-AATGGTCTGGCACGAGGAATCAAC-3′), hsa_miR_760 (CGGCTCTGGGTCTGTGG), CD74 (forward, 5′-CTTTTCCATCCTGGTGACTCT-3′ and reverse, 5′-GAGGTGACTGTCAGTTTGTCC-3′), GAPDH (forward, 5′-CACCCACTCCTCCACCTTTG-3′ and reverse, 5′-CCACCACCCTGTTGCTGTAG-3′) and U6 (forward, 5′-TGGAACGCTTCACGAATTTGCG-3′ and reverse, 5′-GGAACGATACAGAGAAGATTAGC-3′). The size of the PCR product: hsa_circ_0001753: 112 bp/hsa_miR_760: 60 bp/CD74: 99 bp/GAPDH: 138 bp/U6: 68 bp.
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 bisulfite-treated 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.
Results
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.
FIGURE 2. Intragroup data repeatability test for the circRNA miRNA, methylated DNA and mRNA datasets through PCA and Pearson’s correlation analysis. (A–D) PCA of samples from the circRNA, miRNA, methylated DNA and mRNA datasets. The scatter plots are plotted on the x-, y-, and z-axes for PC1, PC2, and PC3, respectively, with each point representing a sample. The further the distance between the two samples was, the larger the difference in gene expression profiles of the two samples would be. (E–H) Pearson’s correlation analysis of the circRNA, miRNA, methylated DNA and mRNA datasets. The colour of heatmaps reflects the intensity of correlation. When 0 < correlation < 1, there is a positive correlation. When −1 < correlation < 0, there is a negative correlation. The greater the absolute value of numbers, the stronger the correlation. (I–L) Scatter diagrams showed Pearson’s correlation. PC1, principal component 1; PC2, principal component 2; PC3, principal component 3; PCA, principal component analysis.
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 downregulated); 81 differentially miRNAs (DEMs, 15 up- and 66 downregulated); 3,025 differentially methylated DNA sites (DMSs), covering 2,445 differentially methylated genes (DMGs, 1,069 hypermethylated and 1,376 hypomethylated) and 2,697 differentially expressed mRNAs (DEGs, 1,184 up- and 1,513 downregulated) were identified. The detailed results of differential expression are displayed in hierarchical cluster heatmaps and volcano plots (Figure 3).
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.
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.
FIGURE 4. The GSEA of DEGs datasets. The top portion of plots shows the enrichment scores for each gene, and the bottom portion shows the ranked genes. Y-axis: ranking metric, X-axis: individual ranks for all genes.
Screening for Methylated DEGs, GO and KEGG Pathway Analysis
Among the DEGs, a total of 95 hypermethylated downregulated mRNAs (hypermethylated DEGs) were screened from 1,069 overlapping hypermethylated mRNAs and 1,513 downregulated mRNAs, while 39 hypomethylated upregulated mRNAs (hypomethylated DEGs) were screened from 1,376 overlapping hypomethylated mRNAs and 1,184 upregulated mRNAs. The 134 methylated DEGs are presented in a Venn diagram (Figure 5A).
FIGURE 5. Identification and enrichment analysis of methylated DEGs. (A) A total of 95 hypermethylated DEGs and 39 hypomethylated DEGs were screened by Venn diagram; (B) GO analysis including CC, MF, and BP biological aspects; (C) BP: The top 15 enriched GO terms of the methylated DEGs; (D) BP: The methylated DEGs of the top 9 enriched GO terms; (E) CC: The top 15 enriched GO terms of the methylated DEGs; (F) CC: The methylated DEGs of the top 9 enriched GO terms; (G) MF: The top 15 enriched GO terms of the methylated DEGs; (H) MF: The methylated DEGs of the top 9 enriched GO terms; (I) KEGG: The top 15 enriched Reactome pathways of the methylated DEGs; (J) KEGG: The methylated DEGs of the top 9 enriched Reactome pathways.
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, 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.
TABLE 1. The top 15 enriched GO-BP terms of the methylated DEGs.
TABLE 2. The top 15 enriched GO-CC terms of the methylated DEGs.
TABLE 3. The top 15 enriched GO-MF terms of the methylated DEGs.
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.
TABLE 4. The top 15 enriched KEGG pathway terms of the methylated DEGs.
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 DEM-hypermethylated 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).
FIGURE 6. View of the circRNA-miRNA-mRNA network and circRNA structure chart. (A) Venn diagram of 152 upregulated DEM-hypermethylated DEG pairs and 284 downregulated DEM-hypomethylated DEG pairs. (B) The ceRNA network consists of 1 circRNA, 1 miRNA, and 20 mRNAs. (C) Has_circ_0001753 structure and hsa_circ_0001753 and hsa_miR_760 interactions identified by the CSCD database.
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.
FIGURE 7. Protein-protein interaction network constructed from methylated DEGs in the ceRNA network.
TABLE 5. Summaries for the function of hub genes.
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.
FIGURE 8. The overall survival Kaplan–Meier and ROC curves of 7 hub genes. Kaplan-Meier curves: (A) CD74, (B) FRZB, (C) FSTL1, (D) GDF5, (E) MYO1C, (F) PRELP, (G) SEMA3F. ROC curve: (H) CD74.
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).
FIGURE 9. Hsa_circ_0001753/has_miR_760/CD74 expression in cell lines and tissues. (A) hsa_circ_0001753 was downregulated in OS cell lines compared to osteoblasts. (B) hsa_miR_760 was upregulated in OS cell lines compared to osteoblasts. (C) CD74 was downregulated in OS cell lines compared to osteoblasts. (D) Representative images indicated the expression of CD74 in 90 normal bone and OS tissues detected with in IHC. (E) Dot distribution graph of CD74 IHC staining scores was shown in 20 normal bone tissues and 70 OS tissues. (F) 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).
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.
FIGURE 10. The DNA methylation pattern of CD74. (A) The CpG island was located between nucleotides 271 to 383 in the promoter region of CD74 gene, and there were 9 CpG sites in total. The short vertical line indicated CpG dinucleotides. (B) BSP primers targeting the CpG islands in the CD74 promoter region were designed and marked with red to show the CpG dinucleotides. (C) BSP primers targeting the CpG islands in the CD74 promoter region were designed and marked with red to show the CpG dinucleotides.
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 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 hypermethylation-modified 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. 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, and HL performed the bioinformatics analysis. BZ and QL performed experiments and data analysis. BW, HL, HS, and WH contributed to discussion and revision of the manuscript. All authors read and approved the final manuscript.
Funding
This work was supported by National Natural Science Foundation of China (61427807), Guangdong Basic and Applied Research Foundation (2020B151512001), Sichuan Applied Basic Research Project (2018JY0402), Luzhou municipal people’s government-southwest medical university science and technology strategic cooperation project (2018LZXNYD-ZK19), Natural Science Foundation of Guangdong Province Science and Technology Department (2020A1515010003), Sanming Project of Medicine in Shenzhen (SZSM201612019), The Young Innovative Talents Project of Guangdong Higher Education Institutions (2021KQNCX023), “Peaking Plan” for the reconstruction of high-level hospital at Affiliated Hospital of Guangdong Medical University and Special Fund for High-level Talents (Shizhen Zhong Team) of the People’s Government of Luzhou-Southwestern Medical University.
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s Note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Acknowledgments
We thank all the people who contributed to this study.
References
Agarwal, V., Bell, G. W., Nam, J.-W., and Bartel, D. P. (2015). Predicting Effective microRNA Target Sites in Mammalian mRNAs. Elife 4, 4. Epub 2015/08/13PubMed PMID: 26267216; PubMed Central PMCID: PMCPMC4532895. doi:10.7554/eLife.05005
Ashburner, M., Ball, C. A., Blake, J. A., Botstein, D., Butler, H., Cherry, J. M., et al. (2000). Gene Ontology: Tool for the Unification of Biology. Nat. Genet. 25 (1), 25–29. Epub 2000/05/10PubMed PMID: 10802651; PubMed Central PMCID: PMCPMC3037419. doi:10.1038/75556
Berner, K., Johannesen, T. B., Berner, A., Haugland, H. K., Bjerkehagen, B., Bøhler, P. J., et al. (2015). Time-trends on Incidence and Survival in a Nationwide and Unselected Cohort of Patients with Skeletal Osteosarcoma. Acta Oncologica 54 (1), 25–33. Epub 2014/06/25PubMed PMID: 24957555; PubMed Central PMCID: PMCPMC4364276. doi:10.3109/0284186X.2014.923934
Chen, L.-L., and Yang, L. (2015). Regulation of circRNA Biogenesis. RNA Biol. 12 (4), 381–388. Epub 2015/03/10PubMed PMID: 25746834; PubMed Central PMCID: PMCPMC4615371. doi:10.1080/15476286.2015.1020271
Chen, X., Mao, R., Su, W., Yang, X., Geng, Q., Guo, C., et al. (2020). Circular RNA circHIPK3 Modulates Autophagy via MIR124-3p-STAT3-Prkaa/ampkα Signaling in STK11 Mutant Lung Cancer. Autophagy 16 (4), 659–671. Epub 2019/06/25PubMed PMID: 31232177; PubMed Central PMCID: PMCPMC7138221. doi:10.1080/15548627.2019.1634945
Cheng, S.-P., Liu, C.-L., Chen, M.-J., Chien, M.-N., Leung, C.-H., Lin, C.-H., et al. (2015). CD74 Expression and its Therapeutic Potential in Thyroid Carcinoma. Endocr. Relat. Cancer 22 (2), 179–190. Epub 2015/01/21PubMed PMID: 25600560. doi:10.1530/erc-14-0269
Cheng, Z., Yu, C., Cui, S., Wang, H., Jin, H., Wang, C., et al. (2019). circTP63 Functions as a ceRNA to Promote Lung Squamous Cell Carcinoma Progression by Upregulating FOXM1. Nat. Commun. 10 (1), 3200. Epub 2019/07/22PubMed PMID: 31324812; PubMed Central PMCID: PMCPMC6642174. doi:10.1038/s41467-019-11162-4
Chicón-Bosch, M., and Tirado, O. M. (2020). Exosomes in Bone Sarcomas: Key Players in Metastasis. Cells 9 (1), 241. Epub 2020/01/23PubMed PMID: 31963599; PubMed Central PMCID: PMCPMC7016778. doi:10.3390/cells9010241
Deng, Z., Li, W., Alahdal, M., Zhang, N., Xie, J., Hu, X., et al. (2019). Overexpression of ClC-3 Chloride Channel in Chondrosarcoma: An In Vivo Immunohistochemical Study with Tissue Microarray. Med. Sci. Monit. 25, 5044–5053. Epub 2019/07/10PubMed PMID: 31281178; PubMed Central PMCID: PMCPMC6637820. doi:10.12659/msm.917382
Ding, W., Chen, G., and Shi, T. (2019). Integrative Analysis Identifies Potential DNA Methylation Biomarkers for Pan-Cancer Diagnosis and Prognosis. Epigenetics 14 (1), 67–80. Epub 2019/01/31PubMed PMID: 30696380; PubMed Central PMCID: PMCPMC6380428. doi:10.1080/15592294.2019.1568178
Ekmekcioglu, S., Davies, M. A., Tanese, K., Roszik, J., Shin-Sim, M., Bassett, R. L., et al. (2016). Inflammatory Marker Testing Identifies CD74 Expression in Melanoma Tumor Cells, and its Expression Associates with Favorable Survival for Stage III Melanoma. Clin. Cancer Res. 22 (12), 3016–3024. Epub 2016/01/20PubMed PMID: 26783288; PubMed Central PMCID: PMCPMC4911309. doi:10.1158/1078-0432.Ccr-15-2226
Gai, J. W., Wahafu, W., Song, L., Ping, H., Wang, M., Yang, F., et al. (2018). Expression of CD74 in Bladder Cancer and its Suppression in Association with Cancer Proliferation, Invasion and Angiogenesis in HT-1376 Cells. Oncol. Lett. 15 (5), 7631–7638. Epub 2018/05/08PubMed PMID: 29731899; PubMed Central PMCID: PMCPMC5920967. doi:10.3892/ol.2018.8309
Gaillard, H., García-Muse, T., and Aguilera, A. (2015). Replication Stress and Cancer. Nat. Rev. Cancer 15 (5), 276–289. Epub 2015/04/25PubMed PMID: 25907220. doi:10.1038/nrc3916
Gil-Yarom, N., Radomir, L., Sever, L., Kramer, M. P., Lewinsky, H., Bornstein, C., et al. (2017). CD74 Is a Novel Transcription Regulator. Proc. Natl. Acad. Sci. USA 114 (3), 562–567. Epub 2016/12/30PubMed PMID: 28031488; PubMed Central PMCID: PMCPMC5255621. doi:10.1073/pnas.1612195114
He, S., Guo, Z., Kang, Q., Wang, X., and Han, X. (2020). Circular RNA Hsa_circ_0000517 Modulates Hepatocellular Carcinoma Advancement via the miR-326/SMAD6 axis. Cancer Cel Int 20, 360. Epub 2020/08/11PubMed PMID: 32774154; PubMed Central PMCID: PMCPMC7397604. doi:10.1186/s12935-020-01447-w
He, X., and Zhang, J. (2006). Why Do Hubs Tend to Be Essential in Protein Networks? Plos Genet. 2 (6), e88. Epub 2006/06/06PubMed PMID: 16751849; PubMed Central PMCID: PMCPMC1473040. doi:10.1371/journal.pgen.0020088
Hsu, S.-D., Lin, F.-M., Wu, W.-Y., Liang, C., Huang, W.-C., Chan, W.-L., et al. (2011). miRTarBase: a Database Curates Experimentally Validated microRNA-Target Interactions. Nucleic Acids Res. 39, D163–D169. Database issueEpub 2010/11/13PubMed PMID: 21071411; PubMed Central PMCID: PMCPMC3013699. doi:10.1093/nar/gkq1107
John, B., Enright, A. J., Aravin, A., Tuschl, T., Sander, C., and Marks, D. S. (2004). Human MicroRNA Targets. Plos Biol. 2 (11), e363. Epub 2004/10/27PubMed PMID: 15502875; PubMed Central PMCID: PMCPMC521178. doi:10.1371/journal.pbio.0020363
Kanehisa, M., and Goto, S. (2000). KEGG: Kyoto Encyclopedia of Genes and Genomes. Nucleic Acids Res. 28 (1), 27–30. Epub 1999/12/11PubMed PMID: 10592173; PubMed Central PMCID: PMCPMC102409. doi:10.1093/nar/28.1.27
Li, J.-H., Liu, S., Zhou, H., Qu, L.-H., and Yang, J.-H. (2014). starBase v2.0: Decoding miRNA-ceRNA, miRNA-ncRNA and Protein-RNA Interaction Networks from Large-Scale CLIP-Seq Data. Nucl. Acids Res. 42, D92–D97. Database issueEpub 2013/12/04PubMed PMID: 24297251; PubMed Central PMCID: PMCPMC3964941. doi:10.1093/nar/gkt1248
Li, L., Xu, C., Liu, P., and Huang, J. (2017). Correlation Study of DNA Methylation of WNT6 Gene with Osteosarcoma in Children. Oncol. Lett. 14 (1), 271–275. Epub 2017/07/12PubMed PMID: 28693164; PubMed Central PMCID: PMCPMC5494854. doi:10.3892/ol.2017.6135
Li, R., Jiang, J., Shi, H., Qian, H., Zhang, X., and Xu, W. (2020). CircRNA: a Rising star in Gastric Cancer. Cell. Mol. Life Sci. 77 (9), 1661–1680. Epub 2019/10/30PubMed PMID: 31659415. doi:10.1007/s00018-019-03345-5
Li, S., Pei, Y., Wang, W., Liu, F., Zheng, K., and Zhang, X. (2019). Circular RNA 0001785 Regulates the Pathogenesis of Osteosarcoma as a ceRNA by Sponging miR-1200 to Upregulate HOXB2. Cell Cycle 18 (11), 1281–1291. Epub 2019/05/23PubMed PMID: 31116090; PubMed Central PMCID: PMCPMC6592237. doi:10.1080/15384101.2019.1618127
Li, X., Yang, L., and Chen, L.-L. (2018). The Biogenesis, Functions, and Challenges of Circular RNAs. Mol. Cel 71 (3), 428–442. Epub 2018/07/31PubMed PMID: 30057200. doi:10.1016/j.molcel.2018.06.034
Lin, H., Wu, T., Peng, L., Su, W., Wang, Y., Li, X., et al. (2020). Lnc-MAP6-1:3 Knockdown Inhibits Osteosarcoma Progression by Modulating Bax/Bcl-2 and Wnt/β-Catenin Pathways. Int. J. Med. Sci. 17 (15), 2248–2256. Epub 2020/09/15PubMed PMID: 32922188; PubMed Central PMCID: PMCPMC7484643. doi:10.7150/ijms.47405
Lin, R.-K., and Wang, Y.-C. (2014). Dysregulated Transcriptional and post-translational Control of DNA Methyltransferases in Cancer. Cell Biosci 4, 46. Epub 2014/01/01PubMed PMID: 25949795; PubMed Central PMCID: PMCPMC4422219. doi:10.1186/2045-3701-4-46
Liu, W., Zhang, J., Zou, C., Xie, X., Wang, Y., Wang, B., et al. (2017). Microarray Expression Profile and Functional Analysis of Circular RNAs in Osteosarcoma. Cell Physiol Biochem 43 (3), 969–985. Epub 2017/09/29PubMed PMID: 28957794. doi:10.1159/000481650
Lu, Y., Gitelis, S., Lei, G., Ding, M., Maki, C., Mira, R. R., et al. (2014). Research Findings Working with the P53 and Rb1 Targeted Osteosarcoma Mouse Model. Am. J. Cancer Res. 4 (3), 234–244. Epub 2014/06/25. PubMed PMID: 24959378; PubMed Central PMCID: PMCPMC4065404.
Luo, Z., Rong, Z., Zhang, J., Zhu, Z., Yu, Z., Li, T., et al. (2020). Circular RNA circCCDC9 Acts as a miR-6792-3p Sponge to Suppress the Progression of Gastric Cancer through Regulating CAV1 Expression. Mol. Cancer 19 (1), 86. Epub 2020/05/11PubMed PMID: 32386516; PubMed Central PMCID: PMCPMC7210689. doi:10.1186/s12943-020-01203-8
Masuda, H., Miller, C., Koeffler, H. P., Battifora, H., and Cline, M. J. (1987). Rearrangement of the P53 Gene in Human Osteogenic Sarcomas. Proc. Natl. Acad. Sci. 84 (21), 7716–7719. Epub 1987/11/01PubMed PMID: 2823272; PubMed Central PMCID: PMCPMC299371. doi:10.1073/pnas.84.21.7716
Mc Gee, M. M. (2015). Targeting the Mitotic Catastrophe Signaling Pathway in Cancer. Mediators Inflamm. 2015, 1–13. Epub 2015/10/23PubMed PMID: 26491220; PubMed Central PMCID: PMCPMC4600505. doi:10.1155/2015/146282
Meazza, C., and Scanagatta, P. (2016). Metastatic Osteosarcoma: a Challenging Multidisciplinary Treatment. Expert Rev. Anticancer Ther. 16 (5), 543–556. Epub 2016/03/22PubMed PMID: 26999418. doi:10.1586/14737140.2016.1168697
Mering, C. v., Huynen, M., Jaeggi, D., Schmidt, S., Bork, P., and Snel, B. (2003). STRING: a Database of Predicted Functional Associations between Proteins. Nucleic Acids Res. 31 (1), 258–261. Epub 2003/01/10PubMed PMID: 12519996; PubMed Central PMCID: PMCPMC165481. doi:10.1093/nar/gkg034
Mirabello, L., Troisi, R. J., and Savage, S. A. (2009). International Osteosarcoma Incidence Patterns in Children and Adolescents, Middle Ages and Elderly Persons. Int. J. Cancer 125 (1), 229–234. Epub 2009/03/31PubMed PMID: 19330840; PubMed Central PMCID: PMCPMC3048853. doi:10.1002/ijc.24320
Ritchie, M. E., Phipson, B., Wu, D., Hu, Y., Law, C. W., Shi, W., et al. (2015). Limma powers Differential Expression Analyses for RNA-Sequencing and Microarray Studies. Nucleic Acids Res. 43 (7), e47. Epub 2015/01/22PubMed PMID: 25605792; PubMed Central PMCID: PMCPMC4402510. doi:10.1093/nar/gkv007
Roomi, M. W., Ivanov, V., Kalinovsky, T., Niedzwiecki, A., and Rath, M. (2006). Effect of Ascorbic Acid, Lysine, Proline, and green tea Extract on Human Osteosarcoma Cell Line MNNG-HOS Xenografts in Nude Mice: Evaluation of Tumor Growth and Immunohistochemistry. Mo 23 (3), 411–418. Epub 2006/10/05PubMed PMID: 17018899. doi:10.1385/mo:23:3:411
Sang, Y., Chen, B., Song, X., Li, Y., Liang, Y., Han, D., et al. (2019). circRNA_0025202 Regulates Tamoxifen Sensitivity and Tumor Progression via Regulating the miR-182-5p/FOXO3a Axis in Breast Cancer. Mol. Ther. 27 (9), 1638–1652. Epub 2019/06/04PubMed PMID: 31153828; PubMed Central PMCID: PMCPMC6731174. doi:10.1016/j.ymthe.2019.05.011
Shannon, P., Markiel, A., Ozier, O., Baliga, N. S., Wang, J. T., Ramage, D., et al. (2003). Cytoscape: a Software Environment for Integrated Models of Biomolecular Interaction Networks. Genome Res. 13 (11), 2498–2504. Epub 2003/11/05PubMed PMID: 14597658; PubMed Central PMCID: PMCPMC403769. doi:10.1101/gr.1239303
Siegel, R. L., Miller, K. D., and Jemal, A. (2018). Cancer Statistics, 2018. CA: A Cancer J. Clinicians 68 (1), 7–30. Epub 2018/01/10PubMed PMID: 29313949. doi:10.3322/caac.21442
Skvortsova, K., Stirzaker, C., and Taberlay, P. (2019). The DNA Methylation Landscape in Cancer. Essays Biochem. 63 (6), 797–811. Epub 2019/12/18PubMed PMID: 31845735; PubMed Central PMCID: PMCPMC6923322. doi:10.1042/ebc20190037
Subramanian, A., Tamayo, P., Mootha, V. K., Mukherjee, S., Ebert, B. L., Gillette, M. A., et al. (2005). Gene Set Enrichment Analysis: a Knowledge-Based Approach for Interpreting Genome-wide Expression Profiles. Proc. Natl. Acad. Sci. 102 (43), 15545–15550. Epub 2005/10/04PubMed PMID: 16199517; PubMed Central PMCID: PMCPMC1239896. doi:10.1073/pnas.0506580102
Tay, Y., Rinn, J., and Pandolfi, P. P. (2014). The Multilayered Complexity of ceRNA Crosstalk and Competition. Nature 505 (7483), 344–352. Epub 2014/01/17PubMed PMID: 24429633; PubMed Central PMCID: PMCPMC4113481. doi:10.1038/nature12986
Vejnar, C. E., Blum, M., and Zdobnov, E. M. (2013). miRmap Web: Comprehensive microRNA Target Prediction Online. Nucleic Acids Res. 41, W165–W168. Web Server issue)Epub 2013/05/30PubMed PMID: 23716633; PubMed Central PMCID: PMCPMC3692044. doi:10.1093/nar/gkt430
Verrecchia, F., and Rédini, F. (2018). Transforming Growth Factor-β Signaling Plays a Pivotal Role in the Interplay between Osteosarcoma Cells and Their Microenvironment. Front. Oncol. 8, 133. Epub 2018/05/16PubMed PMID: 29761075; PubMed Central PMCID: PMCPMC5937053. doi:10.3389/fonc.2018.00133
von Luettichau, I., Segerer, S., Wechselberger, A., Notohamiprodjo, M., Nathrath, M., Kremer, M., et al. (2008). A Complex Pattern of Chemokine Receptor Expression Is Seen in Osteosarcoma. BMC Cancer 8, 23. Epub 2008/01/25PubMed PMID: 18215331; PubMed Central PMCID: PMCPMC2257957. doi:10.1186/1471-2407-8-23
Wang, T., Wang, J., Ren, W., Chen, S., Cheng, Y.-F., and Zhang, X.-M. (2020). CircRNA-0008717 Promotes Cell Proliferation, Migration, and Invasion by Regulating miR-203/Slug in Esophageal Cancer Cells. Ann. Transl Med. 8 (16), 999. Epub 2020/09/22PubMed PMID: 32953799; PubMed Central PMCID: PMCPMC7475474. doi:10.21037/atm-20-5205
Wang, Z., Deng, M., Chen, L., Wang, W., Liu, G., Liu, D., et al. (2020). Circular RNA Circ-03955 Promotes Epithelial-Mesenchymal Transition in Osteosarcoma by Regulating miR-3662/Metadherin Pathway. Front. Oncol. 10, 545460. Epub 2020/12/15PubMed PMID: 33312941; PubMed Central PMCID: PMCPMC7708376. doi:10.3389/fonc.2020.545460
Wang, Z., Yin, J., Zhou, W., Bai, J., Xie, Y., Xu, K., et al. (2020). Complex Impact of DNA Methylation on Transcriptional Dysregulation across 22 Human Cancer Types. Nucleic Acids Res. 48 (5), 2287–2302. Epub 2020/02/01PubMed PMID: 32002550; PubMed Central PMCID: PMCPMC7049702. doi:10.1093/nar/gkaa041
Williams, G. H., and Stoeber, K. (2012). The Cell Cycle and Cancer. J. Pathol. 226 (2), 352–364. Epub 2011/10/13PubMed PMID: 21990031. doi:10.1002/path.3022
Wong, N., and Wang, X. (2015). miRDB: an Online Resource for microRNA Target Prediction and Functional Annotations. Nucleic Acids Res. 43 (Database issue), D146–D152. Epub 2014/11/08PubMed PMID: 25378301; PubMed Central PMCID: PMCPMC4383922. doi:10.1093/nar/gku1104
Yang, P., Liu, Y., Qi, Y. C., and Lian, Z. H. (2020). High SENP3 Expression Promotes Cell Migration, Invasion, and Proliferation by Modulating DNA Methylation of E-Cadherin in Osteosarcoma. Technol. Cancer Res. Treat. 19, 153303382095698. Epub 2020/10/09PubMed PMID: 33030103; PubMed Central PMCID: PMCPMC7549150. doi:10.1177/1533033820956988
Yang, R., Xing, L., Zheng, X., Sun, Y., Wang, X., and Chen, J. (2019). The circRNA circAGFG1 Acts as a Sponge of miR-195-5p to Promote Triple-Negative Breast Cancer Progression through Regulating CCNE1 Expression. Mol. Cancer 18 (1), 4. Epub 2019/01/10PubMed PMID: 30621700; PubMed Central PMCID: PMCPMC6325825. doi:10.1186/s12943-018-0933-7
Yu, H., Zhao, L., Zhao, Y., Fei, J., and Zhang, W. (2020). Circular RNA Circ_0029589 Regulates Proliferation, Migration, Invasion, and Apoptosis in Ox-LDL-Stimulated VSMCs by Regulating miR-424-5p/IGF2 axis. Vasc. Pharmacol. 135, 106782. Epub 2020/08/30PubMed PMID: 32860985. doi:10.1016/j.vph.2020.106782
Yu, J., Yang, M., Zhou, B., Luo, J., Zhang, Z., Zhang, W., et al. (2019). CircRNA-104718 Acts as Competing Endogenous RNA and Promotes Hepatocellular Carcinoma Progression through microRNA-218-5p/TXNDC5 Signaling Pathway. Clin. Sci. (Lond) 133 (13), 1487–1503. Epub 2019/07/07PubMed PMID: 31278132. doi:10.1042/cs20190394
Zeiner, P. S., Zinke, J., Kowalewski, D. J., Bernatz, S., Tichy, J., Ronellenfitsch, M. W., et al. (2018). CD74 Regulates Complexity of Tumor Cell HLA Class II Peptidome in Brain Metastasis and Is a Positive Prognostic Marker for Patient Survival. Acta Neuropathol. Commun. 6 (1), 18. Epub 2018/03/02PubMed PMID: 29490700; PubMed Central PMCID: PMCPMC5831742. doi:10.1186/s40478-018-0521-5
Zhang, X., Wang, S., Wang, H., Cao, J., Huang, X., Chen, Z., et al. (2019). Circular RNA circNRIP1 Acts as a microRNA-149-5p Sponge to Promote Gastric Cancer Progression via the AKT1/mTOR Pathway. Mol. Cancer 18 (1), 20. Epub 2019/02/06PubMed PMID: 30717751; PubMed Central PMCID: PMCPMC6360801. doi:10.1186/s12943-018-0935-5
Keywords: bioinformatics, osteosarcoma, competing endogenous RNA, circRNA, methylation
Citation: Wu T, Wei B, Lin H, Zhou B, Lin T, Liu Q, Sang H, Liu H and Huang W (2021) Integral Analyses of Competing Endogenous RNA Mechanisms and DNA Methylation Reveal Regulatory Mechanisms in Osteosarcoma. Front. Cell Dev. Biol. 9:763347. doi: 10.3389/fcell.2021.763347
Received: 23 August 2021; Accepted: 26 November 2021;
Published: 09 December 2021.
Edited by:
Nejat Dalay, Istanbul University, TurkeyReviewed by:
Sulev Kõks, University of Tartu, EstoniaHui Li, Northeast Agricultural University, China
Copyright © 2021 Wu, Wei, Lin, Zhou, Lin, Liu, Sang, Liu and Huang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Huan Liu, 20016040@163.com; Wenhua Huang, huangwenhua2009@139.com
†These authors have contributed equally to this work