Construction of a transcription factor-miRNA-mRNA interactive network elucidates underlying pathogenesis for osteosarcoma and validation by qRT-PCR

Purpose: Osteosarcoma is characterized by features of rapid growth and early metastasis with a poor prognosis. The aim of our research is to investigate the potential transcription factor (TF)-miRNA-mRNA regulatory mechanism in osteosarcoma utilizing bioinformatics methods and validate by qRT-PCR. Methods: The microRNA (miRNA) expression profiling datasets (GSE28423 and GSE65071) and mRNA expression profiling dataset GSE33382 were collected from the Gene Expression Omnibus (GEO) database. Differentially expressed miRNAs (DEMs) and differentially expressed genes (DEGs) were screened using the limma package. Then, the TransmiR v2.0, miRDB, and Targetscan 7.2 database were applied for the acquisition of TF-miRNA and miRNA-mRNA interaction relationships, respectively. Finally, we built a TF-miRNA-mRNA interactive network. Furthermore, survival analysis was performed to identify sub-network with prognostic value and validate through qRT-PCR. Results: Eight overlapping DEMs and 682 DEGs were identified. Based on bioinformatics methods, 30 TF-miRNA interaction pairs and 25 miRNA-mRNA interaction pairs were screened. Finally, we constructed a TF-miRNA-mRNA regulatory network. Furthermore, laminin subunit gamma 1 (LAMC1) and thrombospondin-1 (THBS1), which involved in the network, were determined to have prognostic value and the corresponding subnetwork was identified. qRT-PCR results showed that LAMC1 mRNA expression was higher in osteosarcoma cells. Conclusion: Based on the survival analysis, a TF-miRNA–mRNA sub-network, that is TFs (SPI1, HEY1, and CEBPB)-hsa-miR-338-3p-target genes (LAMC1 and THBS1) was established. In conclusion, the construction of a potential TF-related regulatory network will help elucidate the underlying pathological mechanisms of osteosarcoma, and may provide novel insights for the diagnosis and treatment of osteosarcoma.


Introduction
Osteosarcoma is a typical representative of common malignant bone tumors, which mainly occur in the mesenchymal tissue of the metaphysis of long bones among children and adolescents. Rapid growth and early metastasis are typical clinical characteristics of osteosarcoma, approximately 11.4% of patients have proven metastasis at the time of diagnosis. [1] In the past 20 years, although comprehensive therapy, including extensive surgical resection combined with chemotherapy, has markedly improved overall survival rates, the 5-year survival rate for non-metastatic tumors remain at 60% to 75%, [2] whereas the metastatic osteosarcoma is only 29%. [1] Hence, it is necessary to investigate the precise mechanism and discover new therapeutic targets for osteosarcoma.
In recent years, transcription factors (TFs) and microRNAs (miRNAs) have been identified as important biomarkers for the diagnosis, targeted therapy and prognosis of osteosarcoma. TFs are a class of proteins that bind to specific DNA sequences and act as cancer suppressor genes or oncogenes. [3] TFs have been reported to perform important functions in modulating gene expression, which may contribute to oncogenic transformation. [4] Medicine miRNAs are a group of short non-coding RNAs that can bind to the 3' untranslated region (UTR) of targeted mRNAs, [5] and further modulate post-transcriptional translation by inhibiting or degrading the target genes. [6] Previous researches have shown that TFs can interact with miRNAs to regulate gene expression. [7] For instance, Salem et al [8] confirmed that miR-590-3p targets the transcription factor FOXA2 3'UTR region to reduce its expression, inhibit the negative regulation of the binding target gene VCAN, and increase cell proliferation, migration, and invasion ability. Additionally, Chandra Mangalhara et al [7] revealed the significant function of the extracellular signal-regulated kinase-2 (ERK2)-zinc finger E-box-binding homeobox 1 (ZEB1)-miR-101-1axis in regulating epithelial-mesenchymal transition (EMT) and cell migration in breast cancer. Mei et al [9] described the key role of the TFs/mir-23b/E26 transformation specific-1 (ETS1)/transcription factor 4 (TCF4) pathway in immune regulation of gastric adenocarcinoma. However, there are few detailed reports on the interactive regulatory network of TFs and miRNAs in osteosarcoma.
In this study, a range of bioinformatics tools were utilized to identify the potential interactive relationships among TFs, miR-NAs, and mRNAs. First, we analyzed differentially expressed miRNAs (DEMs) and differentially expressed genes (DEGs) from public databases. Subsequently, TransmiR v2.0, miRDB, and Targetscan 7.2 database were applied to screen out the upstream TF and downstream target genes of DEMs, respectively. Based on the above predicted outcomes, the TF-miRNA-mRNA interactive network was finally established, which aimed to explore the underlying pathogenesis and provide new insights for the diagnosis and treatment of osteosarcoma.

