This article is mentioned in:

Introduction

Intervertebral disc degeneration (IDD) is among the causes of lower back pain, which leads to disability and societal burden, and has a negative socioeconomic impact worldwide (1). It is estimated that worldwide, lower back pain affects ≤84% of people at a certain point in their lifetime, with 10% becoming chronically disabled (2). Currently, therapies for IDD primarily aim to relieve pain symptoms, offering only temporary benefits rather than a permanent cure (3). Therefore, it is necessary to elucidate the underlying mechanisms of IDD and lower back pain in order to develop novel effective treatment methods for these conditions.

The intervertebral disc is a complex structure that consists of three anatomical sub-structures, including the nucleus pulposus (NP), annulus fibrosus and cartilaginous endplates, which are essential for the maintenance of normal function of the disc (4). Early IDD typically occurs in the central NP and is characterized by a decrease in the number of active cells and a loss of extracellular matrix within the NP region (5). Previous studies have identified certain molecules associated with IDD (6,7). For instance, one study reported that early IDD is associated with the downregulation of caveolin-1 expression and inhibition of canonical Wnt signaling (6). Furthermore, a recent study indicated that long non-coding RNAs (lncRNAs), a subset of non-protein-coding transcripts that are >200 nucleotides in length, are involved in a series of biological processes (8). Additionally, competing endogenous RNAs (ceRNAs) participate in microRNA (miRNA)-dependent crosstalk and contribute to the complexity of gene regulation, and may aid in the understanding of disease pathobiology (9). lncRNAs regulate the expression of neighboring coding genes through various mechanisms and may function as ceRNAs to modulate the derepression of miRNA targets (10). In recent years, the roles of lncRNAs in several diseases and cancers have been investigated (11,12). However, potential roles of lncRNAs in IDD are yet to be elucidated.

Wan et al (13) performed microarray data profiling of NP samples from patients with IDD and identified numerous differentially expressed lncRNAs and mRNAs in degenerative discs. Reanalysis of previous datasets using different, updated methods may lead to the identification of novel molecular markers associated with the analyzed dataset (14,15). In the present study, lncRNA and mRNA microarray data from patients with IDD, including 5 samples with degenerative lumbar NP and 5 normal controls, were downloaded and reanalyzed. Differentially expressed lncRNAs or differentially expressed genes (DEGs) were identified and co-expression network analysis was performed, followed by functional analysis of genes in the network. Additionally, a ceRNA regulatory network was constructed based on DEGs and lncRNAs in the co-expression network. Furthermore, a literature search was performed to identify specific miRNAs associated with IDD and a specific miRNA-associated ceRNA network was extracted from the co-expression network. The present study aimed to identify novel IDD-associated lncRNAs and genes, which may guide further investigations and contribute to the development of novel therapeutic strategies for the treatment of IDD.

Materials and methods

Collection of microarray data

