Gene Expression Omnibus (GEO) (www.ncbi.nlm.nih.gov/geo/) is an international public database that archives high-throughput sequencing gene expression data. In order to explore the different expression profiling in IDD, we searched datasets from the GEO database with the keywords “intervertebral disc degeneration” AND “Homo sapiens” AND “gse”. Expression data generated from the human annulus disc tissue of IDD patients and healthy controls were incorporated into our work. Four mRNA expression datasets (including GSE23130, GSE15227, GSE17077, and GSE70362) and 3 miRNA expression datasets (GSE116726, GSE63492, and GSE45856) were included. The basic information of datasets is shown in Table 1. This study has been approved by the Ethics Committee of Tianjin First Central Hospital.

In order to minimize the heterogeneity among different datasets, raw data were performed for log2 transformation and normalization, and then, the metaMA package was used to combine data from multiple datasets. Individual p values were calculated, and false discovery rate (FDR) was obtained by using the Benjamini–Hochberg method. Differentially expressed mRNAs (DEMs), miRNAs (DEMis), circRNAs (DECs), and lncRNAs (DELs) were investigated in IDD. In our work, mRNAs with p value < 0.05, miRNAs with FDR<0.01, circRNAs with FDR<0.01, and lncRNAs with FDR<0.01 were considered as DEMs, DEMis, DECs, and DELs, respectively.