Acquisition of datasets
The miRNA expression profiling datasets (GSE28423 and GSE65071) and mRNA expression profiling dataset GSE33382 from the Gene Expression Omnibus (GEO) database (http:// www.ncbi.nlm.nih.gov/gds/) were collected in our study. Basic information on the expression profiling datasets from public database were listed in Table 1. Additionally, since the GEO datasets lack clinical data, 85 cases of osteosarcoma in the TARGET dataset with matrix expression and clinical data obtained from the UCSC Xena database (https://xena.ucsc.edu/) were used for subsequent survival analysis.

2.2.The Acquisition of DEMs and DEGs
The expression matrix of miRNA datasets (GSE28423 and GSE65071) and mRNA dataset GSE33382 were normalized by R package "limma". Differentially expressed miRNAs (DEMs) were obtained from GSE28423 and GSE65071 datasets, respectively. Next, the DEMs in these 2 datasets were uploaded to an online Venn diagram platform (http://bioinformatics.psb.ugent. be/webtools/Venn/) to acquire the shared DEMs. Differentially expressed mRNAs (DEGs) were obtained from GSE33382. In our study, |log fold-change (FC)|>1 and adjusted P < .05 were selecting as screening criteria to filter DEMs and DEGs.

Upstream transcription factors prediction for DEMs
TransmiR v2.0 database (http://www.cuilab.cn/transmir) stores the regulatory relationships between TFs and miRNAs, through which we could check the TF-miRNA regulation. In our study, the shared DEMs obtained in the previous step were uploaded to predict upstream TFs and acquire TF-miRNA interactions. Subsequently, the predicted TFs were intersected with DEGs to gain shared TFs and their interactions with miRNAs.

Downstream targeted genes prediction for DEMs
TargetScan 7.2 (http://www.targetscan.org/vert_71/) and miRDB (http://www.mirdb.org/) are the most commonly used online databases to analyze the regulations and functional annotations between miRNAs and target genes. In our study, shared DEMs were uploaded to predict downstream target genes and determine possible miRNA-mRNA regulation. Only those genes verified in both TargetScan 7.2 and miRDB databases were regarded as genuine candidate target genes. The target mRNAs were then overlapped with DEGs to obtain shared mRNAs and their interactions with miRNAs.

Functional enrichment analysis
Gene Ontology (GO) functional analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis are most commonly used to identify the biological functions of shared mRNAs. Functional enrichment analysis was performed utilizing R package "clusterProfiler". In our study, P < .05 was selecting as the screening criterion. The top five terms in the functional enrichment analysis were displayed.

Protein-protein interaction (PPI) network
The STRING online tool was utilized to establish a PPI network of shared mRNAs (http://string-db.org). The confidence score >0.15 was set for PPI pairs in our study. Finally, Cytoscape_ v3.8.0 was used to visualize the PPI network. The node number of each mRNA was calculated to screen for hub genes in the network. Those mRNAs with more than three nodes were selected as the hub genes.

Construction of a TF-miRNA-mRNA interactive network
According to the TF-miRNA and miRNA-mRNA interactions obtained in the previous steps, a TF-miRNA-mRNA interactive network was established and uploaded to Cytoscape_v3.8.0 for visualization.

Survival analysis to identify sub-network
Based on the TARGET dataset, survival analysis, including Cox regression, and Kaplan-Meier method, was utilized to screen out shared mRNAs with prognostic value and identify corresponding sub-network. P < .05 were selected as the screening criteria.

The acquisition of DEMs and DEGs
Based on the miRNA dataset, 223 DEMs (109 upregulated genes and 114 downregulated genes) in GSE28423 and 166 DEMs (78 upregulated genes and 88 downregulated genes) in GSE65071 were acquired. In the mRNA dataset GSE33382, 682 DEGs were screened. Among them, 300 mRNAs were significantly upregulated, whereas the other 382 mRNAs were downregulated. The volcano plot and heatmap of the DEMs and DEGs were shown in Figure 1.

Identification of transcription factors-miRNA Interactions
In the TransmiR v2.0 database, hsa-mir-219-5p had no corresponding upstream TF. Thus, the remaining 7 miRNAs were used in subsequent research. We identified 205 TFs for the remaining 7 miRNAs and 466 pairs of TF-miRNA interactions. After the intersection of 205 predicted TFs and 682 DEGs, 11 shared TFs and 30 pairs of TF-miRNA interactions were eventually obtained.

Identification of miRNA-mRNA interactions
The remaining seven miRNAs obtained in the previous step were uploaded to TargetScan 7.2 and miRDB databases, respectively. The predicted target genes from these 2 databases were then overlapped with DEGs obtained from the GSE33382 dataset. Ultimately, 23 shared mRNAs and 25 miRNA-mRNA interaction pairs were identified.

Functional enrichment analysis
To analyze the biological function of the 23 shared mRNA in the nosogenesis of osteosarcoma, GO, and KEGG enrichment analyses were implemented. For GO enrichment analysis, the negative regulation of cell migration was mainly concentrated in biological process (BP). The transcription regulator complex and collagen-containing extracellular matrix were primarily concentrated in cellular component (CC), while protein heterodimerization activity was mainly concentrated in molecular function (MF). KEGG pathway analysis indicated that the primary enrichment of these genes was concentrated in transcriptional misregulation in cancer (Fig. 3).

PPI network
According to the selecting criterion, those mRNA with more than three nodes were selectied as hub genes. In our research, 8 genes (MEF2C, thrombospondin-1 (THBS1), NT5E, GADD45A, BCL6, LAMC1, SIX4, and SLIT3) were regarded as hub genes, indicating that their roles are of great importance in the entire network (Fig. 4).

Construction of a TF-miRNA-mRNA network
Based on the above interactive relationship among TFs, miR-NAs, and mRNAs, a TF-miRNA-mRNA interactive network was finally built and visualized using Cytoscape_v3.8.0. The interactive network was illustrated in Figure 5.

Survival analysis to identify sub-network
Survival analysis, including Cox regression and Kaplan-Meier method, was applied to screen out shared mRNAs with prognostic value in the TARGET dataset. Eventually, two mRNAs, LAMC1, and THBS1, were considered to have significant prognostic value (Fig. 6A-C). Interestingly, these 2 genes were also considered as hub genes based on the PPI network. According to the survival analysis results, we found that high expression of LAMC1 and THBS1 was strongly related to worse prognosis in osteosarcoma. Based on survival analysis of shared mRNAs with prognostic value, the TF-miRNA-mRNA sub-network, that is TFs (SPI1, HEY1, and CEBPB)-hsa-miR-338-3p -target genes (LAMC1 and THBS1), was established.

of LAMC1 in osteosarcoma cells
Since THBS1 has been reported in osteosarcoma, [10] so we focused on LAMC1. We detected the expression of LAMC1 mRNA in OB, 143 B, HOS, SAOS, MG63, and U2 using qRT-PCR assay. When compared with OB cells, LAMC1 expression was increased in U2, HOS, and MG63 cells, especially in U2 cells (P < .001), but was decreased in 143 B and SAOS cells (Fig. 6D).

Discussion
Osteosarcoma, a highly prevalent malignant bone tumor, has typical clinical characteristics of rapid growth and early metastasis. Adverse clinical features are attributable to poor prognosis, and the exact etiology and pathogenesis are still unclear. Hence, thorough researches on the exact pathological mechanisms for osteosarcoma progression will be helpful in improving prognosis and preventing tumor recurrence. Gene chip technology has already been extensively applied in the research of various tumor types at present, which can screen key genes and therapeutic targets, and provide new insights for further clarifying the carcinogenic mechanisms and optimizing therapeutic strategies. Based on the gene expression profile data of osteosarcoma in the GEO database, 8 DEMs were overlapped from these 2 miRNA chips, which containing 5 upregulated and 3 downregulated genes. In recent years, many researches have documented that TFs can interact with miRNAs to regulate gene expression. Subsequently, a range of bioinformatics tools were used to predict the upstream TF and downstream target genes of DEMs. Generally, not all DEMs have their corresponding TFs. After a series of bioinformatics analyses, 11 TFs, 7 miRNAs, and 23 shared mRNAs, as well as their interactions, were considered to be closely associated with the nosogenesis and tumor progression of osteosarcoma. Eventually, a TF-miRNA-mRNA interactive network was built to investigate the underlying pathogenesis. Further survival analysis showed that TFs (SPI1, HEY1, and CEBPB)hsa-miR-338-3p -target genes (LAMC1 and THBS1) might play a critical part in osteosarcoma. Once the regulatory network was established, the potential roles of key genes were further investigated through functional enrichment analysis. Interestingly, KEGG pathway analysis indicated that the primary enrichment of these genes was concentrated in transcriptional misregulation in cancer, which further demonstrated that these genes were related to transcription factor disorder leading to oncogenesis. are essential for the regulation of target genes related to the malignant biological behavior of cancer cells. [11] In our study, 2 overlapping target genes (LAMC1 and THBS1) regulated  Figure A shows the PPI network of shared mRNAs in STRING. Figure B shows the identification of 8 hub genes. PPI = protein-protein interaction. by miR-338-3p were finally screened using survival analysis. Previous researches have shown that miR-338-3p was involved in the occurrence and progression of a variety of tumors, and acted as a tumor suppressive factor. [12][13][14][15][16][17][18] Xue et al [19] reported that miR-338-3p expression was obviously downregulated in colorectal cancer, and its overexpression significantly reduced cell invasion and migration in vitro. Additionally, miR-338-3p knockdown inhibits cell invasion by targeting smoothing genes in vitro, and could be a new underlying therapeutic target for hepatocellular carcinoma. [20] A recent study concluded that miR-338-3p was apparently downregulated in osteosarcoma, and inhibited cell proliferation, migration, and invasion, as well as epithelial-mesenchymal transition (EMT) in vitro, though targeted activator of 90 kDa heat shock protein ATPase homolog 1 (AHSA1). [21] LAMC1 and THBS1, 2 overlapping target genes regulated by miR-338-3p, have been documented in many studies. They are aberrantly expressed in a wide variety of tumors and strongly associated with the malignant biological behavior of tumors. LAMC1 encodes the laminin γ1 chain, which is essential for the basement membrane assembly. Previous reports have documented that laminin is crucial for tumor progression, invasion and metastasis. [22] Numerous studies have indicated that LAMC1 expression is obviously increased in a variety of tumor types, such as gastric carcinoma, [23] head and neck squamous cell carcinoma, [24] endometrial carcinoma, [25] glioma, [26] hepatocellular carcinoma, [27] and so on. Moreover, high expression of LAMC1 is closely related to poor prognosis, and could act as an underlying therapeutic target and prognostic indicator. Related researches have also been performed to reveal the pathological functions of LAMC1 in tumor progression. A study on esophageal squamous cell carcinoma reported by Fang et al [28] demonstrated that the upregulation of LAMC1 not only increased cell proliferation and migration ability in vitro, but also indirectly induces carcinogenesis by stimulating C-X-C motif Chemokine Ligand 1 (CXCL1) secretion to promote the formation of inflammatory cancer-associated fibroblasts (iCAF). Another study on gastric cancer reported by Han et al [23] suggested that LAMC1 knockdown constrained cell proliferation, migration, invation, and Warburg effect in vitro by deactivating the protein kinase B (AKT) and mitogen-activated protein kinase (MEK)/extracellular regulated protein kinase (ERK) pathway. Moreover, a study on hepatocellular carcinoma found that LAMC1 was strongly associated with tumor progression, possibly by regulating tumor-specific pyruvate kinase M2 (PKM2) expression via the phosphatase and tensin homologue deleted on chromosome 10 (PTEN)/AKT pathway. [29] To date, LAMC1 expression levels in osteosarcoma and its impact on clinical prognosis remain unclear. In our study, we indicated that LAMC1 is highly expressed in osteosarcoma cells, and high expression of LAMC1 predicts poor prognosis. THBS1 is the first member of thrombospondin (THBS) protein family to be discovered. It is a complex macromolecular glycoprotein with a series of different domains and functional domains. Each specific functional domain binds to different types of cytokines and exerts different biological effects. Previously, THBS1 was identified as a primary endogenous angiostatic factors that block tumor progression. [30,31] Whereas, recent studies have revealed its critical role in carcinogenesis. Borsotti et al [32] reported that THBS1 expression in primary melanoma and metastatic tumor was higher than that in normal and dysplastic nevus. Liu et al [33] confirmed that compared with non-cancerous tissues, THBS1 expression in colorectal carcinoma was obviously increased, which was closely related to liver metastasis and unfavorable prognosis of colorectal cancer, and further shown that THBS1 knockdown remarkably reduced cell migration and invasion ability in vitro. Daubon et al [34] verified that compared with low-grade gliomas, the expression of THBS1 in high-grade gliomas is higher, and THBS1 silencing can inhibit tumor cell proliferation and invasion alone or in combination with anti-angiogenic therapy. Recent research on osteosarcoma revealed that the expression of THBS1 in lung metastatic tumors is higher when compared with primary tumors in clinical tissues, and further functional studies have found that THBS1 promotes cell migration, invasion, and lung metastasis through a FAK-dependent pathway. [10] To conclude, we built an underlying TF-miRNA-mRNA interactive network to elucidate the underlying pathological mechanisms and provide novel insights for the diagnosis and treatment of osteosarcoma. However, there are still some limitations to our research. First, we only verified the expression level of LAMC1 in osteosarcoma cells, validation of expression level in the tumor tissue and further functional studies at the cellular level are still lacking. Moreover, another limitation is the lack of experimental validation of key regulatory network. Therefore, further experimental studies are needed to adress these limitations.