The lncRNA and mRNA microarray of IDD was downloaded from the Gene Expression Omnibus database (16) (http://www.ncbi.nlm.nih.gov/geo/) under the accession number GSE56081 (13), which is a subseries of the SuperSeries GSE67567 (13,17,18). The Arraystar Human lncRNA microarray (version 2.0; Agilent_033010 Probe Name version; Arraystar Inc., Rockville, MD, USA) platform was used. A total of 5 degenerative lumbar NP samples derived from patients with IDD and 5 non-degenerative specimens from cadaveric donors were included in this dataset. As described in the original study by Wan et al (13), all specimens were collected within 3 h following disc excision (for discs from patients with IDD) or mortality (for discs from cadaveric donors).

Data preprocessing and DEG screening

Raw data was initially normalized using the robust multi-array average method (19). Subsequently, limma package version 3.28.21 (20) of Bioconductor 3.5 (http://www.bioconductor.org/packages/3.5/bioc/html/limma.html) was adopted to identify DEGs between the IDD samples and normal controls. False discovery rate (FDR) using Benjamini-Hochberg method (21) in multi-test package version 2.2.0 of Bioconductor (22) (http://www.bioconductor.org/packages/2.5/bioc/html/multtest.html) was used to correct multiple testing to avoid false positives. Only genes and lncRNAs with |log2fold change (FC)| >1 and FDR <0.05 were selected as significantly differentially expressed. Sequences and probe IDs were provided by the platform but not gene and lncRNA names or symbols, therefore, the sequences were aligned with hg19 genome annotation (Ensembl) (23) to identify genes corresponding with each ID. Furthermore, heatmaps were constructed to present expression profiles of differentially expressed lncRNAs and DEGs using hierarchical clustering, which was performed using pheatmap package (http://cran.r-project.org/web/packages/pheatmap/index.html) in R based on Euclidean distances.

Co-expression network analysis of differentially expressed lncRNAs and DEGs

To determine the functions of differentially expressed lncRNAs, Pearson correlation analysis of the DEGs and differentially expressed lncRNAs was performed using the cor() function in R. Pearson correlation coefficients between-1 and 0 or 0 and 1 were considered to represent positive and negative correlation, respectively. In the present study, co-expressed lncRNA-mRNA pairs were identified when absolute Pearson correlation coefficient was >0.6. In addition, lncRNA-mRNA co-expression network was visualized using Cytoscape version 2.8.0 (24).

Functional analysis of genes in the network

An ontology-based tool, clusterProfiler, offers three methods for biological term classification and enrichment analyses of genes, including group Gene Ontology (GO), enrichGO and enrich Kyoto Encyclopedia of Genes and Genomes (KEGG) (25). In the present study, in order to elucidate the biological function of genes enriched in the co-expression network, GO biological process (BP) terms and KEGG pathway analysis of genes in the network were performed using enrichGO and enrichKEGG functions in R package clusterProfiler with P<0.05 as a cutoff value (25).

Construction of ceRNA regulatory network

To elucidate the roles of lncRNAs and mRNAs in the co-expression network, data pertaining to lncRNAs and mRNAs was combined with miRNA data in order to construct a miRNA-lncRNA-mRNA ceRNA regulatory network. Initially, miRNAs that had regulatory associations with lncRNAs in the co-expression network were searched for in two databases; miRcode (26) (http://www.mircode.org/) and starBase (27) (http://starbase.sysu.edu.cn/). Subsequently, the targets of the identified miRNAs were determined using the miRTarBase database (http://mirtarbase.mbc.nctu.edu.tw), which includes experimentally validated microRNA-target interactions (28). The overlapping genes between targets of the identified miRNAs and mRNAs contained in the co-expression network were retained for ceRNA construction. ceRNA and disease-specific ceRNA networks were subsequently constructed. Based on the identified miRNAs, miRNA targets in the co-expression network and differentially expressed lncRNAs, a ceRNA regulatory network was constructed and visualized using Cytoscape version 2.8.0 (24). Furthermore, specific miRNAs associated with IDD were identified in previously reported articles by combining search terms ‘intervertebral disc degeneration’ AND ‘human’ AND ‘miRNA’ in PubMed (https://www.ncbi.nlm.nih.gov/pubmed/), and a IDD-specific ceRNA network containing specific miRNAs was extracted from the co-expression ceRNA-regulated network and visualized using Cytoscape version 2.8.0 (24).

Results

Identification of DEGs and differentially expressed lncRNAs

Following data normalization, differential expression analysis of genes and lncRNAs was performed using the dataset that included 5 IDD samples and 5 controls. A total of 967 genes and 137 lncRNAs were significantly differentially expressed within the FDR <0.05 and |log2FC| >1 criteria. Among them, 737 genes were upregulated and 230 genes were downregulated; 60 lncRNAs were upregulated and 77 lncRNAs were downregulated. Additionally, hierarchical clustering of gene expression levels and lncRNA expression levels indicated clustering of normal groups and IDD groups (Fig. 1).

Figure 1. Heatmaps of DEGs and differentially expressed lncRNAs. (A) Heatmap of DEGs between IDD and non-degenerative control samples. (B) Heatmap of differentially expressed lncRNAs between IDD and non-degenerative samples. Red color indicates upregulation and green color indicates downregulation. Normal and degenerated refer to samples obtained from non-degenerative controls and patients with IDD, respectively. DEGs, differentially expressed genes; lncRNA, long noncoding RNA; IDD, intervertebral disc degeneration.

Differential lncRNA-mRNA co-expression network analysis

Differentially expressed lncRNAs and DEGs with Pearson correlation coefficients >0.6 were selected. The lncRNA-mRNA co-expression network, containing 248 nodes and 1,082 edges (448 negative and 634 positive connections), was constructed (Fig. 2). A total of 39 differentially expressed lncRNAs (27 downregulated and 12 upregulated) and 209 DEGs (75 downregulated and 134 upregulated) were included in this co-expression network.

Figure 2. lncRNA-mRNA co-expression network. Triangular and vee-shaped nodes represent upregulated and downregulated genes, respectively. Diamond-shaped and rectangular nodes indicate upregulated and downregulated lncRNAs, respectively. lncRNA, long noncoding RNA. Continuous line represents the connections between lncRNA and miRNA; dotted lines represent connections between miRNA and gene or lncRNA and gene.

Functional analysis of genes and lncRNAs in the co-expression network

To determine the potential functions of the identified DEGs and lncRNAs in the co-expression network, GO and KEGG pathway enrichment analyses were performed in the co-expression network. The analyses revealed the significant enrichment of 17 GO BP terms and 10 KEGG pathways (Table I; Fig. 3). BP and KEGG terms that were significantly associated with members of the co-expression network included ‘skeletal system development’, ‘response to mechanical stimulus’, ‘bone development’ and ‘focal adhesion’ (Table I; Fig. 3).

Figure 3. Functional analysis of members of the long noncoding RNA-mRNA co-expression network. Bar plots presenting enriched (A) GO biological process terms and (B) Kyoto Encyclopedia of Genes and Genomes pathways in the co-expression network. The black line indicates-log10 (P-value). GO, Gene Ontology.

Table I. GO and KEGG pathway enrichment analysis for genes in the lncRNA-mRNA co-expression network.

Table I.

GO and KEGG pathway enrichment analysis for genes in the lncRNA-mRNA co-expression network.

A, GO terms enriched in the lncRNA-mRNA co-expression network
Description (ID, name)	Count	FDR	Genes
GO:0001501, skeletal system development	14	0.000067700	RBP4, AEBP1, LGALS3, TUFT1, KLF10, PRELP, MEN1, MEF2D, SBDS, LECT1, TNFSF11, COL1A2, ACAN, PKD1
GO:0009612, response to mechanical stimulus	5	0.003660942	FOS, BTG2, JUN, PKD1, GFI1
GO:0007179, transforming growth factor beta receptor signaling pathway	5	0.004694322	FOS, LTBP2, KLF10, JUN, COL1A2
GO:0040008, regulation of growth	11	0.005411270	MEN1, EPC1, RBP4, HTRA1, LGMN, BCAR1, NOL8, SIPA1, PEX5, GREM1, GAS6
GO:0001503, ossification	6	0.009830659	MEN1, MEF2D, SBDS, TNFSF11, TUFT1, KLF10
GO:0060348, bone development	6	0.012887071	MEN1, MEF2D, SBDS, TNFSF11, TUFT1, KLF10
GO:0006357, regulation of transcription from RNA polymerase II promoter	16	0.017700242	JDP2, EGR2, KLF9, KLF10, FOXO1, NFYC, MEN1, PLAGL1, FOS, EPC1, SMARCD3, JUN, GFI1, LCOR, KLF4, NR1H3
GO:0045893, positive regulation of transcription, DNA-dependent	12	0.019627483	PLAGL1, MEN1, FOS, EPC1, EGR2, SMARCD3, JUN, FOXF1, NFYC, FOXO1, KLF4, NR1H3
GO:0051254, positive regulation of RNA metabolic process	12	0.020731482	PLAGL1, MEN1, FOS, EPC1, EGR2, SMARCD3, JUN, FOXF1, NFYC, FOXO1, KLF4, NR1H3
GO:0009628, response to abiotic stimulus	10	0.024048392	MEN1, FOS, RBP4, TNFSF11, BTG2, HMGCR, JUN, PKD1, GFI1, GRK1
GO:0007178, transmembrane receptor protein serine/threonine kinase signaling pathway	5	0.029381231	FOS, LTBP2, KLF10, JUN, COL1A2
GO:0030198, extracellular matrix organization	5	0.030296507	LGALS3, LUM, FOXF1, COL1A2, ACAN
GO:0007611, learning or memory	5	0.037173540	FOS, SBDS, EGR2, HMGCR, JUN
GO:0010605, negative regulation of macromolecule metabolic process	15	0.038521124	JDP2, KLF10, CDC26, SIGIRR, TIMP1, MEN1, EPC1, SET, JUN, SERPINA5, BNIP3L, GFI1, LCOR, KLF4, ADAR
GO:0051094, positive regulation of developmental process	8	0.038680685	MEN1, TNFSF11, KLF10, JUN, CLU, GFI1, OTP, ANGPTL4
GO:0042592, homeostatic process	15	0.045211148	CLNS1A, RBP4, EGR2, NCDN, MBP, TIMP1, MEN1, TNFSF11, PRDX6, JUN, SERPINA3, TXNRD1, XCR1, QSOX1, CLN5
GO:0045597, positive regulation of cell differentiation	7	0.045820821	MEN1, TNFSF11, KLF10, JUN, CLU, GFI1, OTP
B, KEGG pathways enriched in the lncRNA-mRNA co-expression network
Description (entry, name)	Count	FDR	Genes
hsa04142, Lysosome	4	0.014272775	CLTA, LGMN, CLN5, CTSF
hsa04670, Leukocyte transendothelial migration	4	0.014534463	BCAR1, SIPA1, CD99, ITGAM
hsa03320, PPAR signaling pathway	3	0.018064527	FADS2, ANGPTL4, NR1H3
hsa04510, Focal adhesion	5	0.018720786	PARVG, PPP1CA, JUN, BCAR1, COL1A2
hsa04640, Hematopoietic cell lineage	3	0.025076738	CD1B, CSF3R, ITGAM
hsa03020, RNA polymerase	2	0.027129221	POLR2C, POLR2J2
hsa00240, Pyrimidine metabolism	3	0.028840787	TXNRD1, POLR2C, POLR2J2
hsa00051, Fructose and mannose metabolism	2	0.031922905	PFKFB2, AKR1B1
hsa04810, Regulation of actin cytoskeleton	4	0.043576076	PPP1CA, BCAR1, FGF13, ITGAM
hsa04514, Cell adhesion molecules (CAMs)	3	0.043833628	CD99, SELPLG, ITGAM

[i] GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; lncRNA, long noncoding RNA; count, number of genes enriched in the GO terms or pathways; FDR, false discovery rate.

Construction of ceRNA regulatory network and disease-specific ceRNA network

A total of 79 associated miRNAs, which were targeted by 9 lncRNAs in the co-expression network, were screened by integrating information from miRcode and starBase databases to determine regulatory associations between miRNAs and lncRNAs. A total of 647 pairs of miRNA-gene interactions were obtained by identifying targets of these miRNAs from genes in the co-expression network using the miRTarBase database. The ceRNA regulatory network was therefore constructed based on combining the identified miRNA-gene and lncRNA-miRNA associations (Fig. 4). The ceRNA network was composed of 236 nodes [79 miRNAs, 9 downregulated lncRNAs, 148 DEGs (including 46 downregulated and 102 upregulated DEGs)] and 1,065 edges (647 miRNA-gene connections, 142 lncRNA-gene connections and 276 lncRNA-miRNAs connections).

Figure 4. Construction of ceRNA network. Triangular and vee-shaped nodes denote upregulated and downregulated genes, respectively. Diamond-shaped nodes indicate downregulated long noncoding RNAs and rectangular nodes indicate miRNAs. ceRNA, competing endogenous RNA; miRNA/miR, microRNA.

Furthermore, IDD-associated miRNAs from previous studies were searched and the collected data indicated that hsa-miR-140 may be implicated in human IDD (29,30). Subsequently, connections between hsa-miR-140 and other DEGs, lncRNAs and miRNAs were extracted from the co-expression ceRNA-regulated network and a hsa-miR-140-associated ceRNA network was constructed (Fig. 5). The disease-specific ceRNA network included 47 nodes (1 miRNA, 3 lncRNAs, 27 downregulated DEGs, and 16 upregulated DEGs) and 55 edges. hsa-miR-140 was demonstrated to be regulated by three lncRNAs in the hsa-miR-140-associated ceRNA network, including KCNQ1 opposite strand/antisense transcript 1 (KCNQ1OT1), OIP5 antisense RNA 1 (OIP5-AS1) and UGDH antisense RNA 1 (UGDH-AS1). The results of the present study indicated that the aforementioned lncRNAs may be associated with IDD and the genes that are co-expressed with these lncRNAs may also be associated with IDD. KCNQ1OT1 was co-expressed with neurochondrin (NCDN) and lon peptidase 2, peroxisomal (LONP2). lncRNAs OIP5-AS1 and UGDH-AS1 were co-expressed with several overlapping genes, including forkhead box F1 (FOXF1) and polycystin 1, transient receptor potential channel interacting (PKD1).

Figure 5. Disease-specific competing endogenous RNA network. Pink and green circular nodes denote upregulated and downregulated genes, respectively. Green diamond-shaped nodes denote downregulated lncRNAs and red square nodes indicate intervertebral disc degeneration-associated miRNAs. Blue lines represent co-expression associations between lncRNAs and genes. Red lines indicate lncRNA-miRNA regulation associations. lncRNA, long noncoding RNA; miRNA, microRNA.

Discussion

In the present study, 967 genes (737 upregulated and 230 downregulated) and 137 lncRNAs (60 upregulated and 77 downregulated) were differentially expressed between the IDD samples and controls. A co-expression network was constructed containing 39 differentially expressed lncRNAs and 209 DEGs that were primarily involved in ‘skeletal system development’, ‘response to mechanical stimulus’, ‘bone development’ and ‘focal adhesion’. Furthermore, a ceRNA network was established, which included 79 miRNAs, 9 downregulated lncRNAs and 148 DEGs. Among the miRNAs, hsa-miR-140 was a previously reported disease-associated miRNA and therefore, a hsa-miR-140-associated ceRNA network was extracted. It was determined that hsa-miR-140 was regulated by three lncRNAs, including KCNQ1OT1, OIP5-AS1 and UGDH-AS1. KCNQ1OT1 was co-expressed with NCDN and LONP2. In addition, lncRNAs OIP5-AS1 and UGDH-AS1 targeted several overlapping co-expressed genes, including FOXF1 and PKD1.

Aberrant expression of miR-140 was reported to contribute to osteoarthritis, a degenerative disease characterized by pathological alterations similar to those identified in patients with IDD (31,32). miRNA-140 has also been demonstrated to serve roles in cartilage development and homeostasis by regulating the expression of matrix metalloprotease-13, insulin-like growth factor binding protein-5, Sp1 transcription factor, interleukin-1 and Smad family member 3 (31,33,34). The results of the present study support the hypothesis that miR-140 serves a role in the pathogenesis of degenerative disorders such as IDD. The present study identified three lncRNAs, KCNQ1OT1, OIP5-AS1 and UGDH-AS1, which were associated with disease-specific miRNAs, including miR-140.

KCNQ1OT1 has been demonstrated to mediate transcriptional silencing (35). In the miR-140-associated ceRNA network, KCNQ1OT1 was co-expressed with NCDN, which was enriched in the GO term ‘positive regulation of developmental process’. A recent study reported that NCDN served as a neuronal target antigen in patients with autoimmune cerebellar degeneration (36). NCDN, a cytoplasmic leucine-rich protein that functions in neurite outgrowth and chondrocyte differentiation, was demonstrated to be NP-specific and age-associated, indicating its potential functional implication in the development of NP (37). Therefore, we hypothesize that the downregulation of KCNQ1OT1 may be associated with the degeneration of the intervertebral disc through the regulation of NCDN expression.

Previous studies have demonstrated that lncRNA OIP5-AS1 acts as a sponge to prevent RNA-binding protein HuR from binding to its target mRNA and reduce cell proliferation (38). Stimulation of disc cell proliferation and accumulation of extracellular matrix may slow or even reverse the degenerative alterations resulting from IDD (7). To the best of our knowledge, the role of UGDH-AS1 has not been previously investigated. In the present study, OIP5-AS1 and UGDH-AS1 shared several overlapping co-expressed genes, including FOXF1 and PKD1. FOXF1 is a member of the forkhead family of transcription factors, which are characterized by a distinct forkhead domain (39). A previous study indicated that FOXF1 promoted the production of extracellular matrix (40). Importantly, IDD is associated with an increased breakdown of extracellular matrix (41). Boulter et al (42) reported that PKD1 may serve a role in skeletal development. Based on the results of the aforementioned studies and data analyzed in the present study, it may be hypothesized that the lncRNAs OIP5-AS1 and UGDH-AS1 are involved in the molecular mechanisms underlying IDD by regulating the expression of FOXF1 and PKD1. The results of the present study require further investigation to clarify and confirm the roles of the identified molecules in IDD.

The present study had several limitations, including the small sample size of the dataset employed; further analyses with an increased sample size are necessary to confirm the results. Furthermore, non-degenerative specimens from cadaveric donors were used as controls in the dataset that was analyzed in the present study. Patients with normal discs or with mild disc degeneration are not typically subjected to surgical intervention. The use of normal discs from cadaveric donors reflected limited access to less degenerated tissues from surgical patients. However, cadaveric tissues were collected within 3 h following mortality to ensure that the tissues remained intact (13). Additionally, the results of the analysis in the present study were not experimentally verified. IDD is a complex and multi-factorial process that is influenced by genetic predispositions, aging, biomechanical loading and physical activities, lifestyle and other health-associated factors, including diabetes (43). Numerous clinical investigations have indicated that cell-based biological therapies, including disc cell augmentation, application of growth factors/cytokines and gene therapy, may be employed for the treatment of IDD (44,45). The development of biomarkers for IDD is valuable in the early stages of disc degeneration (45). lncRNAs have been implicated in multiple pathological processes of IDD (46), and lncRNAs and DEGs in patients with disc degeneration may have potential as diagnostic or prognostic biomarkers (47). In order to determine the potential for the clinical application of lncRNAs and DEGs identified in the present study, an in vivo murine IDD model may be used to confirm the roles of these molecules in IDD in future studies.

In conclusion, the present study identified several potential genes and lncRNAs, including KCNQ1OT1, OIP5-AS1, UGDH-AS1, NCDN, FOXF1 and PKD1, that were associated with IDD. KCNQ1OT1 may be involved in the degeneration of the intervertebral disc by regulating the expression of NCDN, while the lncRNAs OIP5-AS1 and UGDH-AS1 may be implicated in the molecular mechanisms of IDD by affecting the expression of FOXF1 and PKD1. The results of the present study may contribute to elucidating the molecular pathogenesis of IDD and aid in the identification of novel biomarkers associated with IDD.

Acknowledgements

Not applicable.

Funding

No funding was received.

Availability of data and materials

The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.

Authors' contributions

ZGQ, ZQ and YL provided the concept and design of the research. QZ, ZYW and QS acquired the data. XZ, CF and ZW analyzed and interpreted the data. FX, YDL and YW performed the statistical analysis. ZGQ and ZQ drafted the manuscript. YL and ZW revised the manuscript.

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

References