Prediction of a Potential Mechanism of Intervertebral Disc Degeneration Based on a Novel Competitive Endogenous RNA Network

Low back pain which resulted from intervertebral disc degeneration (IDD) is a common health problem that afflicts people all over the world. Due to the lack of an overall understanding of the molecular interactions involved in IDD, we hope to better understand the pathogenetic mechanisms that drive the degenerative process. The purpose of this study is to obtain mRNAs, miRNAs, lncRNAs, and circRNAs associated with IDD gained from public databases and to establish an interaction network. According to the results of microarray analysis and bioinformatics analysis from the contrast of IDD and normal nucleus pulposus tissues, a total of 49 mRNAs, 10 miRNAs, 30 lncRNAs, and 4 circRNAs were obtained and a lncRNA/circRNA–miRNA–mRNA interaction network was constructed. NEAT1–miR-5100–COL10A1 and miR663AHG/HEIH/hsa-circ-0003600–miR-4741–HAS2/HYAL1/LYVE1 might be potential interaction axes of the molecular mechanism in IDD. The increased expression of NEAT1 might inhibit miR-5100 and subsequently upregulate the expression of COL10A1, which leads to IDD, while the increased expression of miR663AHG/HEIH/hsa-circ-0003600 might inhibit miR-4741 and indirectly upregulate HAS2/HYAL1/LYVE1, and leads to the protection from IDD. More interaction axes are to be exploited to provide theoretical bases for further study on IDD.


Introduction
Low back pain (LBP) is a common health problem that afflicts people all over the world, and it is estimated that around 84% of the world's population will experience low back pain in their lives [1]. Not only does it affect patients' quality of life, but it also has a broad socioeconomic impact, with some recent studies showing that related costs in the United States exceed $100 billion a year [2,3]. Importantly, intervertebral disc degeneration (IDD) is one of the major risk factors related to LBP [4]. It has been determined that 40% of chronic LBP is caused by IDD [5,6], and the main factors affecting the pathophysiology of IDD include genetic susceptibility, old age, smoking, alcohol addiction, obesity, and diabetes [7,8].
Currently, there is no specific treatment for IDD-induced LBP. The current treatment strategies of IDD-induced LBP focus only on pharmacological approaches to relieve pain or inflammation, followed by surgery as a last resort for severe disease stages. With the lack of long-term efficacy and increasing drug abuse, it is critical to develop a molecular targeted therapy for IDD, which may fundamentally prevent IDD or restore the function of the intervertebral disc to solve the problem [9,10].
IDD is a chronic pathological process with multiple complex etiologies [11]. The intervertebral disc is composed of the inner-most glycosaminoglycan, notochord-derived nucleus pulpous (NP), fibrocartilaginous annulus fibrosus, and superior and inferior cartilaginous end plates. As NP ages and degenerates, it gradually loses its water-bound matrix and becomes fibrotic and its ability to absorb mechanical loads is also diminished, which is accompanied by the phenotype changes in nucleus pulposus cells (NPC) [12]. NPC regulates the balance of synthesis/catabolism by secreting specific biological factors such as CTGF, Shh, and TGF-1, and the secretion changes of these biological factors may also lead to IDD [13]. Understanding the changes of NPC in IDD is essential for the pathogenesis and treatment of the disease.
So far, there have been many reports on the potential influence of IDD genes, such as COL1A1, COL9A2, COL9A3, COL11A2, IL-6, AGC1, VDR, and MMP-3 [14][15][16][17][18]. However, increasing evidences also show that noncoding RNAs (ncRNAs) including miRNA, lncRNA, and circRNA can exert effects on coding RNAs and influence biological processes such as cell proliferation and apoptosis. miRNAs act by binding to complementary sequences in the 3 ′ -untranslated region (UTR) of their target mRNA to trigger transcriptional inhibition or mRNA degradation [19]. lncRNAs and circRNAs can competitively bind to miRNAs through their miRNA reaction elements, thus acting as competitive endogenous RNAs (ceRNAs) to regulate miRNAtarget mRNAs' expression levels [20]. More and more studies are revealing the influence of the molecular expression in NPC in the process of IDD disease. However, these results are relatively independent, which cannot provide us a more comprehensive understanding of the molecular expression in the NPC of IDD. Meanwhile, there are still many unknown genes to be explored for a more comprehensive explanation of IDD. There have been few reports on the construction of IDD-related lncRNA/circRNA-miRNA-mRNA networks so far. Therefore, the goal of this study is to filtrate differentially expressed RNAs by using multiple microarray datasets collected from public databases to further construct the regulatory mechanism network of lncRNA/circRNA-miRNA-mRNA in IDD. We aim to identify key molecules in the IDD process and their associated interaction axes, which may provide possible targets for further study.

