RETRACTED

This article has been retracted by Hindawi following an investigation undertaken by the publisher [1]. This investigation has uncovered evidence of one or more of the following indicators of systematic manipulation of the publication process:

(1)Discrepancies in scope(2)Discrepancies in the description of the research reported(3)Discrepancies between the availability of data and the research described(4)Inappropriate citations(5)Incoherent, meaningless and/or irrelevant content included in the article(6)Manipulated or compromised peer review

The presence of these indicators undermines our confidence in the integrity of the article’s content and we cannot, therefore, vouch for its reliability. Please note that this notice is intended solely to alert readers that the content of this article is unreliable. We have not investigated whether authors were aware of or involved in the systematic manipulation of the publication process.

Wiley and Hindawi regrets that the usual quality checks did not identify these issues before publication and have since put additional measures in place to safeguard research integrity.

We wish to credit our own Research Integrity and Research Publishing teams and anonymous and named external researchers and research integrity experts for contributing to this investigation.

The corresponding author, as the representative of all authors, has been given the opportunity to register their agreement or disagreement to this retraction. We have kept a record of any response received.

BioMed Research International
Volume 2021 (2021), Article ID 6696335, 33 pages
https://doi.org/10.1155/2021/6696335

Research Article

Transcriptional Profiling Uncovers Biologically Significant RNAs and Regulatory Networks in Nucleus Pulposus from Intervertebral Disc Degeneration Patients

Yuanyuan Chen, Bin Cai, Xiaofeng Lian, Jianguang Xu, and Tao Zhang

Department of Orthopaedics, Shanghai Jiao Tong University Affiliated Sixth People’s Hospital, Shanghai 200233, China

Correspondence should be addressed to Tao Zhang; zhangtao012345678@126.com

Received 20 October 2020; Revised 24 November 2020; Accepted 6 February 2021; Published 22 February 2021

Academic Editor: Parameshachari B. D.

Copyright © 2021 Yuanyuan Chen et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

Objective. This study aimed to uncover biologically significant RNAs in nucleus pulposus tissues of human intervertebral disc degeneration (IVDD) by integrated transcriptional profiling. Methods. From the Gene Expression Omnibus (GEO) database, three IVDD-related microarray profiling datasets were retrieved and assessed by intragroup data repeatability test. Then, differentially expressed circRNAs, lncRNAs, mRNAs, and miRNAs were screened in nucleus pulposus tissues between IVDD and control samples via the limma package. Coexpression networks were separately conducted via weighted gene correlation network analysis (WGCNA). Based on the feature RNAs in the IVDD-related modules, IVDD-related circRNA-miRNA-mRNA and lncRNA-miRNA-mRNA networks were conducted. The differentially expressed mRNAs in the two networks were analyzed by protein-protein interaction (PPI) and functional enrichment analyses. Results. By the intragroup data repeatability test, outlier samples were removed. Abnormally expressed RNAs were separately identified in nucleus pulposus between IVDD and controls. Via WGCNA, IVDD-related coexpression modules were constructed and the feature circRNAs, lncRNAs, mRNAs, and miRNAs were identified. Then, the circRNA- and lncRNA-miRNA-mRNA networks were built for IVDD. These mRNAs in the network exhibited complex interactions. Moreover, they were involved in distinct IVDD-related biological processes and pathways such as transcription, cell proliferation, TNF, TGF-β, and HIF signaling pathways. Conclusion. This study revealed biologically significant noncoding RNAs and their complex regulatory networks for IVDD.

1. Introduction

Intervertebral disc degeneration (IVDD) is a well-known cause of low back pain that is the main cause of disability [1]. IVDD is a multifactorial process, featured by changes in phenotype and genotype [2]. The nucleus pulposus is one of the main components of the intervertebral disc. Excessive apoptosis and extracellular matrix (ECM) degradation of nucleus pulposus cells (NPC) are the main pathological changes of IVDD [3]. The dysfunction of NPC facilitates the overproduction of proinflammatory factors in IVDD [4]. The inhibition of the abnormal activities of NPC can ameliorate the progression of IVDD [5]. Current treatments such as conservative treatment and surgery are mainly focused on pain relief, rather than treatment of the changes in pathology of degeneration [6]. Hence, a clear treatment strategy is urgently required. It is of importance to comprehensively comprehend the pathogenesis of IVDD, thereby accelerating the development of the therapeutic strategy concerning repairing the intervertebral disc injury.

