Integrated miRNA-mRNA network revealing the key molecular characteristics of ossification of the posterior longitudinal ligament.

BACKGROUND: Ossification of the posterior longitudinal ligament (OPLL) refers to an ectopic ossification disease originating from the posterior longitudinal ligament of the spine. Pressing on the spinal cord or nerve roots can cause limb sensory and motor disorders, significantly reducing the patient's quality of life. At present, the pathogenesis of OPLL is still unclear. The purpose of this study is to integrate microRNA (miRNA)-mRNA biological information data to further analyze the important molecules in the pathogenesis of OPLL, so as to provide targets for future OPLL molecular therapy.
METHODS: miRNA and mRNA expression profiles of GSE69787 were downloaded from Gene Expression Omnibus database and analyzed by edge R package. Funrich software was used to predict the target genes and transcription factors of de-miRNA. Gene ontology and Kyoto Encyclopedia of Genes and Genomes enrichment analysis of differentially expressed genes (DEGs) were carried out based on CLUEGO plug-in in Cytoscape. Using data collected from a search tool for the retrieval of interacting genes online database, a protein-protein interaction (PPI) network was constructed using Cytoscape. The hub gene selection and module analysis of PPI network were carried out by cytoHubba and molecular complex detection, plug-ins of Cytoscape software respectively.
RESULTS: A total of 346 genes, including 247 up-regulated genes and 99 down-regulated genes were selected as DEGs. SP1 was identified as an upstream transcription factor of de-miRNAs. Notably, gene ontology enrichment analysis shows that up- and down-regulated DEGs are mainly involved in BP, such as skeletal structure morphogenesis, skeletal system development, and animal organ morphogenesis. Kyoto Encyclopedia of Genes and Genomes enrichment analysis indicated that only WNT signaling pathway was associated with osteogenic differentiation. Lymphoid enhancer binding factor 1 and wingless-type MMTV integration site family member 2 Wingless-Type MMTV Integration site family member 2 were identified as hub genes, miR-520d-3p, miR-4782-3p, miR-6766-3p, and miR-199b-5p were identified as key miRNAs. In addition, 2 important network modules were obtained from PPI network.
CONCLUSIONS: In this study, we established a potential miRNA-mRNA regulatory network associated with OPLL, revealing the key molecular mechanism of OPLL and providing targets for future treatment or prevent its occurrence.

Introduction