Materials and Methods
An overall framework of this study, firstly, is shown in Figure 1. The logical framework is mainly divided into "data collection and screening," "data analysis," and "network construction".

Data Collection and Screening.
Public functional genomics data was obtained through Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/). To obtain microarray expression datasets for IDD mRNA, miRNA, lncRNA, and circRNA, we used the search keywords "intervertebral Disc" and "Intervertebral Disc" and limited the organism to Homo sapiens. By reviewing the summaries, overall designs, and sample sources of these datasets manually, we required that the sample source should be NPC. And in the overall designs, we also required that there were no drug therapy or other interventions in the population to ensure that we could use the selected datasets to screen out the differentially expressed molecules of degenerated NPC.

Data
Analysis. Firstly, the differential expression of each dataset was analyzed by using the GEO2R tool [21]. After analysis, the differentially expressed genes (DEGs), differentially expressed miRNAs (DEMs), differentially expressed lncRNAs (DELs), and differentially expressed circRNAs (DECs) between the patients and the control group were obtained. The statistical significance was set to |log fold change ðlogFCÞ | ≥1, adj. p value < 0.05 (raw p value was selected if there were too few molecules screened under the condition of adj. p value < 0.05). Then, the series matrix files of these datasets were downloaded to do heat map analysis by Morpheus (https://software.broadinstitute.org/morpheus). Intersection analysis was accomplished by the Venn drawing tool of Bioinformatics & Evolutionary Genomics (http:// bioinformatics.psb.ugent.be/webtools/Venn/) to gain common differentially expressed RNAs. Besides, functional enrichment analyses on the common DEGs, including Gene  Figure 1: The logical framework is mainly divided into data collection and screening, data analysis, and network construction. The datasets of mRNA, miRNA, lncRNA, and circRNA correlated to IDD were collected from the public database, and DEMs, DEGs, DELs, and DECs with significant expression differences were screened out. Bioinformatics analysis tools and databases were used to conduct regulatory prediction analysis on DEGs, DELs, and DECs and obtain miRNAs paired with them. DEMs was then matched to establish a complete lncRNA/circRNA-miRNA-mRNA network. Then, the network was used to conduct ceRNA analysis and functional enrichment analysis. IDD: intervertebral disc degeneration; DEM: differentially expressed miRNA; DEG: differentially expressed mRNA; DEL: differentially expressed lncRNA; DEC: differentially expressed circRNA; ceRNA: competitive endogenous RNA.

Protein-Protein Interaction (PPI) Network Construction.
In the construction of PPI relationship, we used STRING as the prediction tool. STRING is a known database for predicting protein-protein interactions, including direct physical and indirect functional connections. The interaction relationship is derived from computer prediction, knowledge transfer between organisms, and interactions aggregated from other databases [24]. STRING extract experimental data from numbers of databases like BIND, DIP, and GRID and extract curated data from Biocarta, BioCyc, Reactome, and so on, which guarantee the reliability of mRNA-mRNA relationships in the interaction network. We mapped the DEGs to STRING and validated experimentally only interactions with a combined score > 0:4, which was selected as significant. Hub genes could be found through PPI analysis of common DEGs with the STRING tool, and the PPI network of DEGs was constructed and visualized using Cytoscape software (version 3.6.1; https://www.cytoscape.org.).
2.4. miRNA-mRNA Relationship Construction. miRNA-mRNA interactions were established by TargetScan Human 7.2, which is a tool that predicts biological targets of miRNAs by searching for the presence of conserved 8mer, 7mer, and 6mer sites that match the seed region of each miRNA [25]. The target mRNAs associated with DEMs could be predicted using TargetScan Human 7.2. We entered the gene names of each DEMs at TargetScan Human 7.2 to search for their predictive binding mRNAs, and among the predicted results, we selected those with low total context score and high aggregate PCT which suggest high possibility. Then, the target mRNAs were overlapped with the origin DEGs to find the very DEGs that are paired with DEMs and construct the miRNA-mRNA relationship using Cytoscape software.
2.5. lncRNA-miRNA Relationship Construction. lncRNA-miRNA relationships were constructed by LncBase Predicted v.2, a tool developed by DIANA Lab which presents an extensive collection of miRNA-lncRNA interactions. This significantly enhanced database includes more than 70000 low-and high-throughput, (in) direct miRNA-lncRNA experimentally supported interactions, derived from manually curated publications and the analysis of 153 AGO CLIP-Seq libraries.
LncBase hosts in silico-predicted miRNA targets on lncRNA, identified with the DIANA-microT algorithm [26]. We used LncBase to predict the target lncRNAs related to DEMs, and then, the target lncRNAs were overlapped with the origin DELs to screen the DELs that are paired with DEMs and to construct the lncRNA-miRNA network using Cytoscape software.
2.6. circRNA-miRNA Relationship Construction. We used two tools, circBase and RegRNA 2.0, to help construct the circRNA-miRNA relationship. circBase can explore public circRNA datasets and export FASTA files containing genomic sequences [27], while RegRNA 2.0 is a tool that could comprehensively identify the functional RNA motifs and sites in an input RNA sequence, as well as predicting miR-NAs that are related to the input circRNA [28]. DEC sequences were obtained from circBase, and then, RegRNA 2.0 was used to predict miRNAs that are related to DECs. After these miRNAs were overlapped with DEMs, DECs paired with DEMs was screened to construct the circRNA-miRNA relationship. The constructed lncRNA-miRNA and circRNA-miRNA relationships were then correlated with the miRNA-mRNA network, and Cytoscape software was used to build the lncRNA/circRNA-miRNA-mRNA network.
Then, by reviewing the summaries, overall designs, and sample sources of these datasets manually, we required that the samples of the selected data sets should be NPC and there were no drug therapy or other interventions in the population in the overall designs. It is worth noting that Lan et al. considered that the discs in scoliosis were not strictly normal discs [29], so they reused the normal discs as the control and replaced GSE19943 with GSE63492. Therefore, we also excluded the inclusion of GSE34095. Finally, we selected GSE70362 as the mRNA dataset, GSE63492 and GSE116726 as the miRNA dataset, GSE56081 as the lncRNA dataset, and GSE67566 as the cir-cRNA dataset. According to the preset threshold (adj. p value/raw p value < 0.05 and |logFC | ≥1), we obtained 80 DEGs from GSE70362 (including 32 upregulated and 48 downregulated), 95 DEMs from GSE63492 (including 57 upregulated and 38 downregulated), and 983 DEMs from GSE116726 (including 527 upregulated and 456 downregulated). Besides, 115 DELs were screened out from GSE56081 (containing 50 upregulated and 65 downregulated) and 636 DECs were gained from GSE67566 (consisting of 354 upregulated and 282 downregulated). The top 20 DEGs, DEMs, DELs, and DECs are presented in Tables 2-4. The hierarchical cluster heat maps (Figure 2(a)) indicate that the expression differences of these RNAs were obvious between IDD NPC and normal NPC. Using the Venn drawing tool, the intersection analysis of the two miRNA datasets was performed and 33 DEMs were obtained, as shown in Figure 2(b). The respective expressions of these 33 miRNAs in the two datasets are shown in Figure 2(c), and 14 DEMs with consistent expression were ultimately selected, among which 9 were upregulated and 5 were downregulated. The 80 DEGs were analyzed using DAVID v6.8 to predict their possible biological functions and mechanism pathways. GO enrichment analysis shows 37 terms in total, and Table 5 Figure 3. DEM-binding target genes predicted by TargetScan Human 7.2 were compared with 80 DEGs for overlapping. Thus, 10 DEMs and 49 DEGs formed a miRNA-mRNA interaction relationship and GO enrichment analysis was also performed on these 49 DEGs. In Table 5, these DEGs may be involved in the biological processes of the extracellular region, extracellular exosome, transcriptional activator activity, cell inhibiting, transcription factor binding, etc.
The direct display of the enrichment analysis results is shown in Figure 3. A PPI network consisting of 49 nodes and 21 interaction pairs was built using the STRING tool. Among these DEGs, CYP1B1, NQO1, FOXF2, FOXQ1, and GATA6 form an interaction network together, COL10A1, IBSP, and DLX3 form an interaction axis, and CD1D and CD80 form another interaction axis. Figure 4 shows the protein-protein interaction (PPI) network built by Cytoscape software. Using LncBase Predicted v.2, a total of 5973 lncRNAs were predicted by 10 DEMs, and after eliminating the duplicate data, they were superimposed with 3073 DELs. Thus, the lncRNA-miRNA interaction network between 30 DELs and 8 DEMs was built (11 DELs with downregulated expression regulate 5 DEMs with upregulated expression, while 19 DELs with upregulated expression regulate 3 DEMs).
Through circBase and RegRNA 2.0 database, 505 miR-NAs predicted to interact with DECs were obtained to overlap with 10 DEMs and eventually 4 down-expressed DECs and 3 up-expressed DEMs form an interaction relationship. Table 6 exhibits the lncRNA/circRNA-miRNA-mRNA relationship and Figure 5 shows the interaction network constructed by Cytoscape software.

Discussion
As a ubiquitous health problem that affects global health and socioeconomics, it is critical to further study the biological mechanisms behind IDD. Except for the potential IDD mRNA, noncoding RNAs including miRNA, lncRNA, and circRNA may also be the targets used to develop potential therapy strategies. We have reviewed some existing research on the mechanisms of IDD ceRNA. Wang et al. suggested that lncRNA TRPC7-AS1 regulates nucleus pulposus cellular senescence and ECM synthesis via competing with HPN for miR-4769-5p binding [30], and Yang et al. found that lncRNA SLC20A1 promotes extracellular matrix degradation in NP cells by targeting the miR-31-5p/MMP3 axis [31]. In addition, a research suggested that LINC00969 promotes    [32]. Except for these researches on lncRNA, there are also literatures reporting the mechanism of circRNA, like circu-larRNA_104670 can directly bind to miR-17-3p and correct the negative regulation of miR-17-3p on MMP-2, thus inhibiting the apoptosis of NP cells [33]. And hsa_circ_4099 can act as a "sponge" by competitive binding of miR-616-5p, thus reversing the inhibitory effect of miR-616-5p on Sox9 [34]. Noncoding RNA plays an important role in the process of IDD. However, there is still a lack of a comprehensive understanding of the molecular interactions involved in IDD, as there may be other important functional regulation axes. A broad regulatory network of the disease is in need. Based on the existing data from public databases, our study integrated and analyzed the genes and ncRNAs and formed a molecular interaction network. The network we built consists of 49 DEGs, 10 DEMs, 30 DELs, and 4 DECs, as well as containing 375 potential lncRNA/circRNA-miRNA-mRNA axes. We look forward to finding the genes and ncRNAs that are most likely related to IDD through this net-work, in order to lay a foundation for further research in the future. GO enrichment analysis results show that the metabolism of hyaluronan may be one of the important biological processes for IDD. The hyaluronan synthase activity, hyaluronan biosynthetic process, and hyaluronan catabolic process are shown repeatedly in the GO analysis result list. Hyaluronic acid (HA) is one of the major glycosaminoglycan components of the extracellular matrix and is thought to be involved in cell proliferation, migration, and differentiation. HA plays a crucial role in retaining water of the intervertebral disc, so it can provide flexibility and shock absorbance in the spine [35]. Seen from the GO enrichment analysis results, HYAL1, HAS2, and LYVE1 are common gene molecules in hyaluronan-related metabolism. We searched these three genes in the constructed PPI interaction network and found that they constitute the interaction axis of HYAL1-HAS2-LYVE1. These findings suggest that they may indeed participate in the metabolism of hyaluronan and play a common regulatory role in IDD. By reviewing literatures, we found that HAS2 is a member of the gene family that encodes HA  9 BioMed Research International synthase, which was elevated in a study of HA-based hydrogel-induced NPC amplification [36]. HYAL1 encodes lysosomal hyaluronidase, which degrades HA in cells. In one study, the gene expression of HYAL1 and the protein expression of HAYL2 were significantly increased in moderate/severe IDD samples compared with those without or with low IDD [37]. The protein encoded by LYVE1 acts as a receptor and binds to both soluble and fixed HA and may also play a role in lymphatic HA transport. It is also been reported in the literature that ingrowth of vascularized fibrous tissue was seen at the edge or within fragments of degenerate disc tissue. Some scattered small vessels were lined by LYVE1+ endothelial cell [38]. In Reactome (https://www.reactome.org), we found a pathway through which HAS2, HYAL1, and LYVE1 coparticipate in HA biological processes ( Figure 6). This pathway indicates that the integral membrane dual-action glycosyltransferase proteins hyaluronan synthases 1-3 (HAS1-3) mediate the polymerization of glucuronic acid (GlcA) with N-acetylglucosamine (GlcNAc) to form HA. The resulting polymer has the arrangement [-4GlcA-1,3GlcNAc-]n [39][40][41]. LYVE1, HMMR, STAB2, and CD44 together form the HAR complex. High molecular weight HA is tethered to the cell surface by HA receptors and the GPI-linked hyaluronidase 2 (HYAL2) to form a HA : HAR : HYAL2 complex in the plasma membrane that localizes to caveolae, continuing to complete the degradation of HA [42][43][44][45]. In the acidic environment of the lysosome, hyaluronidase 1 (HYAL1) is able to hydrolyze large 50 disaccharide unit HA fragments to 2 disaccharide units [46]. These theories suggest that the prediction of interaction relationship in this study has certain effectiveness.
In the interaction network constructed in this study, HAS2, HYAL1, and LYVE1 may all be negatively regulated by miR-4741. Meanwhile, mi663ahg/HEIH (lncRNA) and hsa-circ-0003600 (circRNA) may competitively be bound on miR-4741 to act as ceRNAs, correcting the negative regulation of miR-4741. Through literature searching, we found that the mechanism of miR663AHG/HEIH/hsa-circ-0003600 is still lacking in the study of IDD. The interaction axis of miR663AHG/HEIH/hsa-circ-0003600-miR-4741-HAS2/HYAL1/LYVE1 may be one of the potential mechanisms to be studied in IDD.
We also noticed that COL10A1 was the gene with the greatest difference in expression between degenerative NP and normal NP (  Figure 4: Protein-protein interaction network. Network nodes represent the proteins produced by a single protein-coding gene locus. Blue represents the proteins that jointly contribute to a shared function but do not necessarily mean that they are physically binding each other. Grey represents the proteins without association. Triangle means upregulated expression, whereas arrow means downregulated expression. associated with some existing research conclusions. COL10A1 encodes type X collagen, which is specifically expressed and synthesized by hypertrophic chondrocytes, and it plays an important role in the process of cartilage ossification and may be related to matrix degradation, calcification, and vascular invasion [47]. Type X collagen can be localized in combination with advanced intervertebral disc degeneration, and the positive staining of type X collagen is more obvious in the intervertebral disc matrix at the late stage of IDD in elderly patients [48]. In some experiments on mice that have led to disc degeneration, abnormal expression of type X collagen, increased apoptosis, inappropriate intervertebral disc ossification, and elevated expression of COL10A1, Runx2, MMP-13, etc. were observed in degenerative NP [11,[49][50][51]. Our network analysis also has confirmed that COL10A plays an important role in the pathogenesis of  Figure 5: Competing endogenous RNA interaction network of lncRNA/circRNA-miRNA-mRNA. Yellow represents lncRNA, green represents circRNA, red represents miRNA, and blue represents mRNA. Triangle means upregulated expression, whereas arrow means downregulated expression.

11
BioMed Research International IDD. This not only verifies the credibility of our analysis but also provides an additional theoretical basis for further study on COL10A.
In the lncRNA/circRNA-miRNA-mRNA interaction network constructed in this study, COL10A1 may be negatively regulated by miR-5100, while lncRNA LINC00475, SOX2-OT, NEAT1, SNHG12, GAS5, AGAP2-AS1, and TARID may act as ceRNA for competitively binding miR-5100 and thus correcting the negative regulation of miR-5100. Wang et al. verified the contribution of miR-5100 to osteoblast differentiation and mineralization through functional gain and loss experiments, but there is still a lack of studies on miR-5100 in hypertrophic chondrocytes, type X collagen, or IDD. For lncRNA that may negatively regulate miR-5100, the current studies of LINC00475, Sox2-OT, SNHG12, AGAP2-AS1, and TARID are focused on various types of tumours, such as glioma, lymphoma, ovarian cancer, and renal clear cell sarcoma [52][53][54][55], while Wang et al. sug-gested that the overexpression of lncRNA GAS5 might promote the apoptosis of NPC through downregulation of Bcl-2 and upregulation of caspase-3 [56]. It is worth noting that NEAT1 may be the ceRNA-regulating COL10A1. Ruan et al. suggest that NEAT1 has been demonstrated to participate in ECM remodelling and they found the that increased expression of NEAT1 decreased aggrecan and type II collagen levels and increased ADAMTS4 and MMP-13 levels [57]. Our analysis results showed that NEAT1 expression was increased in IDD patients (logFC = 1:78), further confirming that NEAT1 may play an important role in IDD. Some research has proved that a robust expression of type X collagen and MMP-13 mRNA and protein could occur by hypertrophic chondrocytes during IDD [51,58]. This suggests that NEAT1 and COL10A1 are connected by type X collagen and MMP-13. The ceRNA network shows that both NEAT1 and COL10A1 can interact with miR-5100 to form the interaction axis. Thus, we build a schematic diagram of  Figure 6: Schematic diagram of hyaluronic acid-related IDD molecular mechanism. The integral membrane dual-action glycosyltransferase proteins hyaluronan synthases 1-3 (HAS1-3) mediate the polymerization of glucuronic acid (GlcA) with N-acetylglucosamine (GlcNAc) to form HA. The resulting polymer has the arrangement [-4GlcA-1,3GlcNAc-]n. LYVE1, HMMR, STAB2, and CD44 together form the HAR complex. High molecular weight HA is tethered to the cell surface by HA receptors and the GPI-linked hyaluronidase 2 (HYAL2) to form a HA : HAR : HYAL2 complex in the plasma membrane that localizes to caveolae, continuing to complete the degradation of HA. In the acidic environment of the lysosome, hyaluronidase 1 (HYAL1) could hydrolyze large 50 disaccharide unit HA fragments to 2 disaccharide units. HAS2, HYAL1, and LYVE1 may be negatively regulated by miR-4741, while mi663ahg/HEIH (lncRNA) and hsa-circ-0003600 (circRNA) may competitively be bound on miR-4741 to act as ceRNAs, correcting the negative regulation of miR-4741 to prevent IDD. 12 BioMed Research International the IDD molecular mechanism related to NEAT1 and COL10A1 ( Figure 7). Although it can be hypothesized that in IDD patients, the NEAT1-miR-5100-COL10A1 interaction axis may be unbalanced, leading to the occurrence of IDD; there is still a lack of evidence that NEAT1 can negatively regulate miR-5100 and thereby regulate COL10A1 to realize ECM remodelling. Further elucidation is needed to confirm this hypothesis in the future.

Conclusions
Through bioinformatics analysis of multiple datasets, we constructed a ceRNA interaction network of lncRNA/-circRNA-miRNA-mRNA for IDD. The result of this study verifies that the process of cartilage ossification, matrix degradation, and calcification following the synthesis of type X collagen and the imbalance of hyaluronic acid metabolism are two explanations for the occurrence of IDD. NEAT1-miR-5100-COL10A1 and miR663AHG/HEIH/hsa-circ-0003600-miR-4741-HAS2/HYAL1/LYVE1 might be potential interaction axes of the molecular mechanism in IDD. The increased expression of NEAT1 might inhibit miR-5100 and subsequently upregulate the expression of COL10A1, which leads to IDD, while the increased expression of miR663AH-G/HEIH/hsa-circ-0003600 might inhibit miR-4741 and indirectly upregulate HAS2/HYAL1/LYVE1, and leads to the protection from IDD. This novel ceRNA network could be used to predict potential mechanisms to reach a comprehensive understanding and could provide possible targets for the further study of IDD in the future.

Data Availability
All data generated or analyzed during this study are included in this published article. The datasets analyzed during the current study are available in the Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/).

Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this paper.

Authors' Contributions
Junshen Huang and Yuxi Li contributed equally to this work.