Skip to main content

Identification of key eRNAs for intervertebral disc degeneration by integrated multinomial bioinformatics analysis

Abstract

Background

Intervertebral disc degeneration (IVDD) is a common degenerative condition leading to abnormal stress distribution under load, causing intervertebral stenosis, facet joint degeneration, and foraminal stenosis. Very little is known about the molecular mechanism of eRNAs in IVDD.

Methods

Gene expression profiles of 38 annulus disc samples composed of 27 less degenerated discs (LDs) and 11 more degenerated discs (MDs) were retrieved from the GEO database. Then, differentially expressed enhancer RNAs (DEeRNAs), differentially expressed target genes (DETGs), and differentially expressed transcription factors (DETFs), hallmark of cancer signalling pathways according to GSVA; the types and quantity of immune cells according to CIBERSORT; and immune gene sets according to ssGSEA were analysed to construct an IVDD-related eRNA network. Then, multidimensional validation was performed to explore the interactions among DEeRNAs, DETFs and DEGs in space.

Results

A total of 53 components, 14 DETGs, 15 DEeRNAs, 3 DETFs, 5 immune cells, 9 hallmarks, and 7 immune gene sets, were selected to construct the regulatory network. After validation by online multidimensional databases, 21 interactive DEeRNA-DEG-DETF axes related to IVDD exacerbation were identified, among which the C1S-CTNNB1-CHD4 axis was the most significant.

Conclusion

Based upon the results of our study, we theorize that the C1S-CTNNB1-CHD4 axis plays a vital role in IVDD exacerbation. Specifically, C1S recruits CTNNB1 and upregulates the expression of CHD4 in IVDD, and subsequently, CHD4 suppresses glycolysis and activates oxidative phosphorylation, thus generating insoluble collagen fibre deposits and leading to the progression of IVDD. Overall, these DEeRNAs could comprise promising therapeutic targets for IVDD due to their high tissue specificity.

Peer Review reports

Introduction

Current therapies for intervertebral disc degeneration (IVDD) mainly include conservative and surgical treatments, which are supportive treatments but cannot cure IVDD [1, 2]. Investigations into biological regeneration treatments for IVDD have been increasing in recent years and have mainly focused on intervertebral disc regeneration [3, 4]. Nevertheless, because of the complex mechanisms of IVDD, these therapeutic techniques have not yet been widely utilized in the clinic. Hence, the identification of potential biomarkers that cause IVDD may help reveal the underlying mechanisms, providing novel targets for IVDD treatment.

Enhancer RNAs (eRNAs), long noncoding RNAs of 0.5-5 kb, are transcripts of highly histone-methylated enhancer regions [5]. Compared with general lncRNAs, eRNAs are mostly nonpolyadenylated and have a low stability, low abundance, short length, and short half-life; thus, eRNAs are considered nonfunctional transcriptional noise or transcriptional byproducts [6]. However, as cis-acting elements, eRNAs have been shown to recruit RNA polymerase II and transcription factors after being transcribed by enhancer units to regulate the expression of targeted or nearby genes [7, 8]. eRNAs, including HCG18 [9] (NR_024052.2), SNHG1 [10, 11] (NR_003098.2), H19 [12] (NR_002196.3), NEAT1 [13] (NR_028272.1), and linc-ADAMTS5 [14] (XR_007089974.1), were confirmed to play important roles in the occurrence and exacerbation of IVDD. Regulating the synthesis or degradation of the IVD extracellular matrix (ECM) is the most common method by which these eRNAs affect IVDD. For example, H19 and NEAT1 were upregulated in clinical IVDD specimens compared with normal IVD specimens. Both elevated H19 and NEAT1 promoted the expression of matrix metalloproteinases (MMPs) and ADAMTS5 (an aggrecan-degrading enzyme), while H19 could also elicit oxidative stress in nucleus pulposus (NP) cells, causing cellular senescence. Furthermore, linc-ADAMTS5 downregulation in IVDD led to the epigenetic silencing of ADAMTS5 via the recruitment of a transcription factor called RREB1 to the binding site overlapping the ADAMTS5 promoter. In addition to destabilizing the balance between ECM synthesis and degradation, inhibiting NP cell proliferation and promoting NP cell apoptosis are the mechanisms by which HCG18 and SNHG1 cause IVDD. However, these IVDD-related eRNAs were identified by lncRNA or mRNA microarrays, not systemic eRNA detection, and the functions of eRNAs have not been fully elucidated in IVDD. These findings indicated that eRNAs may play critical roles in IVDD; thus, these molecules could be utilized as biomarkers or therapeutic targets in clinical practice.

In this study, we constructed a regulatory network including 14 DETGs, 15 DEeRNAs, 3 DETFs, 5 immune cells, 9 hallmarks, and 7 immune gene sets. In addition, 21 interactive DEeRNA-DETF-DEG axes related to IVDD were identified. These DEeRNAs are promising therapeutic targets for IVDD due to their high tissue specificity.

Materials and methods

Data acquisition