Gene therapy has received extensive attention on account of its unique capacity to mediate cellular biological processes in IVDD, which provides an excellent synergistic treatment consequence by combining distinct genes [7]. Thus, it may be durable and beneficial to the therapy of IVDD. Targeting the specific noncoding RNAs (such as microRNAs (miRNAs), long noncoding RNA (lncRNA), and circular RNAs (circRNAs)) is a promising therapeutic strategy to lower toxicity and other harmful consequences [8]. miRNAs participate in various biological functions of IVDD cells by targeting different genes related to the IVDD process [9]. As a regulator of gene expression, miRNA has attracted widespread attention in the prevention and treatment of IVDD [10]. Different from linear RNAs, circRNAs with tissue specific and conservative features exhibit a closed circular structure, which are not affected by RNA exonuclease [11]. Thus, circRNAs are more stable and difficult to degrade. They are rich in binding sites of miRNA, as miRNA sponges to abrogate the inhibitory effect of miRNAs concerning their target genes during the progression of IVDD [12]. For instance, circERCC2 improves IVDD through mediating mitochondrial phagocytosis as well as apoptosis in the NPC via the miR-182-5p/SIRT1 axis [13]. lncRNA is abundantly expressed in human genomes, which may regulate key cell biological processes including IVDD cells [14]. lncRNAs can be used as a miRNA sponge to mediate the expression of their target genes in IVDD [15]. However, there is still a lack of systematic analysis of noncoding RNAs and their regulatory mechanisms in IVDD. Hence, this study probed into biologically significant RNAs and conducted their regulatory networks for IVDD, offering underlying therapeutic targets as well as molecular mechanisms concerning IVDD.

2. Materials and Methods

2.1. Microarray Datasets

IVDD-related microarray data were retrieved from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/gds/), including GSE67566, GSE56081, and GSE63492 datasets. The GSE67566 dataset contained circRNA expression profiles in human nucleus pulposus tissues from 5 IVDD patients and 5 normal controls on the GPL19978 platform [16]. The GSE56081 dataset was composed of lncRNA and mRNA profiling of nucleus pulposus tissues from 5 IVDD patients and 5 controls on the GPL15314 platform. The GSE63492 dataset included miRNA expression profiling in nucleus pulposus specimens derived from 5 IVDD patients and 5 controls on the GPL19449 platform.

2.2. Data Preprocessing and the Intragroup Data Repeatability Test

The raw data were preprocessed via the Linear Models for Microarray Data (limma) package v3.34.0 in R [17], which were normalized by quantile normalization as well as log² conversion. When a gene corresponded to multiple probes, the average value was taken as the expression value of the gene. Intragroup data repeatability test was presented via principal component analysis (PCA) and Pearson’s correlation analysis.

2.3. Differential Expression Analysis

Differentially expressed circRNAs (DEcircRNAs), lncRNAs (DElncRNAs), mRNAs (DEmRNAs), and miRNAs (DEmiRNAs) were screened by the limma package in R [17]. The criteria were set as adjusted and . The results were visualized into volcano plots and hierarchical clustering heat maps utilizing the limma or pheatmap packages in R. Furthermore, hierarchical clustering analysis was conducted on the basis of Euclidean distance.

2.4. Weighted Gene Correlation Network Analysis (WGCNA)

Coexpression analysis was presented via the WGCNA package in R [18]. Although the data were subjected to batch effect removal and normalization conversion, the possibility of outliers cannot be ruled out. Outlier samples were prone to adversely affect the analysis results of the network module, so it was necessary to identify and remove outliers before constructing the network. Since the outliers could introduce excessive deviations to WGCNA, they were removed before WGCNA. The pickSoftThreshold function in the WGCNA package was used to establish the soft threshold . The selected soft threshold met the following conditions: (1) the degree of fitting between the constructed network and the scale-free network reached a high level; (2) a smaller value was chosen when the degree of fit was up to a satisfactory level. The adjacency coefficient between gene i and gene j was calculated by exponential adjacency function , as follows: (where was the expression vector of gene and was obtained by the power of the Pearson correlation coefficient between and ). The topological overlap matrix was then calculated. dissTOM was calculated as follows: . The hClust function of the WGCNA package was applied to cluster the samples. The dynamic tree cutting method was utilized to identify and merge coexpressed gene modules with similar expression patterns. According to the genes in each module, module eigengene of each module was calculated, which was defined as the first principal component of the expression levels of all genes in the module, representing the overall level of gene expression in the module. Then, we calculated Pearson’s correlation coefficient between the module and the traits of samples.

2.5. The Competing Endogenous RNA (ceRNA) Network