The DEM-associated target genes were predicted using the miRwalk3 (http://mirwalk.umm.uni-heidelberg.de/), which stores predicted data obtained with a machine learning algorithm including experimentally verified miRNA–target interactions. The target genes were then overlapped with the DEMs, and the negative interaction pairs between DEMis and DEMs (according to their expression levels) were used to construct the DEMi–DEM network using Cytoscape software (version 3.6.1; www.cytoscape.org).

The starBase database (version 2.0; starbase.sysu.edu.cn/index.php) (Cardenas-Gonzalez et al., 2017) was used to screen the interactions between DELs and DEMis, which were then integrated with the miRNA–mRNA interactions to establish the DEL–DEMi–DEM ceRNA network using Cytoscape software (version 3.6.1; www.cytoscape.org). Human sequences of DECs and DEMis were downloaded from the circBase (www.circbase.org) (Li et al., 2020) and miRBase (version 21; www.mirbase.org) (Chen et al., 2019) databases, respectively. miRanda (cbio.mskcc.org/miRNA2003/miranda.html) (Ren et al., 2020) was used to predict the interactions between DECs and DEMis. The interaction pairs between DECs and DEMs were then integrated with the DEMi–DEM interactions to establish the DEC–DEMi–DEM ceRNA network using Cytoscape software (version 3.6.1; www.cytoscape.org). The overlapped DEMi–DEM in the aforementioned two ceRNA networks was also selected to construct the lncRNA/circRNA–miRNA–mRNA network.

In order to explore the diagnostic value of DEMis in IDD, the pROC package in R language curves was used to depict the receiver operating characteristic (ROC), and the area under the curve (AUC) was calculated. The DEMis with AUC≥0.8 were considered as having the performance of distinguishing IDD patients from healthy controls. The diagnostic value of DEMis in GSE116726, GSE63492, and GSE45856 datasets was investigated in the present study.

The biological functions of DEMs in IDD were predicted by both Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) function and pathway through online software Genecodis3 (http://genecodis.cnb.csic.es) as a non-redundant and modular enrichment analysis tool for functional genomics (Tabas-Madrid et al., 2012). The enriched KEGG pathway with FDR <0.05 was the significant enrichment term.

The expression levels of DEC, DEL, and DEM candidates were explored in external datasets, including GSE124272, GSE150408, and GSE153761. Both GSE124272 (8 IDD patients and 8 healthy individuals) and GSE150408 (17 IDD patients and 17 healthy individuals) store transcriptomic profiling data of the whole blood of patients with IDD and healthy individuals. GSE153761 (3 IDD patients and 3 healthy individuals) stores transcriptomic profiling data of the cervical cartilage endplate of patients with IDD and healthy individuals.

The statistical significance between groups was assessed by unpaired Student’s t-test. p < 0.05 was statistically significant. * indicates p < 0.05; ** indicates p < 0.01, and *** indicates p < 0.001.

A total of four mRNA expression datasets including 52 healthy controls (HCs) and 31 IDD patients were used to identify the DEMs in IDD. A total of 1995 DEMs including 907 down-regulated and 1,078 up-regulated DEMs were identified in IDD vs. HC. Three miRNA expression datasets including 11 HCs and 11 IDD patients were used to identify the DEMis in IDD. Eighty-four DEMis including 10 up-regulated DEMis and 74 down-regulated DEMis were identified in IDD compared with HCs. A total of 256 DELs including 107 up-regulated and 149 down-regulated DELs were identified in IDD patients. In addition, 589 DECs including 328 up-regulated and 261 down-regulated DECs were identified in IDD patients. The hierarchical cluster heatmaps indicated that these DECs (Figure 1A), DELs (Figure 1B), DEMis (Figure 1C), and DEMs (Figure 1D) could distinguish IDD from control samples.

We constructed the interaction network between miRNAs and target genes in IDD based on the identified miRNA–target gene interaction pairs of reverse association using Cytoscape software. As shown Figure 2, 8,673 miRNA–target gene pairs of reverse correlation between 80 DEMis and 1,268 DEMs were identified in IDD, which included 7,906 pairs between 70 down-regulated DEMis and 885 up-regulated DEMs (Figure 2A) and 767 pairs between 10 up-regulated DEMis and 383 down-regulated DEMs (Figure 2B). In up-regulated DEMis pairs, miR-3189-3p, miR-3714, miR-1291, and miR-302c-5p had the high connectivity with DEMs, which interacted with 126, 112, 97, and 93 DEMs. In the down-regulated DEMis pairs, miR-4270, miR-3162-5p, miR-320b, miR-3138, miR-3198, miR-3679-5p, miR-3622a-5p, and miR-4257 interacted with more than 200 DEMs, respectively.

The underlying functions of the DEMs in the miRNA–mRNA network were also analyzed using the Genecodis database as described in Materials and Methods. The biological process GO annotation and KEGG pathway of target genes of DEMis with FDR <0.05 were considered as significant enrichment terms. As Figure 3A shows, skeletal system development and cell death biological processes including the apoptotic process and regulation of apoptotic process with FDR <0.05 were significantly enriched. The JAK–STAT signaling pathway, ECM-receptor interaction, Parkinson’s disease, and Alzheimer’s disease with FDR <0.05 were also significantly enriched (Figure 3B).

Using the starBase database, 26 DEMis were predicted to regulate 55 DELs; this was used to establish the lncRNA–miRNA–mRNA ceRNA network via integration with the miRNA–mRNA network (Figure 4). This network comprised 933 nodes (26 DEMis, 55 DELs, and 852 DEMs) and 2,801 interactions. Notably, down-regulated HOTAIR may function as a ceRNA to suppress the inhibitory effects of hsa-miR-454-3p on CRNKL1 and PDGFRB and hsa-miR-642a-5p on MMP13, MAP4K4, and PIK3R1, thus leading to their up-regulated expression. Similarly, up-regulated H19 may regulate the targeted effects of hsa-miR-454-3p on CRNKL1 and PDGFRB, as well as hsa-miR-2355-5p on PBX1 and TFDP2. Functional analysis of genes in the lncRNA-related ceRNA network revealed that they were significantly enriched in leukocyte transendothelial migration, focal adhesion, and ECM-receptor interaction.

Using the miRanda database, 5 DECs were predicted to regulate 17 DEMis; this information was used to establish the circRNA–lncRNA–mRNA ceRNA network via integration with the miRNA–mRNA network (Figure 5). Notably, up-regulated hsa_circRNA_001838 may function as a ceRNA to suppress the inhibitory effects of hsa-miR-4306 on FAM46A, VASH1, and CMTM6, thus resulting in their up-regulated expression. Functional analysis of genes in the lncRNA-related ceRNA network revealed that they were significantly enriched in leukocyte transendothelial migration, focal adhesion, and ECM–receptor interaction.

In our ROC analysis, the AUC of 2 out of 84 up-regulated DEMis was greater than 0.8, as shown in Figures 6A,B. The AUC of miR-1291 and miR-518b was 0.835 and 0.802, respectively. As shown in Figure 7, the AUC of 8 out of 74 down-regulated DEMis was greater than 0.8. The AUC of miR-1273e was greater than 0.9; the AUC of miR-623, miR-890, and miR-584-5p was 0.818; the AUC of miR-155-5p, miR-892b, miR-512-3p was 0.810; and the AUC of miR-454-3p was 0.802.

For DEMs, CRNKL1 was up-regulated in IDD patients in these three datasets with no significant difference (Figures 8A–C); PDGFRB was up-regulated in IDD patients in GSE150408 (Figure 8D) and GSE153761 (Figure 8E) datasets with no significant difference; IL1R1 had a trend of significant up-regulation in IDD patients in GSE124272 (Figure 8F); and CXCL12 (Figures 8H,I) was up-regulated in both GSE124272 and GSE153761 datasets with no significant difference. For DECs, only GSE150408 stores circRNA data of IDD patients. We found that hsa_circ_0001175 was up-regulated in IDD with no significant difference; in addition, hsa_circ_0000200, hsa_circ_0000926, and hsa_circ_0001838 were down-regulated in IDD with no significant difference (Figures 8J–M).