The expression profiles and clinical information of 38 annulus disc samples consisting of 27 less degenerated discs (LDs) (Thompson grade I-III) and 11 more degenerated discs (MDs) (Thompson grade IV and V) samples were acquired from the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geoaccession numbers: GSE15227 [15] and GSE23130 [16]). These two sets of expression profiles were obtained via microarray analysis ([U133_X3P] Affymetrix Human X3P Array) and collected from the same platform, GPL1352. Expression profiles containing eRNAs and corresponding target genes in IVDD were extracted for subsequent analysis based on previous studies [17]. The expression profiles of immune-related genes and transcription factors (TFs) were retrieved from the ImmPort database [18] (http://www.immport.org/) and the Cistrome Cancer database [19] (http://cistrome.org/), respectively. ChIP-seq data for H3K27ac (accession numbers: GSM1195331 [20], GSM935653 [21] and GSM722419 [22]) and ATAC-seq data for eRNAs (accession number: GSE139099 [23]) were downloaded from the GEO database.

Data processing and differentially expressed gene analysis

Patients with incomplete clinical information were excluded from this study. All original microarray data were analysed by using the affy package in R (https://bioconductor.org/biocLite.R). The robust multiarray average (RMA) algorithm was utilized to correct the background of the microarray data and normalize the microarray data. According to the Thompson grading criteria, Grade I and II samples were considered nondegenerative, while Grade III to V samples were considered IVDD. DEeRNAs, DEGs, and DETFs between LDs and MDs were identified using the Linear Models for Microarray Data (limma) package [24] under |log2FC (fold change)| > 1 and a false discovery rate (FDR) < 0.05.

Functional enrichment analysis

Gene Ontology (GO) analysis was conducted by using the clusterProfiler R package, which revealed the functional enrichment of biology process (BP), cell component (CC), and molecular function (MF) terms. Reactome pathway analysis of eRNA-related coding genes was also conducted according to coexpression analysis. To avoid the accumulation of type I errors, enrichment items meeting FDR < 0.05 were considered significant. In addition, Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis was performed using the clusterProfiler R package, and KEGG pathways with an FDR < 0.05 were considered significant.

Immune cell infiltration analysis

The types and quantity of immune-infiltrating cells in LD and MD samples were estimated by cell type identification with the CIBERSORT algorithm [25] based on relative subsets of RNA transcripts. The immune-infiltrating cells with a CIBERSORT output of P < 0.05 were subsequently selected for nonparametric tests.

GSVA and ssGSEA

Significant immune cells and immune-related signalling pathways were evaluated and selected by ssGSEA [26] from 29 immune-associated gene sets, including a total of 707 genes. Potential downstream hallmarks of cancer signalling pathways of DEeRNAs were identified by GSVA [27].

Construction of the eRNA regulatory network

Pearson correlation analysis was performed for DETFs, DEeRNAs, DETGs, immune gene sets, immune cells, and hallmark pathways. Interaction pairs between DEeRNAs and other components were controlled based on a P value < 0.05 and |correlation coefficient| > 0.40. Then, with Cytoscape (v3.7.1) [28], a regulatory network was constructed based on significant interaction pairs between these components.

ATAC-seq validation

ATAC-seq data of eRNAs were downloaded from the GEO database to determine the open chromatin region of all the DEeRNAs on the chromosome where the corresponding enhancers were located.

Multidimensional validation

To decrease inherent deficiencies in silico, based on multivariate online databases, we chose data from tissues similar to the annulus disc, mainly tissues rich in fibroblasts, since there were no annulus disc samples available. Specifically, Pathway Card (https://pathcards.genecards.org/), KEGG, Expression Atlas [29], Genotype-Tissue Expression (GTEx) [30], String [31], and The Human Protein Atlas [32] databases were used.

External single cell RNA-seq validation

In our study, fibroblast data derived from the Single Cell Expression Atlas (https://www.ebi.ac.uk/gxa/sc/experiments/E-HCAD-13/results) were utilized to validate the expression of DEGs and DETFs at the single-cell level.

External ChIP-seq validation

Histone H3K27ac is involved in enhancer-specific modifications, which are critical for enhancers to promote the transcription of corresponding target genes [33]. The distribution of H3K27ac histones in T, B and K562 cells was analysed to determine enhancer-related TFs, and ChIP-seq and transcriptomic analyses were conducted to evaluate the regulatory roles of the TFs. The UCSC Genome Browser (http://genome-asia.ucsc.edu/cgi-bin/hgGateway?redirect=manual&source=genome.ucsc.edu) and the Cistrome Data Browser [19] (http://cistrome.org/) were the main instruments used to validate the chromatin accessibility of DETFs in the regions of DEeRNAs.

hTFtarget validation

Identification of TF-target gene regulatory relationships is essential for revealing the functions of TFs and their relative regulation of the expression of target genes. Here, we used the hTFtarget database [34] (http://bioinfo.life.hust.edu.cn/hTFtarget#!/) to predict the upstream TFs of the identified DETGs.

Hi-C validation

The Hi-C functional module of the 3D Genome Browser [35] (http://3dgenome.fsm.northwestern.edu/) was used to identify the topologically associating domains (TADs) of identified DEeRNAs through analysis of Hi-C data [36], which enabled us to validate the potential regulatory events between identified DEeRNAs and DETGs.

Statistical analysis

In our study, all statistical analysis processes were performed by R software version 3.6.1 (Institute for Statistics and Mathematics, Vienna, Austria). Statistical results were expressed as mean ± standard deviation (M ± SD). Data comparison of the two groups were analyzed using the Wilcoxon rank-sum test. Two-tailed p value < 0.05 was considered statistically significant. GraphPad Prism 7.0 was used to plot line charts.

Results

Identification of DETGs and DEeRNAs

The complete process of this research is illustrated in the flow chart (Fig. 1). The clinical information of the 38 IVDD patients is shown in Table S1. A total of 101 DEeRNAs (89 upregulated, 12 downregulated) and 555 DEGs (457 upregulated, 98 downregulated) were identified between LDs and MDs with a threshold of |log2 FC| > 1 and FDR P value < 0.05, as shown in heatmaps and volcano plots of DEeRNAs (Fig. 2A and B) and DEGs (Fig. 3A and B). The list of DEGs targeted by these DEeRNAs was acquired from the supplementary material of Zhao Zhang et al. [17]. Based on these results, a total of 75 DETGs were identified, and their expression was demonstrated via a heatmap (Figure S1 A) and volcano plot (Figure S1 B).

Fig. 1
figure 1

Flow diagram of the analytical process

Fig. 2
figure 2

Heatmap (A) and volcano plot (B) of eRNAs from the samples. Red and blue indicate the LD samples and MD samples, respectively. The green and red dots represent significantly differential expression

Fig. 3
figure 3

(A) Differentially expressed genes between LDs and MDs of annulus disc samples; red and blue indicate LD samples and MD samples, respectively. Volcano plot for genes from the samples (B); green and red dots represent significantly differential expression. Bubble plot for KEGG (C) and GO (D) enrichment analysis. Abbreviations: Gene Oncology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG)

Functional enrichment analysis

Based on GO and KEGG analyses of the potential functional mechanism of DETGs, the most important GO terms in the BP, CC, and MF of GO analysis were skeletal system development, extracellular matrix, and amide binding, respectively (Fig. 3C). The most significant KEGG pathway was amyotrophic lateral sclerosis (Fig. 3D).

Correlation analysis of DEeRNAs and DETGs

A total of 141 regulatory relationships between DEeRNAs and DETGs were identified by correlation analysis with the abovementioned thresholds of a |correlation coefficient|> 0.300 and a p value < 0.001. Based on the significant regulatory relationships of DEeRNAs and DETGs, complement C1s (C1S) (eRNA, enhancer ID: 12:6989331–6,995,331) and chromodomain helicase DNA binding protein 4 (CHD4) (DETG) were extracted for subsequent analysis (R = 0.915, P < 0.001, positive).

Identification of DETFs and correlation analysis of DEeRNAs and DETFs

Ten DETFs were identified by the limma package with the thresholds of |log2 FC| > 1 and FDR < 0.05, as shown in the heatmap (Fig. 4A) and volcano plot (Fig. 4B). A total of 63 regulatory relationships between DEeRNAs and DETFs were identified by correlation analysis based on the thresholds of a |correlation coefficient|> 0.300 and a P value < 0.001. The most significant DETF related to C1S was transcription factor 7 like 2 (TCF7L2) (R = 0.876, P < 0.001, positive), followed by CTNNB1 (R = 0.831, P < 0.001, positive) and TBL1XR1 (R = 0.825, P < 0.001, positive).

Fig. 4
figure 4

(A) DETFs between LD samples and MD samples; red and blue indicate LD samples and MD samples, respectively. Volcano plot (B) for DETFs from samples; green and red dots represent significantly differential expression. Heatmap (C) and volcano plot (D) showing the 50 hallmark gene sets acquired from MsigDB that were identified by GSVA between LD samples and MD samples. Heatmap (E) showing the GSVA scores of hallmark signalling pathways. Heatmap (F) of significant immune cells and signalling pathways related to immunological characteristics according to ssGSEA. Abbreviations: Quantitative gene set variation analysis (GSVA)

GSVA, ssGSEA, correlation analysis of DEeRNAs and hallmarks of cancer signalling pathways

Based on the gene expression of 27 LDs and 11 MDs, the enrichment levels of 50 hallmark pathways were calculated using GSVA, and the results are shown in a heatmap (Fig. 4C) and volcano plot (Fig. 4D). Specifically, based on the correlation analysis, EPITHELIAL_MESENCHYMAL_TRANSITION (R = 0.888, P < 0.001, positive), OXIDATIVE_PHOSPHORYLAION (R = 0.875, P < 0.001, positive) and MYC_TARGETS_V1 (R = 0.845, P < 0.001, positive) were significantly related to C1S in IVDD samples. The specific t value of the GSVA score of each signalling pathway between IVDD samples and nondegenerated samples is illustrated in Fig. 4E. In addition, 26 hallmark pathways were extracted for further analysis. The enrichment levels of 29 immune-associated gene sets were evaluated by ssGSEA, and the results are shown in Fig. 4F. Specifically, MHC class I, NK cells and iDCs were significantly enriched in IVDD samples.

Immune infiltration analysis and correlation analysis of DEeRNAs and immune cells

The infiltration levels of 22 types of immune cells in each sample were identified through the CIBERSORT algorithm based on relative subsets of RNA transcripts. The levels of immune cells in each sample are shown in the heatmap (Fig. 5A). A significant correlation between immune cells and IVDD was revealed by nonparametric tests (Fig. 5B and C), which indicated that the quantity of resting NK cells and resting mast cells in LDs was greater than that in MDs, and there were significant differences in the quantity of resting NK cells, M0 macrophages, activated myeloid dendritic cells and eosinophils to varying degrees in IVDD. The coexpression patterns among the identified immune cells are shown in Fig. 5D.

Fig. 5
figure 5

(A) The type and composition of immune cells estimated by the CIBERSORT algorithm in the samples. (B) The type of immune cells with significant differences between LD samples and MD samples; red and blue indicate MD samples and LD samples, respectively. (C) The types of immune cells with significant differences between the grades. (D) The coexpression patterns among identified immune cells

Based on the coanalysis of C1S and 22 types of immune cells/pathways, we found that M0 macrophages (R = 0.520, P < 0.01, positive), resting mast cells (R=-0.441, P < 0.05, negative), and resting CD4+ memory T cells (R=-0.375, P < 0.05, negative) were involved in IVDD progression.

Construction of the regulatory network

A total of 15 DEeRNAs, 14 DETGs, 3 DETFs, 5 immune cells, 9 hallmarks of cancer signalling pathways, and 7 immune-associated gene sets showing high correlations were extracted for subsequent analysis. The expression levels of these transcriptome components are illustrated in the heatmap (Fig. 6A). The components in the regulatory network are listed in Table S2. We performed correlation analysis between these eRNAs and the remaining 6 components, and a regulatory network was constructed based on significant correlation pairs (Fig. 6B). The results of the coexpression analysis of all 53 components in the network are displayed in Fig. 6C.

Fig. 6
figure 6

Construction of the DEeRNA network. Heatmap (A) for DEGs, DETFs and DEeRNAs on the final list from the selection process. Coheatmap (B) of all 6 dimensions and network (C) of all 6 dimensions on the final list, including eRNAs (red diamonds), hallmarks of cancer signalling pathways according to GSVA (blue rectangles), immune gene sets according to ssGSEA (green triangles), immune cells according to CIBERSORT (purple ellipses), DEeRNA target genes (pink hexagons) and DETFs (yellow arrowheads)

ATAC-seq validation

Figure 7 shows the open chromatin region of all 15 DEeRNAs on the chromosome where the corresponding enhancer is located, and the presence of colourful peaks indicates chromatin accessibility.

Fig. 7
figure 7

IGV plots of the chromatin accessibility landscape of all 15 DEeRNAs in the samples: C1S (A), COL1A2 (B), COX7A1 (C), DAD1 (D), DYNLL1 (E), FAM20 (F), HIF1A (G), HSPB1 (H), LAPTM4A (I), PAM (J), PLOD2 (K), PMEPA1 (L), PSAP (M), PTMA (N), and YWHAQ (O). The presence of colourful peaks indicates chromatin accessibility

Multidimensional validation

To decrease the bias of different algorithms and databases, we used multivariate external databases. Based on Pathway Card and KEGG analyses, PXN was related to HALLMARK_TGF_BETA_SIGNALING and HALLMARK_MTORC1_SIGNALING, and COX6B1 was involved in the oxidative phosphorylation signalling pathway, which is related to DETFs, CTNNB1 and TCF7L2. CTNNB1 also plays a role in the TGF-β signalling pathway. Moreover, DETFs, CTNNB1, TCF7L2 and TBL1XR1 are vitally important for the Wnt signalling pathway. Genotype-Tissue Expression (GTEx) was used to identify the expression levels of the DEGs and DETFs in multiple normal tissues (Figure S2). DepMap (Figure S3) shows the expression levels of DETGs and DETFs in multiple cancers, including samples homogeneous to the annulus disc, such as chondrosarcoma, and fibroblasts in multiple tissues.

External single cell RNA-seq validation

Based on the abundance of several subclusters of fibroblasts expressing high levels of fibroblast signature genes, such as CEMIP, AKR1C1, MGP, COMP, DNER, and MELTF [37, 38], the data of fibroblasts derived from the Single Cell Expression Atlas were used to validate the expression of DEGs and DETFs at the single-cell level. Importantly, in addition to CTNNB1, C1S and CHD4, the key gene vimentin in EPITHELIAL_MESENCHYMAL_TRANSITION and the key gene superoxide dismutase 1 (SOD1) in OXIDATIVE_PHOSPHORYLAION were coexpressed in fibroblasts (Fig. 8).

Fig. 8
figure 8

External single-cell RNA-seq validation by the Single Cell Expression Altas. A total of 13 key DEGs and 3 DETFs were expressed in the lung fibrosis samples, especially COX6B1, RPLP0, TMEM147, CHD4 and PXN, which were highly expressed

External ChIP-seq validation

The UCSC Genome Browser was used in our study to determine the locations of the DEeRNAs, and the chromatin accessibility (Fig. 9) of CTNNB1 (T lymphocytes), TCF7L2 (K562) and TBL1XR1 (B lymphocytes) was subsequently determined via the Cistrome Data Browser.

Fig. 9
figure 9

Chromatin accessibility of DETFs. The presence of black peaks indicates chromatin accessibility. (A, B, C) CTNNB1 binds to C1S, DYNLL1 and YWHAQ. (D, E) Chromatin accessibility of TBL1XR1 at the chromosomal locations of C1S and YWHAQ. (F, G, H, I) Chromatin accessibility of TCF7L2 at the chromosomal locations of C1S, COX7A1, PMEPA1 and YWHAQ

TF target validation

Using the hTFtarget database for the 3 DETFs and 15 DEGs in the network, we identified that CTNNB1 targets CHD4, and TCF7L2 targets CHD4, COX6B1, LIN37 and ASAP2 (Fig. 10). Moreover, COX6B1 was found to be the common target DEG for TBL1XR1 and TCF7L2. In the hTFtarget database, based on the GSM935574 and 935,678 data, there was a strong association between TBL1XR1 and TCF7L2 in the K562 cell line derived from bone marrow [the relative importance (RI) score = 37.68] (MAX value: 100, a higher score means that the two TFs have a stronger coassociation).

Fig. 10
figure 10

Validation of external TF targets via the hTFtarget database. Regulation and interaction locations of TCF7L2-ASAP2 (A), CTNNB1-CHD4 (B), TCF7L2-CHD4 (C), TCF7L2-COX6B1 (D) and TCF7L2-LIN37. The blue, green and red fishbone lines represent the coding, noncoding and problem states of genes, respectively. The short lines at the bottom of each graph indicate the position at which DETFs bind to genes, and the redness depth of the DETFs indicates the binding intensity

Hi-C validation

IMR90 is a type of fibroblast derived from foetal lung tissue that is histologically homogeneous to the annulus disc. Moreover, flavopiridol is a broad-spectrum inhibitor of cyclin-dependent kinases (CDKs) [39]. This finding is consistent with the significant upregulation of CDK 2-associated protein 1 (CDK2AP1) (logFC = 2.525, adjusted P value = 0.004) between the LD and MD groups in our study. CDK2AP1 is associated with CDK2 and is thought to suppress CDK2 function by sequestering monomeric CDK2 and then targeting CDK2 for proteolysis. Therefore, we used Hi-C sequencing of the IMR90_Flavopiridol cell line to validate the interaction between DETGs and DEeRNAs. We found that C1S and CHD4, VAMP1, ZNF384; DYNLLI and PXN, RNF34, RPLP0; PMEPA1 and RAE1, VAPB; COX7A1 and COX6B1, LIN37, TMEM147, ZNF146; YWHAQ and ASAP2 were close to and interacted with each other in three-dimensions (Fig. 11).

Fig. 11
figure 11

Analysis of Hi-C data from IMR90_Flavopiridol using the 3D Genome Browser. DYNLL1 and PXN, RNF34 and RPLP0 (A); C1S and CHD4, VAMP1, ZNF384 (B); COX7A1 and COX6B1, LIN37, TMEM147 and ZNF146 (C); PMEPA1 and RAE1, VAPB (D); YWHAQ and ASAP2 can interact with each other. The colour of the point where the upper corner of the triangle is located reflects the interaction of the genes through which the vertical line connects the two lower corners of the grey triangle

Discussion

We detected correlations between C1S and CTNNB1, C1S and CHD4, C1S and M0 macrophages, and C1S and OXIDATIVE_PHOSPHORYLAION. A total of 53 components, 14 DEGs, 15 eRNAs, 3 DETFs, 5 immune cells, 9 hallmarks of cancer signalling pathways, and 7 immune gene sets, were selected to construct the regulatory network. In addition to the C1S-CTNNB1-CHD4 axis, 20 other interactive DEeRNA-DETF-DEG axes are involved in IVDD exacerbation (Table 1), among which CTNNB1 targets CHD4 and TCF7L2 targets CHD4, COX6B1, LIN37 and ASAP2, as proven by the hTFtarget database. With respect to cis-acting elements, recent studies have shown that after being transcribed by enhancer units, eRNAs activate their target genes via overtranscription. Specifically, eRNAs located around the TFBS of target genes subsequently recruit RNA polymerase II, transcription factors or coactivators/transcriptional complexes to bind target gene TFBSs, which ultimately leads to the upregulation of their target genes. In our study, based on the functional pattern of eRNAs and the 3D interaction among these DEeRNAs, DETFs and DETGs, we speculated that C1S (eRNA location: 12:6989331–6,995,331) recruits CTNNB1 to increase the expression of CHD4 in IVDD.

Table 1 20 interactive “DEeRNA-DETF-DEG” axises

The intervertebral disc, especially the nucleus pulposus, is an immune-privileged tissue. An autoantigen triggers the immunological cascade once the physiological barrier is disrupted [40]. The impact of the complement system on IVDD was demonstrated in recent work. The complement system plays a role in IVDD by forming the terminal complement complex (TCC), an activator of inflammatory processes and cell lysis that has been shown to be highly deposited in degenerated intervertebral discs [41, 42]. In our study, we found that the gene fragment overlapping the DNA location of C1S, worsening IVDD, occurred not only by forming a TCC but also by acting as an eRNA, regulating some IVDD-related genes.

β-Catenin, encoded by CTNNB1, is an adhesion junction component that, together with cadherin and α-catenin, forms the adherens junction complex and plays an important role in the construction and maintenance of epithelial cell layers by regulating cell growth and intercellular adhesion [43]. Macrophages participate in the degradation and remodelling of the ECM through the production of MMPs (MMP-7, 9, 12) and A disintegrin and metalloproteinase with thrombospondin motifs [4447]. Furthermore, the chondrocytes and chondrocyte-like cells in the IVD secrete MMP-1 [48], MMP-2 [49], MMP-3 [50], MMP-9 [51] and MMP-13 [52]. The decomposition of the ECM is an important change in the development of IVDD. Interestingly, MMP-1, MMP-2, MMP-3, MMP-7, MMP-8, MMP-9, MMP-10, MMP-11, MMP-13, MMP-14, MMP-16, MMP-17, MMP-15, MMP-24, MMP-25, and MMP-27 were related to CTNNB1 (Figure S4), and these genes were verified to be involved in the pathogenesis of IVDD except for the latter four. These findings suggest that the eRNA C1S may recruit CTNNB1 to play a role in the mechanism of IVDD involving macrophages and MMPs.

CHD4 is an ATP-dependent chromatin remodeller involved in the epigenetic regulation of gene transcription, DNA repair, and cell cycle progression [53]. Recent studies have revealed a new function of CHD4 related to cell bioenergetics. CHD4 suppressed the expression of protein arginine deiminase 1 (PADI1) and PADI3, which are able to regulate citrullination of arginine residues of the allosterically regulated glycolytic enzyme pyruvate kinase M2 (PKM2) [54]. As a rate-limiting enzyme in glycolysis, citrullinated PKM2 slows glycolysis and alters the cell-bioenergetic balance between oxidative phosphorylation and glycolysis, which triggers and promotes the progression of fibrosis. Based on our findings on the roles of eRNAs and the C1S-CTNNB1-CHD4 axis in IVDD, we speculated that the eRNA C1S recruited CTNNB1 and upregulated the expression of CHD4 in IVDD and then suppressed glycolysis by modulating PKM2 and activating the oxidative phosphorylation process, thus resulting in insoluble collagen fibre deposits and leading to the genesis and exacerbation of IVDD. Deficience of experimental validation is one of the limitation of our study, and the subsequent validation of interactions among DEeRNAs, DETFs and DEGs, especially the C1S-CTNNB1-CHD4 axis in space in IVDD is worthy to carry out.

Conclusion

Based upon the results of our study, we theorize that the C1S-CTNNB1-CHD4 axis plays a vital role in IVDD exacerbation. Specifically, C1S recruits CTNNB1 and upregulates the expression of CHD4 in IVDD, and subsequently, CHD4 suppresses glycolysis and activates oxidative phosphorylation, thus generating insoluble collagen fibre deposits and leading to the progression of IVDD. Overall, these DEeRNAs could comprise promising therapeutic targets for IVDD due to their high tissue specificity.

Data availability

Publicly available datasets were analyzed in this study. This data can be found here: Gene Expression Omnibus (GSE15227, GSE23130, GSM1195331, GSM935653, GSM722419, GSE139099) (GEO, https://www.ncbi.nlm.nih.gov/geo/), ImmPort database (http://www.immport.org/), and Cistrome Cancer database (http://cistrome.org/).

References

  1. Evans CH, Huard J. Gene therapy approaches to regenerating the musculoskeletal system. Nat Rev Rheumatol. 2015;11(4):234–42. https://doi.org/10.1038/nrrheum.2015.28.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  2. Huang YC, Urban JP, Luk KD. Intervertebral disc regeneration: do nutrients lead the way? Nat Rev Rheumatol. 2014;10(9):561–6. https://doi.org/10.1038/nrrheum.2014.91.

    Article  PubMed  Google Scholar 

  3. Vickers L, Thorpe AA, Snuggs J, Sammon C, Le Maitre CL. Mesenchymal stem cell therapies for intervertebral disc degeneration: consideration of the degenerate niche. JOR Spine. 2019;2(2):e1055. https://doi.org/10.1002/jsp2.1055.

    Article  PubMed  PubMed Central  Google Scholar 

  4. Xia C, Zeng Z, Fang B, Tao M, Gu C, Zheng L, Wang Y, Shi Y, Fang C, Mei S, Chen Q, Zhao J, Lin X, Fan S, Jin Y, Chen P. Mesenchymal stem cell-derived exosomes ameliorate intervertebral disc degeneration via anti-oxidant and anti-inflammatory effects. Free Radic Biol Med. 2019;143:1–15. https://doi.org/10.1016/j.freeradbiomed.2019.07.026.

    Article  CAS  PubMed  Google Scholar 

  5. Chen H, Du G, Song X, Li L. Non-coding transcripts from enhancers: new insights into enhancer activity and gene expression regulation. Genomics Proteom Bioinf. 2017;15(3):201–7. https://doi.org/10.1016/j.gpb.2017.02.003.

    Article  CAS  Google Scholar 

  6. Ding M, Liu Y, Liao X, Zhan H, Liu Y, Huang W. Enhancer RNAs (eRNAs): New insights into Gene transcription and Disease Treatment. J Cancer. 2018;9(13):2334–40. https://doi.org/10.7150/jca.25829.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  7. Andersson R, Gebhard C, Miguel-Escalada I, Hoof I, Bornholdt J, Boyd M, Chen Y, Zhao X, Schmidl C, Suzuki T, Ntini E, Arner E, Valen E, Li K, Schwarzfischer L, Glatz D, Raithel J, Lilje B, Rapin N, Bagger FO, Jørgensen M, Andersen PR, Bertin N, Rackham O, Burroughs AM, Baillie JK, Ishizu Y, Shimizu Y, Furuhata E, Maeda S, Negishi Y, Mungall CJ, Meehan TF, Lassmann T, Itoh M, Kawaji H, Kondo N, Kawai J, Lennartsson A, Daub CO, Heutink P, Hume DA, Jensen TH, Suzuki H, Hayashizaki Y, Müller F, Forrest ARR, Carninci P, Rehli M, Sandelin A. An atlas of active enhancers across human cell types and tissues. Nature. 2014;507(7493):455–61. https://doi.org/10.1038/nature12787.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  8. Kim TK, Hemberg M, Gray JM, Costa AM, Bear DM, Wu J, Harmin DA, Laptewicz M, Barbara-Haley K, Kuersten S, Markenscoff-Papadimitriou E, Kuhl D, Bito H, Worley PF, Kreiman G, Greenberg ME. Widespread transcription at neuronal activity-regulated enhancers. Nature. 2010;465(7295):182–7. https://doi.org/10.1038/nature09033.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  9. Xi Y, Jiang T, Wang W, Yu J, Wang Y, Wu X, He Y. Long non-coding HCG18 promotes intervertebral disc degeneration by sponging miR-146a-5p and regulating TRAF6 expression. Sci Rep. 2017;7(1):13234. https://doi.org/10.1038/s41598-017-13364-6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  10. Sahu D, Hsu CL, Lin CC, Yang TW, Hsu WM, Ho SY, Juan HF, Huang HC. Co-expression analysis identifies long noncoding RNA SNHG1 as a novel predictor for event-free survival in neuroblastoma. Oncotarget. 2016;7(36):58022–37. https://doi.org/10.18632/oncotarget.11158.

    Article  PubMed  PubMed Central  Google Scholar 

  11. Tan H, Zhao L, Song R, Liu Y, Wang L. The long noncoding RNA SNHG1 promotes nucleus pulposus cell proliferation through regulating miR-326 and CCND1. Am J Physiol Cell Physiol. 2018;315(1):C21–7. https://doi.org/10.1152/ajpcell.00220.2017.

    Article  CAS  PubMed  Google Scholar 

  12. Wang X, Zou M, Li J, Wang B, Zhang Q, Liu F, Lü G. LncRNA H19 targets miR-22 to modulate H(2) O(2) -induced deregulation in nucleus pulposus cell senescence, proliferation, and ECM synthesis through wnt signaling. J Cell Biochem. 2018;119(6):4990–5002. https://doi.org/10.1002/jcb.26738.

    Article  CAS  PubMed  Google Scholar 

  13. Ruan Z, Ma H, Li J, Liu H, Jia H, Li F. The long non-coding RNA NEAT1 contributes to extracellular matrix degradation in degenerative human nucleus pulposus cells. Exp Biol Med (Maywood). 2018;243(7):595–600. https://doi.org/10.1177/1535370218760774.

    Article  CAS  PubMed  Google Scholar 

  14. Wang K, Song Y, Liu W, Wu X, Zhang Y, Li S, Kang L, Tu J, Zhao K, Hua W, Yang C. The noncoding RNA linc-ADAMTS5 cooperates with RREB1 to protect from intervertebral disc degeneration through inhibiting ADAMTS5 expression. Clin Sci (Lond). 2017;131(10):965–79. https://doi.org/10.1042/cs20160918.

    Article  CAS  PubMed  Google Scholar 

  15. Gruber HE, Hoelscher G, Loeffler B, Chow Y, Ingram JA, Halligan W, Hanley EN Jr. Prostaglandin E1 and misoprostol increase epidermal growth factor production in 3D-cultured human annulus cells. Spine J. 2009;9(9):760–6. https://doi.org/10.1016/j.spinee.2009.04.024.

    Article  PubMed  Google Scholar 

  16. Gruber HE, Hoelscher GL, Ingram JA, Hanley EN Jr. Genome-wide analysis of pain-, nerve- and neurotrophin -related gene expression in the degenerating human annulus. Mol Pain. 2012;8:63. https://doi.org/10.1186/1744-8069-8-63.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  17. Zhang Z, Lee JH, Ruan H, Ye Y, Krakowiak J, Hu Q, Xiang Y, Gong J, Zhou B, Wang L, Lin C, Diao L, Mills GB, Li W, Han L. Transcriptional landscape and clinical utility of enhancer RNAs for eRNA-targeted therapy in cancer. Nat Commun. 2019;10(1):4562. https://doi.org/10.1038/s41467-019-12543-5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. Bhattacharya S, Andorf S, Gomes L, Dunn P, Schaefer H, Pontius J, Berger P, Desborough V, Smith T, Campbell J, Thomson E, Monteiro R, Guimaraes P, Walters B, Wiser J, Butte AJ. ImmPort: disseminating data to the public for the future of immunology. Immunol Res. 2014;58(2–3):234–9. https://doi.org/10.1007/s12026-014-8516-1.

    Article  CAS  PubMed  Google Scholar 

  19. Zheng R, Wan C, Mei S, Qin Q, Wu Q, Sun H, Chen CH, Brown M, Zhang X, Meyer CA, Liu XS. Cistrome Data Browser: expanded datasets and new tools for gene regulatory analysis. Nucleic Acids Res. 2019;47(D1):D729–35. https://doi.org/10.1093/nar/gky1094.

    Article  CAS  PubMed  Google Scholar 

  20. van Loosdregt J, Fleskens V, Tiemessen MM, Mokry M, van Boxtel R, Meerding J, Pals CE, Kurek D, Baert MR, Delemarre EM, Gröne A, Koerkamp MJ, Sijts AJ, Nieuwenhuis EE, Maurice MM, van Es JH, Ten Berge D, Holstege FC, Staal FJ, Zaiss DM, Prakken BJ, Coffer PJ. Canonical wnt signaling negatively modulates regulatory T cell function. Immunity. 2013;39(2):298–310. https://doi.org/10.1016/j.immuni.2013.07.019.

    Article  CAS  PubMed  Google Scholar 

  21. Pope BD, Ryba T, Dileep V, Yue F, Wu W, Denas O, Vera DL, Wang Y, Hansen RS, Canfield TK, Thurman RE, Cheng Y, Gülsoy G, Dennis JH, Snyder MP, Stamatoyannopoulos JA, Taylor J, Hardison RC, Kahveci T, Ren B, Gilbert DM. Topologically associating domains are stable units of replication-timing regulation. Nature. 2014;515(7527):402–5. https://doi.org/10.1038/nature13986.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  22. Trompouki E, Bowman TV, Lawton LN, Fan ZP, Wu DC, DiBiase A, Martin CS, Cech JN, Sessa AK, Leblanc JL, Li P, Durand EM, Mosimann C, Heffner GC, Daley GQ, Paulson RF, Young RA, Zon LI. Lineage regulators direct BMP and wnt pathways to cell-specific programs during differentiation and regeneration. Cell. 2011;147(3):577–89. https://doi.org/10.1016/j.cell.2011.09.044.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  23. Grubert F, Srivas R, Spacek DV, Kasowski M, Ruiz-Velasco M, Sinnott-Armstrong N, Greenside P, Narasimha A, Liu Q, Geller B, Sanghi A, Kulik M, Sa S, Rabinovitch M, Kundaje A, Dalton S, Zaugg JB, Snyder M (2020) Landscape of cohesin-mediated chromatin loops in the human genome. Nature 583(7818):737–743. https://doi.org/10.1038/s41586-020-2151-x.

  24. Smyth GK. Linear models and empirical bayes methods for assessing differential expression in microarray experiments. Stat Appl Genet Mol Biol. 2004;3:Article3. https://doi.org/10.2202/1544-6115.1027.

    Article  PubMed  Google Scholar 

  25. Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, Hoang CD, Diehn M, Alizadeh AA. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods. 2015;12(5):453–7. https://doi.org/10.1038/nmeth.3337.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  26. Xiao B, Liu L, Li A, Xiang C, Wang P, Li H, Xiao T. Identification and Verification of Immune-related gene prognostic signature based on ssGSEA for Osteosarcoma. Front Oncol. 2020;10:607622. https://doi.org/10.3389/fonc.2020.607622.

    Article  PubMed  PubMed Central  Google Scholar 

  27. Ferreira MR, Santos GA, Biagi CA, Silva Junior WA, Zambuzzi WF. GSVA score reveals molecular signatures from transcriptomes for biomaterials comparison. J Biomed Mater Res A. 2021;109(6):1004–14. https://doi.org/10.1002/jbm.a.37090.

    Article  CAS  PubMed  Google Scholar 

  28. Kohl M, Wiese S, Warscheid B. Cytoscape: software for visualization and analysis of biological networks. Methods Mol Biol. 2011;696:291–303. https://doi.org/10.1007/978-1-60761-987-1_18.

    Article  CAS  PubMed  Google Scholar 

  29. Papatheodorou I, Fonseca NA, Keays M, Tang YA, Barrera E, Bazant W, Burke M, Füllgrabe A, Fuentes AM, George N, Huerta L, Koskinen S, Mohammed S, Geniza M, Preece J, Jaiswal P, Jarnuczak AF, Huber W, Stegle O, Vizcaino JA, Brazma A, Petryszak R. Expression Atlas: gene and protein expression across multiple studies and organisms. Nucleic Acids Res. 2018;46(D1):D246–51. https://doi.org/10.1093/nar/gkx1158.

    Article  CAS  PubMed  Google Scholar 

  30. Human genomics. The genotype-tissue expression (GTEx) pilot analysis: multitissue gene regulation in humans. Science. 2015;348(6235):648–60. https://doi.org/10.1126/science.1262110.

    Article  CAS  Google Scholar 

  31. Snel B, Lehmann G, Bork P, Huynen MA. STRING: a web-server to retrieve and display the repeatedly occurring neighbourhood of a gene. Nucleic Acids Res. 2000;28(18):3442–4. https://doi.org/10.1093/nar/28.18.3442.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  32. Uhlén M, Fagerberg L, Hallström BM, Lindskog C, Oksvold P, Mardinoglu A, Sivertsson Å, Kampf C, Sjöstedt E, Asplund A, Olsson I, Edlund K, Lundberg E, Navani S, Szigyarto CA, Odeberg J, Djureinovic D, Takanen JO, Hober S, Alm T, Edqvist PH, Berling H, Tegel H, Mulder J, Rockberg J, Nilsson P, Schwenk JM, Hamsten M, von Feilitzen K, Forsberg M, Persson L, Johansson F, Zwahlen M, von Heijne G, Nielsen J, Pontén F. Proteomics. Tissue-based map of the human proteome. Science. 2015;347(6220):1260419. https://doi.org/10.1126/science.1260419.

    Article  CAS  PubMed  Google Scholar 

  33. Kang Y, Kim YW, Kang J, Kim A. Histone H3K4me1 and H3K27ac play roles in nucleosome eviction and eRNA transcription, respectively, at enhancers. Faseb j. 2021;35(8):e21781. https://doi.org/10.1096/fj.202100488R.

    Article  CAS  PubMed  Google Scholar 

  34. Zhang Q, Liu W, Zhang HM, Xie GY, Miao YR, Xia M, Guo AY. hTFtarget: a comprehensive database for regulations of human transcription factors and their targets. Genomics Proteom Bioinf. 2020;18(2):120–8. https://doi.org/10.1016/j.gpb.2019.09.006.

    Article  Google Scholar 

  35. Wang Y, Song F, Zhang B, Zhang L, Xu J, Kuang D, Li D, Choudhary MNK, Li Y, Hu M, Hardison R, Wang T, Yue F. The 3D genome browser: a web-based browser for visualizing 3D genome organization and long-range chromatin interactions. Genome Biol. 2018;19(1):151. https://doi.org/10.1186/s13059-018-1519-9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  36. Jin F, Li Y, Dixon JR, Selvaraj S, Ye Z, Lee AY, Yen CA, Schmitt AD, Espinoza CA, Ren B. A high-resolution map of the three-dimensional chromatin interactome in human cells. Nature. 2013;503(7475):290–4. https://doi.org/10.1038/nature12644.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  37. Deroyer C, Charlier E, Neuville S, Malaise O, Gillet P, Kurth W, Chariot A, Malaise M, de Seny D. CEMIP (KIAA1199) induces a fibrosis-like process in osteoarthritic chondrocytes. Cell Death Dis. 2019;10(2):103. https://doi.org/10.1038/s41419-019-1377-8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  38. Marín YE, Seiberg M, Lin CB. Aldo-keto reductase 1 C subfamily genes in skin are UV-inducible: possible role in keratinocytes survival. Exp Dermatol. 2009;18(7):611–8. https://doi.org/10.1111/j.1600-0625.2008.00839.x.

    Article  CAS  PubMed  Google Scholar 

  39. Hah N, Murakami S, Nagari A, Danko CG, Kraus WL. Enhancer transcripts mark active estrogen receptor binding sites. Genome Res. 2013;23(8):1210–23. https://doi.org/10.1101/gr.152306.112.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  40. Ye F, Lyu FJ, Wang H, Zheng Z. The involvement of immune system in intervertebral disc herniation and degeneration. JOR Spine. 2022;5(1):e1196. https://doi.org/10.1002/jsp2.1196.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  41. Teixeira GQ, Yong Z, Goncalves RM, Kuhn A, Riegger J, Brisby H, Barreto Henriksson H, Ruf M, Nerlich A, Mauer UM, Ignatius A, Brenner RE, Neidlinger-Wilke C. Terminal complement complex formation is associated with intervertebral disc degeneration. Eur Spine J. 2021;30(1):217–26. https://doi.org/10.1007/s00586-020-06592-4.

    Article  PubMed  Google Scholar 

  42. Teixeira GQ, Yong Z, Kuhn A, Riegger J, Goncalves RM, Ruf M, Mauer UM, Huber-Lang M, Ignatius A, Brenner RE, Neidlinger-Wilke C. Interleukin-1β and cathepsin D modulate formation of the terminal complement complex in cultured human disc tissue. Eur Spine J. 2021;30(8):2247–56. https://doi.org/10.1007/s00586-021-06901-5.

    Article  PubMed  Google Scholar 

  43. Arbore C, Sergides M, Gardini L, Bianchi G, Kashchuk AV, Pertici I, Bianco P, Pavone FS, Capitanio M. α-catenin switches between a slip and an asymmetric catch bond with F-actin to cooperatively regulate cell junction fluidity. Nat Commun. 2022;13(1):1146. https://doi.org/10.1038/s41467-022-28779-7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  44. Burke B, Giannoudis A, Corke KP, Gill D, Wells M, Ziegler-Heitbrock L, Lewis CE. Hypoxia-induced gene expression in human macrophages: implications for ischemic tissues and hypoxia-regulated gene therapy. Am J Pathol. 2003;163(4):1233–43. https://doi.org/10.1016/s0002-9440(10)63483-9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  45. De Palma M, Venneri MA, Galli R, Sergi Sergi L, Politi LS, Sampaolesi M, Naldini L. Tie2 identifies a hematopoietic lineage of proangiogenic monocytes required for tumor vessel formation and a mesenchymal population of pericyte progenitors. Cancer Cell. 2005;8(3):211–26. https://doi.org/10.1016/j.ccr.2005.08.002.

    Article  CAS  PubMed  Google Scholar 

  46. Giraudo E, Inoue M, Hanahan D. An amino-bisphosphonate targets MMP-9-expressing macrophages and angiogenesis to impair cervical carcinogenesis. J Clin Invest. 2004;114(5):623–33. https://doi.org/10.1172/jci22087.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  47. Sun Q, Tian FM, Liu F, Fang JK, Hu YP, Lian QQ, Zhou Z, Zhang L. Denosumab alleviates intervertebral disc degeneration adjacent to lumbar fusion by inhibiting endplate osteochondral remodeling and vertebral osteoporosis in ovariectomized rats. Arthritis Res Ther. 2021;23(1):152. https://doi.org/10.1186/s13075-021-02525-8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  48. Grivas TB, Vasiliadis ES, Kaspiris A, Khaldi L, Kletsas D. Expression of matrix metalloproteinase-1 (MMP-1) in Wistar rat’s intervertebral disc after experimentally induced scoliotic deformity. Scoliosis. 2011;6(1):9. https://doi.org/10.1186/1748-7161-6-9.

    Article  PubMed  PubMed Central  Google Scholar 

  49. Mauth C, Bono E, Haas S, Paesold G, Wiese H, Maier G, Boos N, Graf-Hausner U. Cell-seeded polyurethane-fibrin structures–a possible system for intervertebral disc regeneration. Eur Cell Mater. 2009;18:27–38. https://doi.org/10.22203/ecm.v018a03. discussion 38 – 29.

    Article  CAS  PubMed  Google Scholar 

  50. Wang G, Huang K, Dong Y, Chen S, Zhang J, Wang J, Xie Z, Lin X, Fang X, Fan S. Lycorine suppresses endplate-chondrocyte degeneration and prevents intervertebral disc degeneration by inhibiting NF-κB signalling pathway. Cell Physiol Biochem. 2018;45(3):1252–69. https://doi.org/10.1159/000487457.

    Article  CAS  PubMed  Google Scholar 

  51. Zigouris A, Alexiou GA, Batistatou A, Voulgaris S, Kyritsis AP. The role of matrix metalloproteinase 9 in intervertebral disc degeneration. J Clin Neurosci. 2011;18(10):1424–5. https://doi.org/10.1016/j.jocn.2011.01.036.

    Article  CAS  PubMed  Google Scholar 

  52. Guo Z, Su W, Zhou R, Zhang G, Yang S, Wu X, Qiu C, Cong W, Shen N, Guo J, Liu C, Yang SY, Xing D, Wang Y, Chen B, Xiang H. (2021) Exosomal MATN3 of Urine-Derived Stem Cells Ameliorates Intervertebral Disc Degeneration by Antisenescence Effects and Promotes NPC Proliferation and ECM Synthesis by Activating TGF-β. Oxid Med Cell Longev 2021, 5542241. https://doi.org/10.1155/2021/5542241.

  53. Weiss K, Terhal PA, Cohen L, Bruccoleri M, Irving M, Martinez AF, Rosenfeld JA, Machol K, Yang Y, Liu P, Walkiewicz M, Beuten J, Gomez-Ospina N, Haude K, Fong CT, Enns GM, Bernstein JA, Fan J, Gotway G, Ghorbani M, van Gassen K, Monroe GR, van Haaften G, Basel-Vanagaite L, Yang XJ, Campeau PM, Muenke M. De Novo mutations in CHD4, an ATP-Dependent chromatin remodeler Gene, cause an intellectual disability syndrome with distinctive dysmorphisms. Am J Hum Genet. 2016;99(4):934–41. https://doi.org/10.1016/j.ajhg.2016.08.001.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  54. Coassolo S, Davidson G, Negroni L, Gambi G, Daujat S, Romier C, Davidson I. Citrullination of pyruvate kinase M2 by PADI1 and PADI3 regulates glycolysis and cancer cell proliferation. Nat Commun. 2021;12(1):1718. https://doi.org/10.1038/s41467-021-21960-4.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

Download references

Acknowledgements

Not applicable.

Funding

This study was funded by Natural Science Foundation of Shanghai (No.22ZR1450500); Excellence Program of Shanghai Municipal Health Commission (No.20234Z0020); Shanghai Rising-Star Program (Sailing Special Program) (No. 23YF1458400).

Author information

Authors and Affiliations

Authors

Contributions

RH designed the manuscript. YL wrote the manuscript. JY, DS, XH, TM, and HY revised the manuscript. RH and YL downloaded and analyzed data. RH and HY sponsored and provided the funding. All authors contributed to the article and approved the submitted version.

Corresponding authors

Correspondence to Tong Meng, Dianwen Song or Huabin Yin.

Ethics declarations

Ethical approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Competing interests

The authors declare no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Electronic supplementary material

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Li, Y., Huang, R., Ye, J. et al. Identification of key eRNAs for intervertebral disc degeneration by integrated multinomial bioinformatics analysis. BMC Musculoskelet Disord 25, 356 (2024). https://doi.org/10.1186/s12891-024-07438-6

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12891-024-07438-6

Keywords