DEmiRNA-DEmRNA and DElncRNA-DEmiRNA relationships were predicted via the StarBase v2.0 database (http://starbase.sysu.edu.cn/) [19]. Based on circBase (http://www.circbase.org/cgi-bin/downloads.cgi) [20], DEcircRNA-DEmiRNA pairs were analyzed. After integration of DEcircRNA-DEmiRNA, DElncRNA-DEmiRNA, and DEmiRNA-DEmRNA pairs, IVDD-related circRNA-miRNA-mRNA and lncRNA-miRNA-mRNA networks were separately conducted via the Cytoscape v3.7.2 software (http://www.cytoscape.org/) [21].

2.6. The Protein-Protein Interaction (PPI) Network

The Search Tool for the Retrieval of Interacting Genes (STRING) (string db.org) v10.0 database has been widely applied to analyze the interactions between proteins [22]. The mRNAs in the ceRNA networks were imported into the STRING database to study the interactions between the proteins encoded by genes. Their relationships were visualized by the Cytoscape v3.7.2 software (http://www.cytoscape.org/) [21].

2.7. Functional Annotation Analysis

Gene ontology (GO), covering biological process, cellular component, molecular function, and Kyoto Encyclopedia of Genes and Genomes (KEGG) were analyzed by using the Database for Annotation, Visualization, and Integrated Discovery (DAVID) online tool database v6.8 (http://david.abcc.ncifcrf.gov) [23]. was significantly enriched.

3. Results

3.1. The Intragroup Data Repeatability Test

This study retrieved the circRNA, lncRNA/mRNA, and miRNA expression datasets from human nucleus pulposus tissues separately derived from 5 IVDD patients and 5 normal control in the GEO database. To avoid the existence of abnormal samples from adversely affecting subsequent analysis, PCA was presented to visualize the clustering of samples of microarrays. The results showed that no abnormal samples were found in the GSE67566 (Figure 1(a)) and GSE56081 (Figure 1(b)). All samples from the two datasets were obviously divided into control and IVDD groups. In the GSE63492 dataset, there was one outlier sample (GSM1551029) that was removed (Figure 1(c)). Next, we performed correlation analysis on 10 samples in the three datasets. In the GSE67566 circRNA expression dataset, there were high consistencies between the samples in the IVDD group or the control group (Figure 1(d)). Furthermore, the correlation of samples between the IVDD group was higher than the correlation between the IVDD group and the control group in the GSE56081 lncRNA/mRNA dataset (Figure 1(e)). In Figure 1(f), we found that there was a case of an IVDD sample (GSM1551029) that was highly correlated with the samples in the control group from the GSE63492 miRNA expression dataset. Thus, the GSM1551029 sample was eliminated in the subsequent analysis. Figure 1(g) showed the correlation between the samples in the IVDD and control groups.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

Figure 1: PCA and Pearson’s correlation analysis. (a) PCA results for samples of circRNA expression profiles in the GSE67566 dataset. (b) PCA results for samples of mRNA and lncRNA expression profiles in the GSE56081 dataset. (c) PCA results for samples of miRNA expression profiles in the GSE63492 dataset. Blue is indicative of control samples and yellow is indicative of IVDD samples. Each point represents a sample. Sample correlation results in the (d) GSE67566, (e) GSE56081, and (f) GSE63492 datasets. (g) Sample correlation results in the GSE63492 dataset after excluding the GSM1551029 sample. The more the color tends to be redder, the stronger the correlation; on the contrary, the more the color tends to be greener, the weaker the correlation.

3.2. Identification of Differentially Expressed circRNAs, lncRNAs, mRNAs, and miRNAs for IVDD Nucleus Pulposus Tissues

Under the criteria of adjusted and , we carried out differential expression analysis in nucleus pulposus tissues between IVDD and control groups. In Figure 2(a), there were 105 DEcircRNAs (47 up- and 58 downregulated) between IVDD and control groups (Supplementary Table 1). Furthermore, 131 DElncRNAs (96 up- and 35 downregulated) and 928 DEmRNAs (687 up- and 241 downregulated) were identified in nucleus pulposus tissues between IVDD and control groups (Figure 2(b); Supplementary Table 2). 62 DEmiRNAs (31 up- and 31 downregulated) were screened in nucleus pulposus tissues between IVDD and control groups (Figure 2(c); Supplementary Table 3). Hierarchical clustering analysis revealed that DEcircRNAs (Figure 2(d)), DElncRNAs/mRNAs (Figure 2(e)), and DEmiRNAs (Figure 2(f)) could distinctly distinguish IVDD samples from control samples.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 2: Identification of differentially expressed circRNAs, lncRNAs, mRNAs, and miRNAs for IVDD nucleus pulposus tissues. Volcano plots depicting all differentially expressed (a) circRNAs, (b) lncRNAs/mRNAs, and (c) miRNAs in nucleus pulposus tissues between IVDD and control groups. Blue is indicative of downregulation and red signifies upregulation in the IVDD group. Hierarchical clustering analysis results for differentially expressed (d) circRNAs, (e) lncRNAs/mRNAs, and (f) miRNAs between IVDD and control groups. Green is indicative of downregulation and red signifies upregulation in the IVDD group.

3.3. An IVDD-Related Coexpression Network in the GSE67566 Dataset

5 IVDD nucleus pulposus samples and control samples in the GSE67566 dataset were used for WGCNA. In Figure 3(a), no outlier samples were detected among them, which was consistent with PCA results. Using the pickSoftThreshold function, we calculated the degree of scale-free topology model fit (Figure 3(b)) and mean connectivity (Figure 3(c)) when the value of was from 1 to 30. Herein, was chosen. Via the adjacency function, the adjacency coefficient between genes was calculated. TOM was calculated via the TOMsimilarity function, followed by the calculation of dissTOM. Based on the dissTOM, the hierarchical clustering graphs of genes were generated. In this study, the dynamic tree cutting method was utilized to identify and merge 3 coexpressed gene modules with similar expression patterns (Figure 3(d)). As shown in Figure 3(e), the genes in the turquoise module were significantly related to IVDD (). There were 2726 feature genes in the turquoise module (Supplementary Table 4), which contained all 105 DEcircRNAs.

(a)

(b)

(c)

(d)

(e)

Figure 3: An IVDD-related coexpression network in the GSE67566 dataset. (a) Sample hierarchical clustering tree and corresponding sample types (IVDD or control). The height of the clustering tree was set as 100. Dark-blue represents the IVDD group and cyan-blue suggests the control group. (b) Scale-free topology model fit under different values. (c) Mean connectivity under a series of values. (d) Coexpressed network construction. (e) Correlation between coexpression modules and IVDD. The darker the color, the greater the correlation. Red indicates positive correlation and green is indicative of negative correlation.

3.4. An IVDD-Related Coexpression Network in the GSE56081 Dataset

Consistent with our PCA results, there was no outlier sample in the GSE56081 dataset (Figure 4(a)). The optimal soft threshold was determined under different values (1–20). Combining the scale-free topology model fit (Figure 4(b)) and mean connectivity (Figure 4(c)), was determined. In the hierarchical clustering tree, 9 coexpression modules were merged via the dynamic cutting tree method (Figures 4(d)–4(f)). Pearson’s correlation between modules and IVDD was analyzed. In Figure 4(g), the genes in the blue ( and ) and turquoise () modules had a significant correlation with IVDD. The feature genes in the blue () and turquoise () modules were listed in Supplementary Tables 5 and 6.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

Figure 4: An IVDD-related coexpression network in the GSE56081 dataset. (a) Sample dendrogram and corresponding sample types (IVDD or control). The height of the clustering tree was set as 200. Dark-blue represents the IVDD group and cyan-blue suggests the control group. (b) Scale-free topology model fit with distinct values (1–20). The red line indicates an value of 0.90. (c) Mean connectivity with distinct values. (d–f) A coexpressed network containing 9 modules. (g) Pearson’s correlation between coexpression modules and IVDD. The darker the color, the greater the correlation. Red suggests positive correlation and green suggests negative correlation.

3.5. An IVDD-Related Coexpression Network in the GSE63492 Dataset

9 samples in the GSE63492 dataset were used for coexpression analysis. No outliers were found, as shown in Figure 5(a). On the basis of scale-free topology model fit (Figure 5(b)) and mean connectivity (Figure 5(c)), 10 was the optimal threshold of . In Figure 5(d), 4 coexpression modules were conducted by the dynamic cutting method. Among them, the genes in the blue module were distinctly correlated to IVDD ( and ) in Figure 5(e). The blue module contained a total of 537 feature genes (Supplementary Table 7).

(a)

(b)

(c)

(d)

(e)

Figure 5: An IVDD-related coexpression network in the GSE63492 dataset. (a) Sample hierarchical clustering tree and corresponding sample types (IVDD or control). The height of clustering tree was set as 100. Dark-blue represents the IVDD group and cyan-blue suggests the control group. (b) Scale-free topology model fit under different values. (c) Mean connectivity under a series of values. (d) Coexpressed network construction. (e) Correlation between coexpression modules and IVDD. The darker the color, the greater the correlation. Red indicates positive correlation and green is indicative of negative correlation.

3.6. Construction of a circRNA-miRNA-mRNA Network

After overlapping feature genes in the blue and turquoise modules and DEmRNAs in the GSE56081 dataset, 1456 feature DEmRNAs were obtained (Figure 6(a)). Furthermore, a total of 35 feature DEmiRNAs were obtained by the intersection of feature genes in the blue module and DEmiRNAs in the GSE63492 dataset (Figure 6(b)). 105 DEcircRNAs were all included in the feature circRNAs. Through the StarBase database, the target mRNAs of these feature DEmiRNAs were predicted. As a result, we only obtained the target mRNAs of 8 miRNAs (hsa-miR-519d-3p, hsa-miR-489-3p, hsa-miR-486-5p, hsa-miR-431-5p, hsa-miR-4306, hsa-miR-3196, hsa-miR-193a-5p, and hsa-miR-155-5p), as shown in Figure 6(c). Then, we took the intersection between the target mRNAs of 8 miRNAs and DEmRNAs. Consequently, 462 target DEmRNAs were identified, as shown in Figure 6(d). Moreover, circRNAs that targeted the 8 miRNAs were predicted using the circBase database. In Figure 6(e), 68 DEcircRNAs were intersected, which could be as the sponge of the 8 miRNAs. Based on the circRNA-miRNA and miRNA-mRNA pairs, we conducted a circRNA-miRNA-mRNA ceRNA network, composed of 68 DEcircRNAs, 8 DEmiRNAs, and 462 DEmRNAs (Figure 6(f)).

(a)

(b)

(c)

(d)

(e)

(f)

Figure 6: Construction of a circRNA-miRNA-mRNA network. (a) Venn diagram showing the overlap between feature genes in the blue and turquoise modules and DEmRNAs in the GSE56081 dataset. (b) Venn diagram depicting the overlap between feature genes in the blue module and DEmiRNAs in the GSE63492 dataset. (c) The number of targeted mRNAs of 8 DEmiRNAs. (d) The intersection between targeted mRNAs of 8 DEmiRNAs and DEmRNAs. (e) The intersection between circRNAs that targeted DEmiRNAs and DEcircRNAs. (f) The circRNA-miRNA-mRNA network for IVDD. The hexagon represents circRNA, the triangle represents miRNA, and the circle represents mRNA.

3.7. A PPI Network and Function Enrichment Analysis for DEmRNAs in the circRNA-miRNA-mRNA Network

To explore the intersections between proteins, we conducted a PPI network via the STRING database based on the DEmRNAs in the circRNA-miRNA-mRNA network. There were 212 nodes in the network (Figure 7(a)). Among them, 167 DEmRNAs were upregulated and 45 were downregulated in IVDD in comparison to controls. Furthermore, we probed into the biological functions of these DEmRNAs in the ceRNA network. In Figure 7(b), these genes were significantly related to key biological processes such as regulation of transcription, protein ubiquitination, and cell proliferation. They had a close relationship with multiple cellular components like nucleus, cytoplasm, and cytosol (Figure 7(c)). As shown in Figure 7(d), these genes exhibited distinct molecular functions such as protein binding, DNA binding, RNA binding, and transcription binding. For KEGG enrichment analysis results, critical pathways were markedly enriched like ubiquitin-mediated proteolysis, ribosome, FoxO signaling pathway, cell cycle, TNF signaling pathway, TGF-β signaling pathway, and HIF signaling pathway (Figure 7(e)).

(a)

(b)

(c)

(d)

(e)

Figure 7: A PPI network and function enrichment analysis for DEmRNAs in the circRNA-miRNA-mRNA network. (a) A PPI network based on DEmRNAs in the circRNA-miRNA-mRNA network. The red circle is indicative of upregulated mRNA in IVDD and blue is indicative of downregulated mRNA in IVDD. (b) The top ten biological processes of DEmRNAs in the ceRNA network. (c) The top ten cellular components of DEmRNAs in the ceRNA network. (d) The top ten molecular functions of DEmRNAs in the ceRNA network. (e) The top ten KEGG pathways of DEmRNAs in the ceRNA network. The more the color tends to be redder, the smaller the value. The size of the circle is proportional to the number of enriched genes.

3.8. Construction of a lncRNA-miRNA-mRNA Network for IVDD

To construct an IVDD-related lncRNA-miRNA-mRNA network, we analyzed the DElncRNA-DEmiRNA and DEmiRNA-DEmRNA pairs. Firstly, we obtained 15 feature DElncRNAs by overlapping DElncRNAs in the GSE56081 dataset and predicted lncRNAs of 8 DEmiRNAs (Figure 8(a)). In Figure 8(b), a total of 437 DEmRNAs that were targeted by 8 DEmiRNAs were identified by intersection of the targeted mRNAs of 8 DEmiRNAs and DEmRNAs in the GSE56081 dataset. Following the integration of lnRNA-miRNA and miRNA-mRNA pairs, a lncRNA-miRNA-mRNA network was built for IVDD (Figure 8(c)), composed of 15 lncRNAs, 7 miRNAs, and 437 mRNAs. A PPI network was then constructed through DEmRNAs in the lncRNA-miRNA-mRNA network (Figure 8(d)). In the network, there were 206 nodes, including 163 up- and 43 downregulated mRNAs in IVDD. The underlying functions of these DEmRNAs in the ceRNA network were then analyzed. Consequently, several critical biological processes were markedly enriched including the nodal signaling pathway, cell-cell junction organization, protein ubiquitination, and cell proliferation and transcription (Figure 8(e)). Moreover, they could be involved in the regulation of various cell components like cytoplasm, nucleoplasm, and nucleus (Figure 8(f)). Their molecular functions were predicted such as protein binding, transcription factor activity, and RNA binding (Figure 8(g)). KEGG enrichment analysis results suggested that these genes were distinctly related to IVDD-related pathways including the FoxO signaling pathway, TNF signaling pathway, cell cycle, TGF-β signaling pathway, and HIF signaling pathway (Figure 8(h)).

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Figure 8: Construction of a lncRNA-miRNA-mRNA network for IVDD. (a) Venn diagram showing the intersection between targeted lncRNAs of 8 DEmiRNAs and DElncRNAs in the GSE56081 dataset. (b) The overlap between predicted mRNAs of 8 DEmiRNAs and DEmRNAs in the GSE56081 dataset. (c) A lncRNA-miRNA-mRNA network for IVDD. (d) A PPI network based on the DEmRNAs in the ceRNA network. The red circle suggests upregulated mRNA in IVDD and the blue circle indicates downregulated mRNA in IVDD. (e) The top ten biological processes of DEmRNAs in the ceRNA network. (f) The top ten cellular components of DEmRNAs in the ceRNA network. (g) The top ten molecular functions of DEmRNAs in the ceRNA network. (h) The top ten KEGG pathways of DEmRNAs in the ceRNA network. The more the color tends to be redder, the smaller the value. The size of the circle is proportional to the number of enriched genes.

4. Discussion

The pathological mechanisms of IVDD remain ill-defined [24]. There is currently no effective therapy for reversal of IVDD progression [25]. Thus, it is of necessity to develop new therapeutic targets for IVDD. This study comprehensively illustrated the biologically significant RNAs in nucleus pulposus for IVDD, which could be potential targets for IVDD. Their regulatory networks were conducted, revealing underlying molecular mechanisms during the progression of IVDD.

The progression of IVDD is an extremely complex biological process. Studies have shown that mRNAs, lncRNAs, miRNAs, and circRNAs are all involved in the development of IVDD [12]. This study explored differentially expressed RNAs in nucleus pulposus tissues between IVDD and control samples based on transcriptomics [12]. With the criteria of adjusted and , 105 DEcircRNAs, 131 DElncRNAs, DEmRNAs, and 62 DEmiRNAs were identified for IVDD, which represented key genes during the progression of IVDD.

The construction of coexpression modules through the WGCNA algorithm to find potential key genes has been widely applied in various diseases [26–28]. For example, through WGCNA, key modules and genes related to non-small-cell lung cancer have been identified [29]. The genes that are divided into the same module exhibiting the same molecular functions. By calculation of the correlation between the module and the trait, we can determine the modules related to the trait and lock the feature genes. In this study, we applied WGCNA to find feature circRNAs, lncRNAs, mRNAs, and miRNAs related to IVDD development. By integration of differentially expressed circRNAs, lncRNAs, mRNAs, and miRNAs, feature RNAs have been explored for IVDD. It has been confirmed that circRNAs and lncRNAs could serve as the sponge of miRNAs, thereby regulating the expression of downstream target genes [30]. In the circRNA-miRNA-mRNA ceRNA network, there were 68 DEcircRNAs, 8 DEmiRNAs, and 462 DEmRNAs. For the lncRNA-miRNA-mRNA network, there were 15 DElncRNAs, 7 DEmiRNAs, and 437 DEmRNAs. Their regulatory mechanisms should be validated in cellular levels in depth.

We further explored the biological functions of DEmRNAs in the circRNA- and lncRNA-miRNA-mRNA networks. Our results showed that these mRNAs were mainly enriched in transcription, cell proliferation, TNF, TGF-β, and HIF signaling pathways. The inflammatory process exacerbated by TNF-α is a key mediator of IVDD [31]. TNF-α participates in distinct pathological processes of IVDD like inflammation, apoptosis, and proliferation [32]. TGF-β pathway activation is a promising treatment strategy for IVDD. However, its excessive activation could accelerate IVDD progression. Our findings offered more clues on its specific molecular mechanisms [33]. Furthermore, hypoxia affects the synthesis of intervertebral disc cells and weakens the ability to support the extracellular matrix of the intervertebral disc [34]. Hence, these DEmRNAs may be involved in the progression of IVDD via various biological processes and pathways.

However, several limitations should be pointed out in this study. Firstly, the sample size was not insufficient. In our future studies, the above findings of biologically significant noncoding RNAs (circRNAs, lncRNAs, and miRNAs) and mRNAs will be validated in a larger IVDD cohort. Secondly, our study is a preliminary screening study. Hence, we will verify the biological functions and interactions of these noncoding RNAs and mRNAs in IVDD by a series of experiments. Collectively, our findings identified IVDD-related RNAs and constructed two complex ceRNA networks, which offered potential therapeutic targets as well as molecular mechanisms for IVDD.

5. Conclusion

Taken together, this study identified biologically significant noncoding RNAs (circRNAs, lncRNAs, and miRNAs) and mRNAs. The regulatory networks involved in them were conducted for IVDD, indicating the complex molecular mechanisms. These DEmRNAs participated in distinct IVDD-related biological processes or pathways. Hence, these RNAs could be promising therapeutic targets for IVDD.

Abbreviations

IVDD: | Intervertebral disc degeneration
GEO: | Gene Expression Omnibus
WGCNA: | Weighted gene correlation network analysis
PPI: | Protein-protein interaction
ECM: | Extracellular matrix
NPC: | Nucleus pulposus cells
miRNAs: | MicroRNAs
lncRNA: | Long noncoding RNA
circRNAs: | Circular RNAs
limma: | Linear Models for Microarray Data
PCA: | Principal component analysis
DEcircRNAs: | Differentially expressed circRNAs
DElncRNAs: | Differentially expressed lncRNAs
DEmRNAs: | Differentially expressed mRNAs
DEmiRNAs: | Differentially expressed miRNAs
ceRNA: | Competing endogenous RNA
STRING: | Search Tool for the Retrieval of Interacting Genes
GO: | Gene ontology
KEGG: | Kyoto Encyclopedia of Genes and Genomes
DAVID: | Database for Annotation, Visualization, and Integrated Discovery.

Data Availability

The data used to support the findings of this study are included within the supplementary information files.

Conflicts of Interest

The authors declare no conflicts of interest.

Acknowledgments

This work was funded by Shanghai Science and Technology Commission Project (19JC1411702) and Retrospective Clinical Research Project of Shanghai Sixth People’s Hospital (ynhg201901).

Supplementary Materials

Supplementary 1. Supplementary Table 1. 105 differentially expressed circRNAs in nucleus pulposus tissues between IVDD and control groups.

Supplementary 2. Supplementary Table 2. 131 differentially expressed lncRNAs and 928 differentially expressed mRNAs in nucleus pulposus tissues between IVDD and control groups.

Supplementary 3. Supplementary Table 3. 62 differentially expressed miRNAs in nucleus pulposus tissues between IVDD and control groups.

Supplementary 4. Supplementary Table 4. 2726 genes in the turquoise module for the GSE67566 dataset.

Supplementary 5. Supplementary Table 5. 2726 genes in the blue module for the GSE56081 dataset.

Supplementary 6. Supplementary Table 6. 2726 genes in the turquoise module for the GSE56081 dataset.

Supplementary 7. Supplementary Table 7. 537 feature genes in the blue module in the GSE63492 dataset.

References

W. Hua, S. Li, R. Luo et al., “Icariin protects human nucleus pulposus cells from hydrogen peroxide-induced mitochondria-mediated apoptosis by activating nuclear factor erythroid 2-related factor 2,” Biochimica et Biophysica Acta - Molecular Basis of Disease, vol. 1866, p. 165575, 2020.

P. Sampara, R. R. Banala, S. K. Vemuri, G. R. AV, and S. GPV, “Understanding the molecular biology of intervertebral disc degeneration and potential gene therapy strategies for regeneration: a review,” Gene Therapy, vol. 25, no. 2, pp. 67–82, 2018.

G. Zheng, Z. Pan, Y. Zhan et al., “TFEB protects nucleus pulposus cells against apoptosis and senescence via restoring autophagic flux,” Osteoarthritis and Cartilage, vol. 27, no. 2, pp. 347–357, 2019.

W. Fang, X. Zhou, J. Wang et al., “Wogonin mitigates intervertebral disc degeneration through the Nrf2/ARE and MAPK signaling pathways,” International Immunopharmacology, vol. 65, pp. 539–549, 2018.

Y. Zhang, F. He, Z. Chen et al., “Melatonin modulates IL-1β-induced extracellular matrix remodeling in human nucleus pulposus cells and attenuates rat intervertebral disc degeneration and inflammation,” Aging (Albany NY), vol. 11, no. 22, pp. 10499–10512, 2019.

Z. Buser, A. S. Chung, A. Abedi, and J. C. Wang, “The future of disc surgery and regeneration,” International Orthopaedics, vol. 43, no. 4, pp. 995–1002, 2019.

S. Chen, M. Luo, H. Kou, G. Shang, Y. Ji, and H. Liu, “A review of gene therapy delivery systems for intervertebral disc degeneration,” Current Pharmaceutical Biotechnology, vol. 21, no. 3, pp. 194–205, 2020.

Z. Li, X. Li, C. Chen et al., “Long non-coding RNAs in nucleus pulposus cell function and intervertebral disc degeneration,” Cell Proliferation, vol. 51, no. 5, article e12483, 2018.

M. L. Ji, H. Jiang, X. J. Zhang et al., “Preclinical development of a microRNA-based therapy for intervertebral disc degeneration,” Nature Communications, vol. 9, no. 1, p. 5051, 2018.

X. Cheng, L. Zhang, K. Zhang et al., “Circular RNA VMA21 protects against intervertebral disc degeneration through targeting miR-200c and X linked inhibitor-of-apoptosis protein,” Annals of the Rheumatic Diseases, vol. 77, no. 5, pp. 770–779, 2018.

Z. Li, X. Chen, D. Xu, S. Li, M. T. Chan, and W. K. Wu, “Circular RNAs in nucleus pulposus cell function and intervertebral disc degeneration,” Cell Proliferation, vol. 52, article e12704, 2019.

J. Zhu, X. Zhang, W. Gao, H. Hu, X. Wang, and D. Hao, “lncRNA/circRNA-miRNA-mRNA ceRNA network in lumbar intervertebral disc degeneration,” Molecular Medicine Reports, vol. 20, no. 4, pp. 3160–3174, 2019.

L. Xie, W. Huang, Z. Fang et al., “CircERCC2 ameliorated intervertebral disc degeneration by regulating mitophagy and apoptosis through miR-182-5p/SIRT1 axis,” Cell Death & Disease, vol. 10, no. 10, p. 751, 2019.

S. Zhan, K. Wang, Y. Song et al., “Long non-coding RNA HOTAIR modulates intervertebral disc degenerative changes via Wnt/β-catenin pathway,” Arthritis Research & Therapy, vol. 21, no. 1, p. 201, 2019.

L. Yu, Y. Hao, C. Xu, G. Zhu, and Y. Cai, “LINC00969 promotes the degeneration of intervertebral disk by sponging miR-335-3p and regulating NLRP3 inflammasome activation,” IUBMB Life, vol. 71, pp. 611–618, 2018.

P. H. Lan, Z. H. Liu, Y. J. Pei et al., “Landscape of RNAs in human lumbar disc degeneration,” Oncotarget, vol. 7, no. 39, pp. 63166–63176, 2016.

M. E. Ritchie, B. Phipson, D. I. Wu et al., “limma powers differential expression analyses for RNA-sequencing and microarray studies,” Nucleic acids research, vol. 43, article e47, 2015.

P. Langfelder and S. Horvath, “WGCNA: an R package for weighted correlation network analysis,” BMC Bioinformatics, vol. 9, no. 1, p. 559, 2008.

J. H. Li, S. Liu, H. Zhou, L. H. Qu, and J. H. Yang, “StarBase v2.0: decoding miRNA-ceRNA, miRNA-ncRNA and protein-RNA interaction networks from large-scale CLIP-Seq data,” Nucleic Acids Research, vol. 42, pp. D92–D97, 2013.

P. Glažar, P. Papavasileiou, and N. Rajewsky, “circBase: a database for circular RNAs,” RNA, vol. 20, no. 11, pp. 1666–1670, 2014.

N. T. Doncheva, J. H. Morris, J. Gorodkin, and L. J. Jensen, “Cytoscape StringApp: network analysis and visualization of proteomics data,” Journal of Proteome Research, vol. 18, pp. 623–632, 2018.

D. Szklarczyk, A. L. Gable, D. Lyon et al., “STRING v11: protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets,” Nucleic Acids Research, vol. 47, no. D1, pp. D607–d613, 2019.

G. Dennis Jr., B. T. Sherman, D. A. Hosack et al., “DAVID: Database for Annotation, Visualization, and Integrated Discovery,” Genome Biology, vol. 4, no. 5, p. P3, 2003.

J. Long, X. Wang, X. Du et al., “JAG2/Notch2 inhibits intervertebral disc degeneration by modulating cell proliferation, apoptosis, and extracellular matrix,” Arthritis Research & Therapy, vol. 21, no. 1, p. 213, 2019.

L. Y. Jin, Z. D. Lv, K. Wang et al., “Estradiol alleviates intervertebral disc degeneration through modulating the antioxidant enzymes and inhibiting autophagy in the model of menopause rats,” Oxidative Medicine and Cellular Longevity, vol. 2018, Article ID 7890291, 12 pages, 2018.

J. Deng, W. Kong, X. Mou, S. Wang, and W. Zeng, “Identifying novel candidate biomarkers of RCC based on WGCNA analysis,” Personalized medicine, vol. 15, no. 5, pp. 381–394, 2018.

M. Ding, F. Li, B. Wang, G. Chi, and H. Liu, “A comprehensive analysis of WGCNA and serum metabolomics manifests the lung cancer-associated disordered glucose metabolism,” Journal of Cellular Biochemistry, vol. 120, no. 6, pp. 10855–10863, 2019.

L. Yin, Z. Cai, B. Zhu, and C. Xu, “Identification of key pathways and genes in the dynamic progression of HCC based on WGCNA,” Genes (Basel), vol. 9, no. 2, p. 92, 2018.

M. Niemira, F. Collin, A. Szalkowska et al., “Molecular signature of subtypes of non-small-cell lung cancer by large-scale transcriptional profiling: identification of key modules and genes by weighted gene co-expression network analysis (WGCNA),” Cancers, vol. 12, 2020.

J. J. Chan and Y. Tay, “Noncoding RNA:RNA regulatory networks in cancer,” International Journal of Molecular Sciences, vol. 19, no. 5, p. 1310, 2018.

Z. I. Johnson, Z. R. Schoepflin, H. Choi, I. M. Shapiro, and M. V. Risbud, “Disc in flames: roles of TNF-α and IL-1β in intervertebral disc degeneration,” European Cells & Materials, vol. 30, pp. 104–117, 2015.

C. Wang, X. Yu, Y. Yan et al., “Tumor necrosis factor-α: a key contributor to intervertebral disc degeneration,” Acta biochimica et biophysica Sinica, vol. 49, pp. 1–13, 2016.

S. Chen, S. Liu, K. Ma, L. Zhao, H. Lin, and Z. Shao, “TGF-β signaling in intervertebral disc health and disease,” Osteoarthritis and Cartilage, vol. 27, no. 8, pp. 1109–1117, 2019.

Y. Huang, Y. Wang, C. Wu, and W. Tian, “Elevated expression of hypoxia-inducible factor-2alpha regulated catabolic factors during intervertebral disc degeneration,” Life Sciences, vol. 232, p. 116565, 2019.