Ossification of the posterior longitudinal ligament (OPLL) refers to the abnormal thickening of the posterior longitudinal ligament of the spine and the formation of bone tissue,[ about 70% occur in the cervical spine and more common in Asian populations.[ Most patients have no clinical symptoms, which are only seen by imaging examination. However, when calcified posterior longitudinal ligament compresses the spinal cord or nerve root, patients will have limb sensory and motor disorders, and even autonomic nerve dysfunction, which will greatly reduce the quality of life of patients.[ Surgery is needed if conservative treatment does not work,[ although now the improvement of surgical technique, but can’t seem to stop the progress of OPLL.[ In recent years, many medical researchers on the molecular structure characteristics were studied, the existence of a drug non-surgical treatment method and prevention of OPLL is now the focus of attention.[

MicroRNA (miRNA) is a kind of endogenous small RNA with a length of about 20 to 24 nucleotides, which has a variety of important and complex regulatory functions in cells.[ The effect of microRNA on the target gene mRNA mainly depends on the degree of complementarity with the target gene transcribed sequence, which can be expressed in three ways. The first is to cut off the mRNA molecule of the target gene. The miRNA binds to the target gene in a completely complementary manner, the action mode and function are very similar to siRNA (small interfering RNA), and finally cut off the target mRNA, resulting in degradation of the mRNA. The second is to inhibit the translation of target genes. During the action, it does not completely complement and bind to the target gene, thereby inhibiting translation without affecting the stability of the mRNA. The third is binding inhibition. It has the above two modes of action: when it is complementary to the target gene, it directly targets the cleaved mRNA; when it is not completely combined with the target gene, it plays a role in regulating gene expression.[ It has been reported that miRNA plays a key role in the occurrence and development of OPLL. For example, miR-10a regulates OPLL development by regulating the ID3 / RUNX2 axis.[ miR-563 can significantly promote osteogenic differentiation of posterior longitudinal ligament cells by down-regulating SMURF1 in vitro.[ There are significant differences in the expression of MiR-10a-5p, miR-563, and miR-210-3p in OPLL patients and normal population.[

In this study, we conducted multiple bioinformatics analysis methods on the data provided by Chen et al,[ with the purpose of revealing important molecular mechanisms in the pathogenesis of OPLL and providing a possible basis for future targeted and preventive treatment.

Materials and methods

Microarray data

The GSE69787 gene expression data set is derived from GPL17303 Ion Torrent Proton (Homo sapiens) (http://www.ncbi.nlm.nih.gov/geo/) which contains 12 samples in 2 groups of different RNA profiles, mRNA group and miRNA group. In mRNA group, 3 protein–protein interaction (PLL) cell specimens and 3 OPLL cell specimens were included. Meanwhile, there are 3 PLL cell specimens and 3 OPLL cell specimens in miRNA group. This study was based on data from an open database. Ethics and patient consent are not applicable.

Identification of differentially expressed miRNA and mRNA

The differential expression analysis of mRNA and miRNA was conducted with the edgeR package of R software, an excellent tool for statistical calculation and statistical mapping.[ The results of the analysis are saved in a file format, the cut-off criteria of differentially expressed genes (DEGs) is set as: FDR < 0.01, |log fold-change |> 2. Venn diagrams were produced by FunRich software, an independent software tool for functional enrichment and interaction network analysis of genes and proteins.[

Prediction potential transcription factors

The upstream transcription factors of DE-miRNAs were predicted using FunRich software. The screened DE-miRNAs were respectively up-regulated and down-regulated into FunRich, and the results showed only the top 10 predictive transcription factors.

Identification of DEGs

FunRich software was used to identify DEGs. We took the Venn intersection of the up-regulated DE-miRNAs targeted genes predicted by FunRich and the previously obtained down-regulated DE-mRNAs as the down-regulated DEGs. Similarly, the Venn intersection of the down-regulated DE-miRNAs targeted genes predicted by FunRich and the up-regulated DE-mRNAs as the up-regulated DEGs.

GO and KEGG enrichment analysis

DEGs’ gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis were performed using the plugin ClueGO[ (version 2.5.4) in the Cytoscape software[ (3.7.2) and the collated results are then imported into imageGP to adjust the bubble image visualization (http://www.ehbio.com/ImageGP). The standard setting with statistically significant difference was P < .05 and k value was 1.0.

PPI network construction and molecular complex detection (MCODE)

Integrate the PPI network. The interaction between DEGs was evaluated using the search tool for the retrieval of interacting genes database[ (version 11.0); a composite score > 0.4 was considered a statistically significant interaction. Subsequently, the analysis results of PPI network were loaded into Cytoscape software for visual adjustment. After that, we detected molecular complex of PPI network by using MCODE (version 1.5.1), which is a plugin of Cytoscape software.

Identification of hub genes and hub miRNAs

CytoHubba[ (version 0.1) is a plugin for identifying key genes from PPI network in Cytoscape. Considering the CytoHubba calculation ordering will count the genes for signaling pathways unrelated to osteogenic differentiation that would confuse our purpose and results. We took the top 10 genes after degree calculation method and the osteogenic differentiation pathway related genes to take Venn intersection, and the genes coexisting with the KEGG osteogenic differentiation related signaling pathway were considered as hub genes. Meanwhile, we took Venn intersection between target miRNAs of hub genes predicted by FunRich and De-miRNAs, the result was taken as hub miRNAs.

miRNA-mRNA network construction

Based on the acquisition of hub genes and hub miRNAs, as well as the differential expression result files analyzed by R software, Cytoscape was used to construct the mRNA-miRNA network according to their fold change values.

Results

Identification of DE-miRNAs, DE-mRNAs, and DEGs

According to the analysis results of edgeR package of R software, there were 344 de-miRNAs and 1631 de-mRNAs, among which 80 miRNAs were up-regulated and 264 miRNAs were down-regulated, 918 mRNAs were up-regulated and 713 mRNAs were down-regulated. According to FunRich prediction results, 2101 targeted genes were predicted by up-regulated de-miRNAs and 4338 targeted genes were predicted by down-regulated de-miRNAs, and 346 DEGs were obtained by Venn intersection, among which 247 DEGs were up-regulated and 99 DEGs were down-regulated (Fig. 1).

Figure 1

Screen of candidate differentially expressed genes. (A) The intersection of target genes of down-regulated DE-microRNAs and up-regulated DE-mRNAs; (B) The intersection of target genes of up-regulated DE-microRNAs and down-regulated DE-mRNAs; (C) The intersection of genes in osteogenesis related signaling pathways and top 10 genes with degree score.

Transcription factors of DE-miRNAs predicted

In our study, upstream transcription factors of DE-miRNAs were predicted by FunRich software, and the first 10 predicted results of up-regulated and down-regulated DE-miRNAs were shown in Figure 2, A and B, respectively. Interestingly, it is obvious that SP1 may regulate most DE-miRNAs (Fig. 2).

Figure 2

Predicted transcription factors of DE-microRNAs. (A) Top 10 transcription factors of up-regulated DE-microRNAs; (B) Top 10 transcription factors of down-regulated DE-miRNAs.

GO functional enrichment analysis

GO functional enrichment analysis of differentially expressed genes by “ClueGO” (Fig. 3). The results indicate that both up- and down-regulated DEGs are mainly enriched in biological processes. In Figure 3, A and B are the bubble maps of top 5 items of GO functional enrichment analysis of up-regulated and down-regulated genes in “ClueGO”, respectively (standard cut-off P < .05). As is evident from Figure 3, the related biological process of up-regulated DEGs are mainly enriched in skeletal system development, skeletal system morphogenesis, embryonic morphogenesis, regionalization, and embryonic organ morphogenesis. In cellular components analysis, the up-regulated DEGs are mainly enriched in lateral plasma membrane, voltage-gated potassium channel complex, potassium channel complex, presynaptic active zone cytoplasmic component, and cell cortex region. Besides, in molecular function (MF) analysis, up-regulated DEGs are mainly enriched in acetylgalactosaminyltransferase activity, UDP-glycosyltransferase activity, estrogen receptor binding, HMG box domain binding, voltage-gated potassium channel activity involved in cardiac muscle cell action potential repolarization. However, for down-regulated DEGs, there were no results about CC and MF enrichment analysis, and the main BP enrichment was tertiary alcohol biosynthetic process, cortisol biosynthetic process, aldosterone metabolic process, positive regulation of osteoblast differentiation, animal organ formation.

Figure 3

Gene ontology function enrichment bubble diagram. (A) Top 5 of up-regulate differentially expressed genes. (B) Top 5 of down-regulate differentially expressed genes (P < .05).

KEGG pathway enrichment analysis

To more intuitively show the enrichment analysis of DEGs in the KEGG pathways, we concentrated all DEGs together for analysis (P < .05) (Fig. 4). It can be clearly seen that 12 signaling pathways are enriched, and only one of the key pathways associated with osteogenic differentiation is WNT signaling pathway.

Figure 4

Kyoto Encyclopedia of Genes and Genomes enrichment analysis results of all differentially expressed genes (P < .05).

16 Identification of hub genes and hub miRNAs

The top 10 genes were gotten by the CytoHubba (Table 1). Next, we derived the data from the KEGG analysis in Cluego and used the relevant genes in the signaling pathways associated with osteogenic differentiation (Table 1) for the Venn crossover (Fig. 1). As a result, 2 hub genes were identified: wingless-type MMTV integration site family member 2 (Wnt2) and lymphoid enhancer binding factor 1 (LEF1) (Table 2). Based on these 2 hub genes, miR-520d-3p, miR-4782-3p, miR-6766-3p, and miR-199b-5p were identified as hub miRNAs through FunRich miRNA prediction function.

Table 1

Associated genes in the osteogenesis differentiation pathway and degree top 10 genes.

Table 2

Function of 2 hub genes.

PPI network construction and Module analysis

The search tool for the retrieval of interacting genes database is used for all DEGs PPI network construction. After that, the data was entered into Cytoscape to adjust the visualization, which has 266 nodes and 572 edges (Fig. 5). As can be seen in the figure, red represents up-regulated DEGs, blue represents down-regulated DEGs, and the color depth represents the absolute value of fold change. The larger the fold change value is, the darker the color is. Interestingly, the 2 rings in the middle were obtained by analysis of MCODE plug-in from PPI network. However, GO and KEGG enrichment analysis were further performed on these 2 modules (Table 3). Surprisingly, the genes in module 1 are mainly enriched in BP like positive regulation of osteoblast differentiation, positive regulation of animal organ morphogenesis, but KEGG just enriched in Basal cell carcinoma signaling pathway. For module 2, the genes are enriched in collagen fibril organization (BP), extracellular matrix component (CC), and extracellular matrix structural constituent conferring tensile strength (MF).

Figure 5

Protein–protein interaction network. The red circle represents up-regulated differentially expressed genes, while the blue circle represents down-regulated differentially expressed genes. The larger the absolute value of fold change is, the darker the color is, and the smaller the absolute value of fold change is, the lighter the color is. The middle 2 rings are module 1 and module 2 respectively.

Table 3

Modules analysis of the PPI network category.

Based on the screened hub miRNAs, we constructed the network interacting with miRNAs in DEGs (Fig. 6).

Figure 6

MicroRNA-mRNA network based on hub microRNAs.

Discussion

Ossification of the posterior longitudinal ligament of the spine is a degenerative disease. Calcification occurs in normal ligament tissues, which compress the spinal cord and nerve roots, leading to sensory and motor and other autonomic nerve dysfunction.[ This will greatly reduce the quality of life of patients.[ Currently, there is no other treatment except surgery, but surgical treatment is accompanied by the possibility of complications.[ In addition, the mechanism of OPLL occurrence and development is still poorly understood at present. Exploring the key molecular mechanism of OPLL occurrence provides the possibility of future molecularly targeted therapy, which will greatly improve the quality of life of patients, avoid the possibility of surgery, and even prevent the occurrence of OPLL.

Molecular targeted therapy is a new type of therapy. Targeted therapy drugs inhibit the occurrence of diseases and reverse the progression of diseases by inhibiting specific molecules or signaling pathways.[ Happily, with the progress of scientific research, more and more mechanisms of tumor occurrence and development have been discovered, and molecular targeted therapy has achieved good efficacy in the field of tumor. Such as targeted therapy in lung cancer[ and molecular targeted therapy in ovarian cancer.[

In our study, we identified 346 differential genes and 344 differentially expressed miRNAs, including 247 up-regulated DEGs and 99 down-regulated DEGs, 80 up-regulated DE-miRNAs and 264 down-regulated DE-miRNAs. SP1 has been identified as the upstream transcription factor of de-miRNAs. Interestingly, recent studies have shown that SP1 regulated mir-545-3p inactivates Wnt/β-catenin signaling by targeting LRP5, thereby acting as an osteo-suppressive factor.[ GO analysis showed that both up and down regulated DEGs are mainly enriched in biological processes like skeletal structure morphogenesis, skeletal system development and animal organ morphogenesis. KEGG enrichment analysis indicated that pathways related to osteogenic differentiation include MAPK signaling pathway[ and WNT signaling pathway,[ while the rest are non-osteogenic pathways such as pathways in cancer, breast cancer signaling pathway, arrhythmogenic right ventricular cardiomyopathy (ARVC) signaling pathway. However, a review of the relevant genes enriched in MAPK signaling pathway suggests that these DEGs-enriched MAPK signaling pathways are related to tumors and not involved in osteogenic differentiation,[ while the related genes of WNT signaling pathway are involved in osteogenic differentiation. So, we thought that the hub gene was in the WNT signaling pathway. Considering the cytoHubba computational ranking will count the genes of the other non-osteogenic pathways, which will confuse our research. Therefore, we took the top 10 genes in Degree and the genes in WNT signaling pathway to take Venn intersection, and identified 2 Hub genes: Wnt2 and LEF1. Besides, based on the 2 hub genes, we found in our de-miRNAs analysis results that there were 4 miRNAs correlated with them, namely miR-520d-3p, miR-4782-3p, miR-6766-3p, and miR-199b-5p, respectively. Strikingly, in our PPI complex module analysis, GO functional analysis of module 1's genes mainly enriched in the positive regulation of osteoblast differentiation and positive regulation of animal outraged morphogenesis, which are the osteogenesis related functions. However, the KEGG enrichment analysis has only enriched in Basal cell carcinoma, which is considered because the modular analysis is based on the whole PPI network, which includes other non-osteogenic genes, leading to confusion of our results. For module 2, GO functions are enriched in collagen fibril organization, extracellular matrix structural constituent conferring and tensile strength and extracellular matrix component, KEGG no enrichment results. Once again, consider that module 2 comes from the entire PPI network.

Interestingly, the GeneCards database indicates that Wnt2 and LEF1 play important roles in WNT signaling pathway. Wnt2 is a ligand for members of the frizzled family of 7 transmembrane receptors. Functions in the canonical Wnt signaling pathway that results in activation of transcription factors of the TCF/LEF family. Functions as upstream regulator of FGF10 expression. Besides, it plays an important role in embryonic lung development and may contribute to embryonic brain development by regulating the proliferation of dopaminergic precursors and neurons (by similarity).[ Lymphoid enhancer binding factor 1 (LEF1) participates in the Wnt signaling pathway. Activating transcription of target genes in the presence of CTNNB1 and EP300 and may plays a role in hair cell differentiation and follicle morphogenesis.[ Surprisingly, further investigation showed that Wnt2 and LEF1 regulated the osteogenic differentiation of mesenchymal stem cells through the WNT signaling pathway.[ As shown in Figure 7, SP1 transcription factors down-regulated miR-520d-3p, miR-4782-3p, and miR-6766-3p stimulated LEF1 overexpression, which led to OPLL through WNT signaling pathway. On the other hand, SP1 can also up-regulate mir-199b-5p inactivation of Wnt2, leading to OPLL through the WNT signaling pathway. These 2 mechanisms may work independently or cooperatively to cause OPLL.

Figure 7

The candidate transcription factor-microRNA-hub gene regulatory network in ossification of the posterior longitudinal ligament.

In this study, although an integrated bioinformatics analysis has been performed and a potential miRNA-mRNA regulatory network has been constructed, there are still some limitations. For example, the samples of gene expression omnibus data set are not large enough, and there is no relevant in vivo and in vitro experimental support that conforms to our conclusion. Although a large amount of evidence can be obtained from the literature review to support our conclusion, further experimental verification is still necessary.

Conclusion

In conclusion, we constructed a miRNA-mRNA regulatory network that may be involved in the pathogenesis of OPLL (Fig. 7). Through this study, we will further understand the potential pathogenesis of OPLL, provide targets for future OPLL therapy and improve patients’ quality of life.

Acknowledgments

We are thankful to Dr Xin Li Zhan (Spine and Osteopathy Ward, The First Affiliated Hospital of Guangxi Medical University) for his kindly assistance in all stages of the present study.

Author contributions

Conceptualization: Guoyong Xu, Chong Liu.

Data curation: Guoyong Xu, Tuo Liang, Zhaojie Qin

Formal analysis: Guoyong Xu, Chong Liu.

Funding acquisition: Guoyong Xu, Xinli Zhan.

Investigation: Guoyong Xu, Chong Liu, Chaojie Yu.

Methodology: Guoyong Xu, Chong Liu, Jie Jiang.

Project administration: Guoyong Xu, Chong Liu.

Resources: Guoyong Xu, Zide Zhang.

Software: Guoyong Xu, Tuo Liang, Jie Jiang, Xinli Zhan.

Supervision: Guoyong Xu, Xinli Zhan.

Validation: Guoyong Xu, Chong Liu, Zide Zhang, Jiarui Chen.

Visualization: Guoyong Xu, Tuo Liang, Zhaojie Qin, Chaojie Yu.

Writing – original draft: Guoyong Xu.

Writing – review & editing: Guoyong Xu, Chong Liu.