Potential Role of lncRNAs in Contributing to Pathogenesis of Intervertebral Disc Degeneration Based on Microarray Data

Background

Intervertebral disc degeneration (IDD) is an important cause of low back pain, which affects patient quality of life and causes heavy socioeconomic burdens [1,2]. Currently, low back pain is mainly treated by medication combined with physiotherapy. Despite progress in surgical therapy, it is not curative and is related with several complications [3]. Fortunately, gene therapy has attracted more attention in recent years [4].

Several studies have unveiled a large number of genes or proteins that are implicated in IDD. For example, a former study has demonstrated that axin2 and Caveolin-1 are down-regulated in IDD [5]. Another study also suggests that expression of IL-1β is significantly higher in IDD [6], and IL-1β provokes the expression of brain-derived neurotrophic factor (BDNF) and vascular endothelial growth factor (VEGF), bringing about angiogenesis in human IDD [7]. Moreover, growth factors including BMP2 as well as TGF-β1 have recently been suggested to enhance the generation of proteoglycans by nucleus pulposus (NP) cells [8]. Degradation of metalloproteinases (MMPs) has been indicated to be connected with IDD [9]. Using a combination of a variety of growth factors has also been proposed to be appropriate to treat IDD, and platelet-rich plasma (PRP) is the source of these growth factors and has inhibitory effects on IDD in a rabbit model [10,11]. Targeting dysregulated genes has been suggested in animal IDD models, with inspiring results [12]. Nevertheless, dysregulation of gene expression is very complicated. Therefore, it is of great importance to elucidate the mechanisms underlying IDD at molecular levels.

Recently, aberrant expression of long non-coding RNAs (lncRNAs) is shown to cause disordered gene expression. LncRNAs are mRNA-like transcripts ranging in length from 200 nt to 100 kb, and they lack significant open reading frames [13,14]. In spite of this, accumulating studies have indicated that lncRNAs play crucial roles in many biological processes and human diseases [15,16]. For example, Wan et al. [17] demonstrated that over-expressed lncRNA RP11-296A18.3 induced the up-regulation of FAF1, which ultimately enhanced the abnormal apoptosis of disc cells. Wan et al. [17] also suggested the crucial roles of differentially expressed lncRNA in the nucleus pulposus of IDD, in which lncRNA-mRNA microarray profile GSE56081 was used. However, there are few studies investigating the role of lncRNAs in IDD. As a result, identifying the aberrant lncRNAs in IDD is still quite a challenge.

In our study, we intended to use the same microarray profile to further identify the differentially expressed lncRNA and mRNAs. Moreover, gene ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway, construction of lncRNA-mRNA weighted coexpression network, and protein-protein interaction (PPI) network were used to elucidate the molecular mechanism in the process of IDD. We believe that different results obtained from the microarray profile of Wan et al. might provide a profound understanding of IDD.

Material and Methods

TISSUE SAMPLES AND DATA ACQUISITION:

The gene expression data of GSE56081 deposited by Wan et al. [17] was downloaded from the National Center of Biotechnology Information (NCBI) Gene Expression Omnibus (http://www.ncbi.nlm.nih.gov/geo/) database, which was on the basis of GPL15314 platform of Arraystar Human LncRNA microarray V2.0 (Agilent_033010 Probe Name version). A total of 10 samples (5 human control nucleus pulposus tissues and 5 degenerative nucleus pulposus tissues) were included for the development of this microarray.

DATA PREPROCESSING:

Expression data of probes were converted to corresponding genes symbols according to the annotation of GPL15314 platform. Aggregate() function of R was used to compute the average expression value, which was used for genes corresponding to multi-probes, and KNN method [18] of Impute package of R [19] was used to add the missing values of probes.

Then, quantile normalization was carried out using preprocessCore package of R [20]. Subsequently, the expression matrix was obtained.

IDENTIFICATION OF DIFFERENTIALLY EXPRESSED LNCRNAS AND MRNAS:

Limma (linear models for microarray data) package [21] in R was used to identify differentially expressed lncRNAs and mRNAs between the 2 groups through use of the t test. Multiple testing correction was implement by calculating the Benjamini-Hochberg [22] false discovery rate (FDR). |log2 fold-change| >1 and an FDR <0.05 were regarded as the criteria for differential expression.

GENE ONTOLOGY (GO) AND PATHWAY ENRICHMENT ANALYSIS:

GO analysis is frequently used in functional enrichment studies of large-scale genes [23]. Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis was performed to analyze the biological pathways, involving the differentially expressed mRNAs. In the current study, TargetMine (http://targetmine.nibio.go.jp) [24] was used to investigate the functional enrichment condition for the up- and down-regulated differentially expressed mRNAs. FDR less than 0.05 was selected as the threshold.

THE CONSTRUCTION OF LNCRNA-MRNA WEIGHTED COEXPRESSION NETWORK:

WGCNA package [25] in R was used to construct the lncRNA/mRNA coexpression network. The network construction steps included the following: (1) Network construction: the lncRNA-mRNA weighted coexpression network is fully specified by its adjacency matrix amn, and amn encodes the network connection strength between nodes m and n. In order to calculate the adjacency matrix, the default approach defines the co-expression similarity Smn as the absolute value of the correlation coefficient between nodes of m and n: Smn=|cor(am, an)|. The weighted adjacency amn between 2 genes is proportional to their similarity on a logarithmic scale with β ≥0.8, log(amn)=β×log(Smn). Adjacency functions were obtained by approximate scale-free topology criterion. Then, the adjacency matrix was converted into topology matrix; (2) Module detection: we applied Dynamic Tree Cut method and Static Tree Cut method to detect modules with a minimum of 30 lncRNA/genes; and (3) Correlation analysis between disease information and gene module using the correlation coefficient method and gene significance (GS) value.

Then, we performed the GO and KEGG pathway enrichment analysis of genes in the module.

PPI NETWORK CONSTRUCTION:

The interaction of the proteins were analyzed by means of the online STRING database (Search Tool for the Retrieval of Interacting Genes, http://string-db.org/) [26] and the combined score >0.9 was used as the cut-off criterion. Subsequently, Cytoscape (http://cytoscapeweb.cytoscape.org/) [27] was used to display the PPI network.

Results

DATA PREPROCESSING:

On the basis of the annotation information of GPL56081, a total of 18 208 probes were identified, including 1970 lncRNAs and 16238 mRNAs. The data before and after normalization are shown in Figure 1A and 1B.

IDENTIFICATION OF DIFFERENTIALLY EXPRESSED LNCRNAS AND MRNAS:

Based on the cut-off criteria, 135 significantly up- and 170 down-regulated lncRNAs were screened out. Moreover, 2133 significantly up- and 1098 down-regulated mRNAs were observed, as shown in Figure 2.

GO AND PATHWAY ENRICHMENT ANALYSIS:

TargetMine was used to identify GO enriched functions for significant differentially expressed mRNAs. Top 5 GO and KEGG terms are listed in Table 1. Up-regulated genes were significantly enriched in macromolecule metabolic process (FDR=1.41E-03) and cellular macromolecule metabolic process (FDR=1.28E-03) of BP, membrane-bounded organelle (FDR=1.13E-08) and organelle (FDR=1.74E-07) of CC, RNA binding (FDR=1.31E-05), and nucleic acid binding (FDR=9.35E-04).

Significantly enriched KEGG pathways of up-regulated differentially expressed mRNAs were mainly on translation (FDR=2.24E-15) and cellular responses to stress (FDR=6.40E-10). However, the significant functions and pathways of down-regulated genes were not identified.

THE CONSTRUCTION OF LNCRNA-MRNA WEIGHTED COEXPRESSION NETWORK:

In the current study, a total of 7 weighted coexpression sub-networks were identified. The results of correlation analysis demonstrated that in addition to the grey module, the other 6 modules were highly relevant in intervertebral disc degeneration, as shown in Table 2 and Figure 3. Then, lncRNA-mRNA weighted coexpression network was constructed on the basis of the genes with top 30 connectivity degrees in 6 modules and P value <0.01. Among these, 8 lncRNAs were involved in the 4 modules (green module: CTC-523E23.5, RP4-639J15.1 and RP11-363G2.4; green-yellow module: CTD-2246P4.1 and LINC00917; magenta module: AC005082.12; purple module: MIR132 and RP11-38F22.1, Table 3) and were up-regulated (Figure 4A–4D). Notably, in these modules, red diamond nodes represented up-regulated lncRNA, pink square nodes represented up-regulated genes, and pastel green square nodes represented down-regulated genes (Figure 4).

Then, we performed the GO and KEGG pathway enrichment analysis of genes in the module. As shown in Table 4, the genes in the green module were significantly enriched in regulation of cell adhesion and response to axon injury. The genes in the green-yellow module were significantly enriched in positive regulation of cell migration and positive regulation of cell motion. The genes in the magenta module were significantly enriched in nucleotide catabolic process as well as nucleobase, nucleoside, nucleotide, and nucleic acid catabolic process. The genes in the purple module were significantly enriched in contractile fiber part and striated muscle thin filament.

PPI NETWORK:

The PPI network of significantly up-regulated genes was constructed using the STRING database. As shown in Figure 5, several PPI nodes were high in connectivity degrees, as follows: UBA52 (ubiquitin A-52 residue ribosomal protein fusion product 1, degree=51), RPS3 (ribosomal protein S3, degree=41), RPL9 (ribosomal protein L9, degree=38), RPL19 (ribosomal protein L19, degree=37), RPL23 (ribosomal protein L23, degree=36), RPL30 (ribosomal protein L30, degree=36), and RPS14 (ribosomal protein S14, degree=36).

Discussion

In this study, we identified 135 significantly up- and 170 down-regulated lncRNAs and 2133 significantly up- and 1098 down-regulated mRNAs. Additionally, up-regulated UBA52 had the highest connectivity degree in PPI network and was remarkably enriched in the pathway of metabolism of proteins. lncRNA-mRNA weighted coexpression network analysis showed that a total of 8 key lncRNAs were identified and were remarkably higher in IDD samples relative to normal tissues, for example, lncRNA LINC00917 and CTD-2246P4.1 were identified in the green-yellow module. These lncRNAs and their interacting genes may be responsible for IDD development and merit further consideration.

Significantly, up-regulated UBA52 had the highest connectivity degrees in PPI network and was remarkably enriched in the pathway of metabolism of proteins in our study. The intervertebral disc (IVD), a specialized fibrocartilaginous connective tissue, is mainly composed of the cell and extracellular matrix (ECM), and the molecular components of ECM are collagen, elastin, and proteoglycans [28]. Generally, the destruction of ECM is believed to be an important cause of IDD because the degradation of ECM help endothelial cells to migrate to further induce the generation of neovascularization [29,30]. MMPs are demonstrated to be the major proteolytic enzymes responsible for destruction of the ECM in the IVD [31]. UBA52 is a fusion protein composed of ubiquitin and ribosomal protein. Moreover, ubiquitin proteins are encoded by some genes, including UBA52 and UBA80, and are typically related with tagging proteins for proteosome degradation [32,33]. A previous study has demonstrated that extracellular ubiquitin increases MMP-2 protein level and enhances the migration of endothelial cells [34]. In light of these results, we infer that inhibiting the expression of UBA52 may be a good strategy in the treatment of IDD by reducing the ubiquitin proteins levels to further suppress the destruction of ECM.

In the green-yellow module, lncRNA LINC00917 and CTD-2246P4.1 were remarkably higher in IDD samples relative to normal tissues. Moreover, their interacted genes, such as up-regulated SPHK1 and down-regulated EGFR, were significantly enriched in positive regulation of cell migration. Several growth factors have been demonstrated to increase endothelial cell migration to further induce the generation of neovascularization [30]. Significantly, the loss of ECM and infiltration of blood vessels are 2 characteristics of IDD [28]. SPHK1 is a members of the SPHKs family, which has been indicated to play a crucial role in cell migration [35]. Moreover, growing evidence has demonstrated that SPHK1 participates in the production of angiogenesis. SPHK1 overexpression promoted endothelial cells survival in the absence of extracellular matrix and played a crucial role in vascular generation [36]. In the current study, the level of SPHK1 was up-regulated. Accordingly, SPHK1 might be a potential biomarker in the process of IDD through increasing endothelial cell migration to further induce the generation of neovascularization. Several previous studies have indicated that the abnormal expression of lncRNAs may dysregulate the expression of their neighboring genes [37]. Our findings suggest that the aberrant expression of SPHK1 may be caused by the up-regulated LINC00917 or CTD-2246P4.1, to further influence the occurrence of IDD. Therefore, lncRNAs LINC00917 and CTD-2246P4.1 may play a key role in the development of IDD via interacting with SPHK1. However, this assumption needs to be further investigated.

Furthermore, our study also identified the other 6 differentially expressed lncRNAs as well as hundreds of differentially expressed genes in different modules. CTC-523E23.5, RP4-639J15.1, and RP11-363G2.4 were identified and could interact with Integrin, Alpha 11 (ITGA11) in the green module. ITGA11 is shown to be differentially expressed and is associated with cell proliferation in IVD [38]. AC005082.12 was screened out and could interact with Ephrin-A3 (EFNA3) in the magenta module. MIR132 and RP11-38F22.1 could interact with Cathepsin L (CTSL) in the purple module. CTSL is demonstrated to degrade ECM and basement membrane components [39]. The destruction of ECM has been regarded to be an important cause of IDD [29,30]. Although the underlying mechanisms of these key lncRNAs in the development of IDD have not been fully discussed in our study, our results indicate that they may contribute to IDD development to some extent.

However, there were still some limitations in the current study. First, sample size was small which might limit the validity of the results. Second, no animal experiments were performed in our study to further validate these results. Thus, large sample sizes and more experiments, such as RNA interference and fluorescence in situ hybridization, are required in future work to more comprehensively elucidate the mechanisms of lncRNAs in the development of IDD.

Conclusions

lncRNAs LINC00917, CTD-2246P4.1, CTC-523E23.5, RP4-639J15.1, RP11-363G2.4, AC005082.12, MIR132, and RP11-38F22.1 were differentially expressed and might play important roles in the development of IDD. Key genes, such as UBA52 and SPHK1, may be pivotal biomarkers for IDD. Nevertheless, further experiments are required to validate their effects and mechanisms in IDD.