Identication of Key Genes and Pathways Associated With Pathogenesis of Intervertebral Disc Degeneration by Integrated miRNA-mRNA Network Analysis

Background: Intervertebral disc degeneration (IDD) is one of the most common cause of low back pain. Previous studies have suggested that miRNAs are associated with the pathogenesis of IDD. However, the underlying mechanisms remain unclear based on inconsistent results of available literatures. In addition, integrated miRNA-mRNA comprehensive analysis is limited. Material/Methods: In this study, we investigated the proles of differentially expressed miRNAs (DEMIs) and mRNAs (DEGs) and constructed a miRNA-mRNA regulatory network. First, transcription factors and target genes of DEMIs were predicted. Then, an intersection between DEMIs predicted genes and DEGs were performed to screen out the most signicant differential expressed common genes. Results: A total of 65 DEMIs and 61 common target genes were identied from datasets. Functional enrichment analysis showed that most genes were mainly involved in extracellular matrix organization and extracellular structure organization. Furthermore, DEGs were primarily enriched in PI3K−Akt signaling pathway, ECM−receptor interaction, focal adhesion and p53 signaling pathway, indicating that these pathways may be the critical pathways. Conclusions: In summary, several important miRNAs, as well as their related target genes and transcription factors in the pathogenesis of IDD were identied from our bioinformatic analysis, which may provide insights into underlying mechanisms and offer potential target genes for the treatment of IDD. and remodeling of collagen by annulus brosus cells of the intervertebral disc, that both hsa-miR-106b-5p and MMP-2 may serve as potential therapeutic target of IDD. Hsa-miR-98-5p, a downregulated miRNA which targets CHUK, COL4A1, COL4A2, GHR, GLRX, ITGB8, and THBS1 in the network, has that contributes to extracellular matrix degradation by targeting IL-6/STAT3 pathway in human intervertebral disc degeneration A preclinical that miR-141 IDD progression by interacting with SIRT1/NF-κB pathway and inhibition of miR-141 in may as a potential therapeutic approach in the of IDD downregulated in IDD mice while upregulation of miR-181a protected against inammatory response by inactivating the ERK pathway via suppression of tumor necrosis factor-related apoptosis-inducing ligand (TRAIL) in IDD mice result is consistent with nding of our miRNA-miRNA network. miR-29a has been demonstrated with protective effect by silencing the expression of MMP-2, inhibiting the brosis process, and reversing IDD in animal models via blocking the β-catenin translocation pathway from the cytoplasm to the nucleus [23] the LBP, low back pain; IVD, intervertebral disc; IDD, intervertebral disc degeneration; AF, annulus brosus; NP, nucleus pulposus; CEP, cartilage endplate; miRNA, microRNA; ECM, extracellular matrix; MMP, matrix metalloproteinase; DEGs, differentially expressed genes; BP, biological process; CC, cellular component; MF, molecular function; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Intervertebral disc degeneration (IDD) is considered as a main source of discogenic low back pain (LBP), which may cause disabilities in adults along with an enormous socioeconomic burden [1][2][3] . The incidence of LBP has risen dramatically with the changing of lifestyle [4] . Nowadays, LBP is highly recognized as a global health problem and it is estimated that over two-thirds of the population experience at least one episode of LBP, with roughly 10% turning to be chronically disabled [5] . Despite years of effort and intensive investigation, the currently available treatment options remain unsatisfactory. Several major theories have been put forward to explain the biology of IDD, including apoptosis [6] , in ammation [7] , aging [8,9] and biomechanical loading [10,11] . However, the genetic factors are regarded as the most signi cant contributor. Therefore, it is necessary to elucidate the possible molecular mechanism underlying IDD for developing effective treatments.
MicroRNAs are a class of non-coding single-stranded RNA of approximately 22 nucleotides in length encoded by endogenous genes. Increasing evidence has shown that miRNAs participate in RNA silencing and post-transcriptional regulation of gene expression. Recent reports have demonstrated that miRNAs are closely associated with pathogenesis of IDD and even act as a new therapeutic target or clinical biomarker for early diagnosis. For example, miR-24-3p leads to IDD by targeting insulin-like growth factor binding protein 5 and the ERK signaling pathway [12] . Inhibition of miR-129-5p results in IDD by promoting the apoptosis of nucleus pulposus cells via targeting BMP2 [13] .
To date, integrated miRNA-mRNA expression analysis regarding underlying mechanisms of IDD is still limited. The aim of this study is to investigate reliable miRNAs and their target genes related to IDD. Meanwhile, the study focuses on the functions of the target genes regulated by the differentially expressed miRNAs (DEMIs) during the pathogenesis of IDD, which may provide new insights into pathogenesis of IDD and potential target candidates for new therapy. offs to screen out DEG. Both heatmap and volcano plots were constructed to present expression pro les of differentially expressed genes, which was performed using R software.
2.3 Prediction of miRNA-targeted gene and DEMI-DEG regulatory network Both transcription factor (TF) and target genes of differentially expressed miRNAs were predicted using FunRich software (http://www.targetscan.org/), which is a powerful functional enrichment analysis tool that not only predicts TF and target genes of miRNAs but also provides GO enrichment analysis. The intersection of predicted target genes of DEMIs and DEGs were regarded as signi cantly differential expression target genes, which was performed by R software. These signi cant differential expression target genes and corresponding miRNAs were used to construct the miRNA-mRNA regulatory network using the Cytoscape software.
2.4 Functional enrichment analysis of target DEGs. GO is a major bioinformatics tool to annotate genes, analyze gene products and sequences to underlying biological phenomena, including biological process (BP), molecular function (MF), and cellular component (CC) [14] . Kyoto Encyclopedia of Genes and Genomes (KEGG) is a knowledge base for systematic analysis, annotation or visualization of gene functions and critical biological pathways closely related to intervertebral disc degeneration [15] . Both GO and KEGG enrichment analysis were conducted by Bioconductor. A P < 0.05 was considered to have statistical signi cance and to achieve signi cant enrichment.

Results
3.1 Differential Expression Analysis. For the expression pro les of miRNAs in IDD, a total of 414 DEMIs were identi ed from intervertebral disc samples in GSE116726. Among them, 189 were downregulated and 225 were upregulated within the P value < 0.01 and |log2FC| >4 criterion.
Meanwhile, for the expression pro les of genes in IDD, A total of 207 genes were identi ed to be signi cantly differently expressed (P < 0.05), including 116 downregulated and 91 upregulated DEGs. Heatmap plots and volcano plots of both DEMIs and DEGs are shown in Fig. 2 (A and B) and Fig. 2 (C and D) respectively, in which the red represents upregulated genes and the green represents downregulated genes.

3.2
The transcription factor analysis of the differential expressed miRNAs in IDD. We also analyzed the TF genes corresponding to the differential expressed miRNAs in IDD and compared the similarities and differences. The top 10 most signi cantly differential TFs included EGR1, SP1, SP4, 3.3 Differential expression miRNA predicted target genes (DEMIGs) and miRNA-mRNA regulatory network. Target genes of miRNA were predicted by FunRich software. In order to nd out the most signi cant differential expressed common genes, an intersection between DEMI predicted genes of GSE116726 and DEGs of GSE70362 were performed by R software and subsequently screened out 61 common genes, including 13 downregulated and 48 upregulated. These 61common DEGs and DEMIs were used to constructed miRNA-mRNA network, which consisted of 126 nodes and 218 edges, including 65 miRNAs (11 up-regulated and 54 down-regulated) (Fig. 4). The top 10 hub genes were identi ed by CytoHubba plugin using the Maximal Clique Centrality (MCC) method, including EDEM3, ITGB8, JARID2, GHR, THBS1, GATA6, hsa-miR-19b-3p, COL4A2, hsa-let-7f-5p, and CHUK.
To be mentioned, hub DEGs were all up-regulated while both hsa-miR-19b-3p and hsa-let-7f-5p were downregulated in the network. The relationship between miRNA and mRNA of the regulatory network was summarized in Table 1.

Discussion
IDD is a complex and multifactorial pathophysiological process which is characterized by excessive apoptosis of nucleus pulposus (NP) cells and over-degradation of extracellular matrix (ECM) components, leading to reduced hydration, loss of disc height, and decreased ability to absorb load [16] . To date, genetic factors have been considered as main contributor of IDD while the exact underlying mechanism has not been fully elucidated.
In the present study, we used mRNA and miRNA expression pro les to identify differential expression mRNAs (DEGs) and miRNAs (DEMIs). A total of 414 DEMIs and 207 DEGs were identi ed to be signi cantly differently expressed. In addition, the top 10 most signi cantly differential TFs included EGR1, SP1, SP4, POU2F1, MEF2A, NFIC, ZFP161, NKX6-1, FOXA1 and TCF3. Furthermore, 61 common genes via an intersection between DEMIs predicted genes and DEGs were used to build miRNA-mRNA regulatory network. Functional and pathway enrichment analysis of common DEGs were associated with extracellular matrix organization and PI3K−Akt signaling pathway.
Increasing evidence has shown that many cellular processes, including cell proliferation, apoptosis, and cytokine release, are regulated by miRNAs.
Recent studies further demonstrated that miRNA functioned by suppressing protein production from targeted genes which was closely associated with development of IDD [17] . Hsa-let-7b-5p is the most signi cant downregulated miRNA, which belongs to the let-7 family and plays an important regulatory role in the cell cycle and differentiation. Inhibition of let-7b-5p is associated with elevated expression of TNF-α, IL-1β, and IL-6, leading to activation of in ammation cascade [18] . Hsa-miR-29c-3p is the most signi cant upregulated miRNA. Recent study showed that overexpression of miR-29c-3p might inhibit proliferation and promote apoptosis and differentiation by inhibiting the expression of Akt3 [19] . We found that hsa-miR-106b-5p was downregulated while target gene MMP-2 was upregulated in our network. MMP-2 is a member of MMP family and MMP-2 mediates local degradation and remodeling of collagen by annulus brosus cells of the intervertebral disc, indicating that both hsa-miR-106b-5p and MMP-2 may serve as potential therapeutic target of IDD. Hsa-miR-98-5p, a downregulated miRNA which targets CHUK, COL4A1, COL4A2, EDEM3, GHR, GLRX, ITGB8, and THBS1 in the network, has been demonstrated that contributes to extracellular matrix degradation by targeting IL-6/STAT3 signaling pathway in human intervertebral disc degeneration [20] . A preclinical study showed that miR-141 promoted IDD progression by interacting with SIRT1/NF-κB pathway and inhibition of miR-141 in vivo may serve as a potential therapeutic approach in the treatment of IDD [21] . miR-181a-5p was downregulated in IDD mice while upregulation of miR-181a protected against in ammatory response by inactivating the ERK pathway via suppression of tumor necrosis factor-related apoptosis-inducing ligand (TRAIL) in IDD mice [22] . This result is consistent with nding of our miRNA-miRNA network. miR-29a has been demonstrated with protective effect by silencing the expression of MMP-2, inhibiting the brosis process, and reversing IDD in animal models via blocking the β-catenin translocation pathway from the cytoplasm to the nucleus [23] .
Transcription factors (TFs) regulate gene expression by activating or repressing target genes at the transcriptional level. Hence, we explore some important TFs related to miRNA target genes expression. SP1 is a ubiquitous zinc nger transcription factor that binds to GC-rich regions in the promoters of genes to activate transcription. Inhibition of Sp1 and Sp1 knockdown has been previously demonstrated to decrease expression of MMP3, ADAMTS4 and ADAMTS5 and subsequently suppress TNF-α-induced catabolic activity in nucleus pulposus cells [24] . EGR1 belongs to the early growth response family, with a wide range of activities as transcription factor ranging from regulation of cell growth to differentiation. EGR1 plays an important role that allows NP cells to adapt to anabolic or catabolic stimuli [25] . Both FOXA1 and FOXA2 are necessary for the formation of intervertebral disc and inhibition of FOXA1 or FOXA2 leads to aberrant differentiation of NP cells from notochord cells [26] .
In terms of GO enrichment analysis, we found that most of DEGs were mainly involved in extracellular matrix organization, extracellular structure organization, response to acid chemical, endoderm development, and cellular response to acid chemical. It has been reported that the increased expression of MMPs in NPCs leads to inadequate synthesis and excessive degradation of extracellular matrix (ECM), resulting in progression of IDD.
We found that MMP-2 was signi cantly up-regulated in the degenerative samples compared with normal samples in our analysis. Previous nding also revealed that the expression of MMP-2 was increased in intervertebral disc and MMP-2 was associated with local degradation and remodeling of collagen tissue [27] .
Of the signi cantly enriched pathways in the KEGG pathway analysis, PI3K/AKT signaling pathway signaling was of interest as it played an important role in protective effects on human nucleus pulposus under different pathological conditions based on literatures. The activation of PI3K/AKT has been showed protective effect against IDD by increase of ECM content, prevention of cell apoptosis and induction or prevention of cell autophagy. For example, under oxidative damage, Resveratrol increased ECM synthesis of NPCs by enhancing autophagy through PI3K/AKT pathway [28] . Further study con rmed that PI3K/AKT pathway had protective effects on human nucleus pulposus-derived mesenchymal stem cells under hypoxia and nutrition de ciency [29] . Conversely, PTEN, the only known lipid phosphatase counteracting the PI3K/AKT pathway, promoted intervertebral disc degeneration by regulating nucleus pulposus cell behaviors [30] . miR-21 exosomes prevented NPCs from apoptotic process and alleviated IVD degeneration via restraining PTEN and activating PI3K/AKT pathway in apoptotic NPCs [31] .
In conclusion, several important miRNAs, as well as their related target genes and transcription factors in the pathogenesis of IDD were identi ed from our bioinformatic analysis, which may provide insights into underlying mechanisms and offer potential target genes for the treatment of IDD.

Figure 1
Microarray data collection. Both miRNA (GSE116726) and mRNA (GSE70362) were downloaded from the National Center of Biotechnology Information (NCBI) Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) database. GSE116726 is a miRNA microarray containing 3 nucleus pulposus (NP) tissues from IDD patients and 3 normal NP tissues from fresh traumatic lumbar fracture patients. Meanwhile, mRNA expression pro le of NP tissue was extracted from GSE70362, including 14 non-degenerative NP samples and 10 degenerative NP samples. The whole analysis process is summarized in Figure 1.