Identification of microRNAs and target genes associated with the development of human intervertebral disc degeneration by bioinformatic analysis and verification

Micro RNAs (miRNAs) are widely recognized to play an essential role via target genes in the development of intervertebral disc degeneration( IDD), but the molecular mechanisms remain unclear. To identify the key microRNAs and potential targets during IDD, the Gene Expression Omnibus datasets (GSE19943, GSE63492, and GSE116726) were downloaded.An R package was used to identify differentially expressed miRNAs (DEMs) and four online tools(TargetScan, miRDB, miRTarBase, and DIANA-TarBase) were performed to predict their target genes. Functional enrichment analysis revealed that DEMs gene targets were highly enriched in cell development, cell differentiation, and the p53 and Wnt signaling pathways. we identified 13 hub genes with node degree≥10 through established a protein-protein interaction (PPI) network. Among them,MAPK8, BMP4, and GSK3B were top 3 highest degree. After constructing the miRNA-target gene-functional analysis network, we found that most hub genes could be regulated by miR-557 and were mainly enriched in cell development, cell differentiation, and Wnt signaling pathway. Further in vitro experiment by qRT-PCR confirmed that miR-577 was significantly downregulated than the control,whereas miR-516-3p was significantly upregulated. Together, the key microRNA and their target genes identified in this study help us understand the underlying pathogenesis mechanisms in the development of IDD, and provide diagnostic biomarkers and new therapeutic strategies for the treatment of IDD.


Introduction
The Search Tool for the Retrieval of Interacting Genes (STRING; http://string-db.org/) [20]database is a pre-computed global resource, which was designed to predict functional associations between proteins. These target genes were mapped to STRING to evaluate the PPI information and then visualized with the Cytoscape software (http://www.cytoscape.org/). A combined score > 0.4 was considered statistically significant. To screen the hub genes, node degree ≥ 10 was set as the cut-off criterion. Moreover, functional and pathway enrichment analysis for target genes was performed using DAVID to create miRNA-target genes functionalanalysis network. The interaction network was visualized using Cytoscape software.

Clinical Specimens
The human lumber NP specimens(3 normal discs and 3 degenerative discs) were obtained from patients with lumbar disc herniation undergoing spinal surgery at our hospital. The severe degree of enrolled patients with IDD was assessed according to the Pfirrmann classification (normal: grade 1-2, degeneration: grade 3-5) [21].This study was approved by the Ethics Committees of the Yangpu Hospital Affiliated to Tongji University and all patients written informed consent before operation.

Statistical analysis
All data are presented as the means ± SEM of three independent experiments. Statistical analysis were carried out using GraphPad prism 6.0 (GraphPad Software Inc.LaJolla,CA,USA).Differential expression levels of DEMs were compared by using Student's t test. A two-tailed value of P < 0.05 was considered as statistically significant.

GO and KEGG pathway enrichment analysis
We used GO analysis to investigate the underlying biological function of these target genes predicted by DEMs. The top 30 of the GO enrichment analysis of target genes were shown in Fig. 3A. In the BP category, these target genes were mainly enriched in cell development, cell differentiation, positive regulation of the transforming growth factor beta (TGF-beta) receptor signaling pathway, and Notch receptor processing. In the MF category, these target genes were mainly enriched in H4 histone acetyltransferase activity, transmembrane receptor protein tyrosine phosphatase activity, R-SMAD binding, and histone acetyltransferase binding. Finally, of the top 30 GO terms only Cul4−RING ubiquitin ligase complex related to a cellular component. The KEGG pathway analysis revealed that these target genes were enriched in the p53 signaling pathway, the Wnt signaling pathway and so on ( Fig.3B).

PPI network construction
To further explore the molecular mechanisms of IDD, we mapped the target genes of the two DEMs onto the PPI network using the STRING database. Based on the information from this database, 267 PPI pairs were obtained as shown in Fig. 4. Then target genes from the PPI network with a node degree ≥ 10 were selected as hub genes (Table 1). Of these genes, MAPK8 (degree, 23), BMP4 (degree, 20) and GSK3B(degree, 19) were the nodes with the highest degrees. This suggests that MAPK8, BMP4 and GSK3B may be key targets correlated with IDD.

miRNA -target genes functional analysis network construction and analysis
The miRNA -target genes functional analysis network was constructed with the cytoscape software ( Fig. 5A). We then investigated the potential interaction between miRNA and hub genes (Fig. 5B). The network showed that 9 hub genes (MAPK8, BMP4, GSK3B, CUL1, STAT5B, UBE3A, LEF1, PAX6, and RPS6KB1) predicted to be regulated bymiR-577,While 4 hub genes (PTPN11, CBL, SMURF1, and SIAH1) predicted to be regulated by miR-516a-3p. Notably, the 3 highest degree hub genes (MAPK8, BMP4, and GSK3B) were modulated by miR-577. This network also showed that these hub genes were highly correlated with cell development, cell differentiation, and were mainly involved in the Wnt signaling pathway. Taken together, our results indicated that miR-577 and miR-516a-3p may be participate inthe Wnt signaling pathway to regulate intervertebral disc cells development, and differentiation through targeting these hub genes.
Compared with the normal disc tissue samples, miR-577 was decreased significantly in IDD tissue samples, and miR-516a-3p expression was increased significantly.

Discussion
IDD is a major cause of low back pain and various degenerative spinal disorders. IDD has been a global health issue, which places a heavy burden on the healthcare system and results in severe In this study,3 miRNA microarray datasets were downloaded from GEO to obtain DEMs between normal and degenerative disc samples. Two DEMs, miR-557 and miR-516a-3p, were present in all the datasets. In the three microarray datasets, miR-557 was down-regulated. miR-516a-3p was downregulated in datasets. However, it was up-regulated in GSE116726 dataset. Our qRT-PCR data showed that miR-577 was significantly downregulated,and miR-516a-3p was significantly upregulated,which consistent with the GSE116726 dataset.MiR-557 is a recently identified miRNA. Previous studies have been reported that it is related to the progression of various human cancers, such as lung cancer, By constructing a PPI network, we identified 13 hub genes with degrees ≥ 10, including MAPK8, BMP4, and GSK3B, which were the 3 nodes with the highest degrees. To analyze the potential role of hub genes, an miRNA-target gene GO and KEGG analysis network was created and it was found that the hub genes, MAPK8, BMP4 and GSK3B, were mainly modulated by miR-577. Further analysis indicated that the hub genes were significantly associated with cell development, cell differentiation, and the Wnt signaling pathway. A literature search showed that the interaction among IDD and hub genes has rarely been reported. MAPK8, which is also known as JNK1, is a member of the MAP kinase and JNK family. Its functions include cell proliferation, cell differentiation, and apoptosis [38,39].Hua et al.found that MAPK8 might be candidate biomarker for apoptosis in Melanosis coli [40]. Wang et al. In conclusion, this study identified two novel DEMs (miR-557 and miR-516a-3p) and 13 hub genes via a comprehensive bioinformatics analysis approach. Validation experiments by qRT-PCR suggested that miR-557 and miR-516a-3p may have a certain influence on the pathological processes of IDD.
Notably, the roles of miR-557 through its targeting of MAPK8, BMP4 and GSK3B may improve our       The relative expression level of miR-516a-3p and miR-557 between normal disc tissues (n=3) and IDD tissue samples(n=3) by qRT-PCR. *P<0.05.

Supplementary Files
This is a list of supplementary files associated with this preprint. Click to download.