CircRNA expression profile and potential role of hsa_circ_0040039 in intervertebral disc degeneration.

PURPOSE: Circular RNAs (circRNAs) play an critical role in the pathological processes associated with IDD. However, the potential roles of circRNAs in IDD remain largely unclear. Here, we identify the circRNAs expression profiles and elucidate the potential role of candidate circRNAs in the pathogenesis of intervertebral disc degeneration (IDD) through microarray data and bioinformatics analyses.
METHODS: We obtained the datasets of microarrays (GSE67566 and GSE116726) from the Gene Expression Omnibus database. The differentially expressed circRNAs and miRNAs were identified using the Limma R package. The target miRNAs and target genes of the candidate circRNAs were predicted using an online tool. Functional enrichment analyses of the target genes were performed using the clusterProfiler R package. A protein-protein interaction (PPI) network was constructed using STRING.
RESULTS: A total of 104 differentially expressed circRNAs were identified between the IDD and the control groups, including 41 upregulated circRNAs and 63 downregulated circRNAs (cutoff criteria (|log2 fold change| > 2, P < .05)). Hsa_circ_0040039, which was the most upregulated circRNA (log2 fold change = 2.95), was selected for further analysis. The regulatory circRNA-miRNA-mRNA network comprised hsa_circ_0040039, 2 target miRNAs (hsa-miR-424-5p and hsa-miR-15b-5p), and 77 target genes. Functional enrichment analysis showed that the 77 promising target genes are mainly enriched in the ubiquitin proteasome system and Wnt signaling pathway. Further, the PPI network showed that the top 3 hub genes are BRTC, SIAH1, and UBE2V1.
CONCLUSIONS: A total of 104 differentially expressed circRNAs were identified between the IDD and control groups. Hsa_circ_0040039 may serve as a sponge of hsa-miR-424-5p and hsa-miR-15b-5p, to regulate the expression of downstream genes (such as BRTC, SIAH1, and UBE2V1); thus, it may be involved in IDD-associated pathological processes via the Wnt/β-catenin signaling pathway. Further studies are required to confirm the potential roles of hsa_circ_0040039 in IDD.

1. Introduction

Intervertebral disc degeneration (IDD) is a physiological and pathological process of natural aging and degeneration of lumbar spine. IDD is considered as one of the major causes of low back pain and spine-related disorders, and it severely affects the health of human beings and leads to tremendous economic burden to the whole society.[ The etiology of IDD is complicated; IDD can be attributed to multiple factors, including aging, environmental factors, and genetic factors.[ Of these, genetic factors are considered to be the dominant contributing factors to IDD.[ As documented previously, nucleus pulpous (NP) cells are the main type of cells residing in the NP; they play a vital role in maintaining disc health by synthesizing type II collagen, aggrecan, and other components of the extracellular matrix (ECM).[ Furthermore, the dysfunction of NP cells, including aberrant cell differentiation, proliferation, apoptosis, and ECM synthesis, markedly promote the occurrence and progression of IDD.[ Hence, it is meaningful to obtain deeper insights into the molecular mechanisms associated with IDD in NP cells.

Circular RNAs (circRNAs) are a special novel type of endogenous noncoding RNAs produced by the back-splicing of precursor mRNAs in eukaryotes; circRNAs are characterized by covalently closed loop structures without 5′ to 3′ polarity and polyadenylated tails.[ With the development of high-throughput sequencing and bioinformatics technologies, it has been revealed that circRNAs are widespread, conserved, and stable in mammalian cells, and are involved in transcriptional or posttranscriptional regulation.[ Specifically, circRNAs can function as miRNA sponges to influence the expression of target mRNAs, and they may be involved in the pathological processes associated with many diseases, including IDD.[ Cheng et al[ revealed that circVMA21 acted as a sponge of miR-200c, and alleviated NP cell apoptosis and ECM imbalance by targeting miR-200c and X-linked inhibitor-of-apoptosis protein (XIAP). Song et al also revealed that circRNA_104670 is upregulated in human IDD tissues and regulates MMP-2 expression by directly sponging miRNA-17-3p.[ However, the potential roles of circRNAs in IDD remain largely unclear. Thus, in this preliminary study, we attempted to identify the circRNA expression profiles and elucidate the potential role of candidate circRNAs in the pathogenesis of IDD using microarray data (GSE67566 and GSE116726) and bioinformatics analyses. We hope that this study can provide a feasible avenue for research regarding IDD-associated circRNAs.

2. Methods

2.1. Dataset selection and data processing

The circRNA expression profiling dataset (GSE67566) and the miRNA expression profiling dataset (GSE116726) were obtained from the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) for the study regarding the potential genetic factors that may contribute to human lumbar IDD. Briefly, 5 normal NP tissue samples and 5 IDD NP tissue samples were enrolled in the GSE67566 dataset (platform: GPL19978 Agilent-069978 Arraystar Human CircRNA microarray) and 3 normal NP tissue samples and 3 IDD NP tissue samples were enrolled in the GSE116726 dataset (platform: GPL20712 Agilent-070156 Human miRNA). The NP tissues were grouped by T2-weighted magnetic resonance images according to the Pfirrmann classification.[

After the datasets were downloaded, the probe IDs were converted into gene symbols based on matched platform annotation files using ActivePerl (Version 5.24.3.2404) (https://www.activestate.com/). Subsequently, we used Impute package (http://www.bioconductor.org/packages/) to impute missing the expression data in R (Version 3.5.3); then, we used the Limma package to the normalize the expression data and to identify differentially expressed circRNAs or miRNAs from the 2 datasets between patients with or without IDD in R (Version 3.5.3). The following cutoff criteria were considered for screening circRNAs or miRNAs: |log2 Fold Change| > 2 and the P value <.05. Further, a Volcano Plot was constructed to visualize the significant circRNAs, and a heatmap was constructed to visualize the significant circRNAs and perform a hierarchical clustering analysis with the aid of the Pheatmap package (https://cran.r-project.org/web/packages/) in R (Version 3.5.3).

2.2. Reconstruction of circrna-mirna-mrna regulatory network

The most upregulated circRNA was selected for further analysis, and its sequence was obtained from circBase (http://www.circbase.org/). After predicting the target miRNAs of the circRNA using RegRNA 2.0 (http://regrna2.mbc.nctu.edu.tw/) and starBase 3.0 (http://starbase.sysu.edu.cn/), the overlapping miRNAs were acquired based on the predicted target miRNAs and downregulated miRNAs. The Venn diagram was constructed using a web tool (http://bioinformatics.psb.ugent.be/webtools/Venn/). The distribution of differentially expressed circRNAs in human chromosomes was analyzed by GraphPad Prism 7 (https://www.graphpad.com/). Furthermore, the miRWalk web tool (http://mirwalk.umm.uni-heidelberg.de/), which comprises various algorithms, including Targetscan, miRDB, and miRTarBase, was used to predict the potential target mRNAs of the common miRNAs. The circRNA-miRNA-mRNA regulatory network was constructed to visualize the interrelationships between them using Cytoscape 3.7.1 (https://cytoscape.org/).

2.3. Functional annotation

The clusterProfiler package (http://www.bioconductor.org/packages/) was used to perform the functional and pathway enrichment analysis of the target mRNAs using R (Version 3.5.3), including the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses. These analyses provide a comprehensive set of functional annotation tools for the visualization and further understand the biological significance of large lists of target mRNAs.

2.4. Protein-protein interaction network construction

The functional protein-protein interaction (PPI) analysis is essential to interpret the molecular mechanisms underlying the key cellular activities in IDD. The STRING online tool (https://string-db.org) was used to analyze and visualize the interactions between the target mRNAs. An interaction score of 0.7 (high confidence) was regarded as the cutoff criterion, and disconnected nodes in the network were hidden.

3. Results

3.1. Identification of differentially expressed circrnas and mirnas

Based on the cutoff criteria (|log2 Fold Change| > 2, P < .05), a total of 104 circRNAs, including 41 upregulated circRNAs and 63 downregulated circRNAs, were differentially expressed in the IDD group compared with the control group. The top 5 upregulated and the top 5 downregulated circRNAs are ranked in Table 1 based on the fold changes in their expression levels. The hierarchical clustering (Fig. 1A) and the Volcano Plot (Fig. 1B) demonstrated the differential expression profiles of circRNAs in the GSE67566 dataset. The distribution of differentially expressed circRNAs in human chromosomes is presented in Figure 1C. In addition, a total of 840 differentially expressed miRNAs, including 485 upregulated miRNAs and 355 downregulated miRNAs, were identified between the IDD group and the control group (Table 1, Supplemental Digital Content, http://links.lww.com/MD/H77).

Table 1

The top 5 upregulated and downregulated circRNAs ranked by fold changes.

circRNA	log2FC*	P	Position	Strand	Sequence length	Best transcript	Gene symbol	Regulation
hsa_circ_0040039	2.953	<.00001	chr16:69279504-69318147	+	765	NM_006750	SNTB2	Up
hsa_circ_0092342	2.808	<.00001	chr11:8706439-8707219	+	780	NM_000990	RPL27A	Up
hsa_circ_0004354	2.799	<.00001	chr16:69317950-69318147	+	197	NM_006750	SNTB2	Up
hsa_circ_0028173	2.795	<.00001	chr12:110764194-110780253	+	1774	NM_170665	ATP2A2	Up
hsa_circ_0047151	2.760	<.00001	chr18:21001330-21010641	–	125	NM_032933	TMEM241	Up
hsa_circ_0082686	–3.319	<.00001	chr7:139741443-139746807	–	320	NM_022750	PARP12	Down
hsa_circ_0036763	–3.318	<.00001	chr15:90760670-90763123	+	326	NM_020210	SEMA4B	Down
hsa_circ_0003258	–3.179	<.00001	chr17:47388673-47389404	–	261	NM_001145365	ZNF652	Down
hsa_circ_0003183	–3.089	<.00001	chr8:61496766-61504528	+	288	NM_002865	RAB2A	Down
hsa_circ_0072464	–3.086	<.00001	chr5:53409031-53450469	–	269	NM_019087	ARL15	Down

Figure 1.

The expression profiles of circRNAs and their parental gene distribution in chromosomes. (A) Hierarchical clustering analysis showed 104 differentially expressed circRNAs (|log2 fold change| > 2, P < .05) between the IDD group and the control group. The red color represents the upregulated circRNAs and the green color represents the downregulated circRNAs. (B) Volcano plot of differentially expressed circRNAs (|log2 fold change| > 2, P < .05). The red plots represent the upregulated circRNAs with log2 (fold change) > 2, P < .05. The green plots represent the downregulated circRNAs with log2 (fold change) < –2, P < .05. logFC: log2 fold change. (C) The distribution of the parental genes of differentially expressed circRNAs in human chromosomes. IDD = intervertebral disc degeneration.

3.2. The circrna-mirna-mrna regulatory network

Acting as miRNA sponges, circRNAs may indirectly upregulate target mRNAs by inhibiting the negative regulation of miRNAs. Hsa_circ_0040039, which was the most upregulated circRNA (log2 Fold Change = 2.95), was selected for further analysis in our study. First, we obtained 4 overlapping miRNAs from the predicted target miRNAs of hsa_circ_0040039 (accessed using RegRNA 2.0 and starBase 3.0) and the 355 downregulated miRNAs (Fig. 2A), including hsa-miR-424-5p, hsa-miR-15b-5p, hsa-miR-5010-5p, and hsa-miR-574-5p. Then, a total of 77 target mRNAs of the 4 overlapping miRNAs were identified using the miRWalk online tool; these were predicted based on hsa-miR-424-5p and hsa-miR-15b-5p, rather than hsa-miR-5010-5p or hsa-miR-574-5p. The potential binding sites of hsa_circ_0040039 on hsa-miR-15b-5p and hsa-miR-424-5p are illustrated in Figures 2B, C, respectively. Finally, using the Cytoscape software, the hsa_circ_0040039 regulatory network was constructed using the 2 miRNAs and 77 mRNAs (Fig. 3).

Figure 2.

Venn diagrams of overlapping miRNAs and potential binding sites of hsa_circ_0040039 on hsa-miR-424-5p and hsa-miR-15b-5p. (A) The Venn diagrams show that 4 overlapping miRNAs were identified from the downregulated miRNAs and predicted target miRNAs (predicted from RegRNA 2.0 and starBase 3.0) of hsa_circ_0040039. (B) Three potential binding sites common to hsa_circ_0040039 and has-miR-15b-5p (obtained from starBase 3.0). (C) Three potential binding sites common to hsa_circ_0040039 and has-miR-424-5p (obtained from starBase 3.0).

Figure 3.

circRNA-miRNA-mRNA regulatory network. This network is composed of 1 circRNA (yellow octagon), 2 target miRNAs (purple tetragons), and 77 indirect target mRNAs (cyan nodes).

3.3. Enrichment analysis of the target mRNAs

Based on network described above, we performed GO and KEGG analyses of all the target mRNAs to provide potential guidance for the further functional verification of the roles of hsa_circ_0040039 and the processes it is associated with. In the GO analysis, the significantly enriched and meaningful terms were molecular function, including binding and activity. The top 15 significantly enriched terms of biological process (in the GO analysis) are shown in Figure 4A and Table 2, and the 7 significantly enriched pathways (in the KEGG analysis) are shown in Figure 4B and Table 3.

Figure 4.

GO and KEGG analyses of 77 target mRNAs. The length of the column represents the number of mRNAs. The color of the column represents the adjusted P value. Red: high-level enrichment; blue: low-level enrichment. (A) The top 15 significantly enriched GO terms. (B) The 7 significantly enriched KEGG signaling pathways. GO = gene ontology; KEGG = Kyoto Encyclopedia of Genes and Genomes.

Table 2

The top 15 enriched GO terms of target mRNAs.

ID	Term	Gene symbol	Count	P
GO:0046332	SMAD binding	ACVR2A, AXIN2, FKBP1A, SKI, SMAD3, TGFBR3	6	9.64E-07
GO:0019199	Transmembrane receptor protein kinase activity	ACVR2A, FGFR1, INSR, KDR, TGFBR3	5	2.12E-05
GO:0008013	Beta-catenin binding	AXIN2, BTRC, CARM1, SMAD3, TBL1XR1	5	2.25E-05
GO:0031625	Ubiquitin protein ligase binding	ATF6, AXIN2, BAG4, CUL2, GABARAPL1, PRKAR2A, SKI, SMAD3	8	2.34E-05
GO:0048185	Activin binding	ACVR2A, FKBP1A, TGFBR3	3	2.40E-05
GO:0044389	Ubiquitin-like protein ligase binding	ATF6, AXIN2, BAG4, CUL2, GABARAPL1, PRKAR2A, SKI, SMAD3	8	3.37E-05
GO:0005160	Transforming growth factor beta receptor binding	ACVR2A, FKBP1A, SMAD3, TGFBR3	4	5.83E-05
GO:0004713	Protein tyrosine kinase activity	FGF2, FGFR1, INSR, KDR, MAP2K3, WEE1	6	9.17E-05
GO:0004674	Protein serine/threonine kinase activity	ACVR2A, AKT3, CASK, CHEK1, MAP2K3, STK38, TAOK1, TGFBR3, WNK3	9	.000105
GO:0061630	Ubiquitin protein ligase activity	BTRC, COP1, CUL2, SIAH1, ZNRF2, ZNRF3	6	.000221
GO:0061659	ubiquitin-like protein ligase activity	BTRC, COP1, CUL2, SIAH1, ZNRF2, ZNRF3	6	.000251
GO:0019838	Growth factor binding	ACVR2A, FGFR1, INSR, KDR, TGFBR3	5	.000259
GO:0034713	Type I transforming growth factor beta receptor binding	ACVR2A, FABP1A	2	.000908
GO:0004842	Ubiquitin protein transferase activity	BTRC, COP1, CUL2, SIAH1, UBE2V1, ZNRF2, ZNRF3	7	.001246
GO:0005024	Transforming growth factor beta-activated receptor activity	ACVR2A, TGFBR3	2	.001714

Table 3

The 7 enriched KEGG signaling pathways of target mRNAs.

Pathway ID	Pathway name	Gene symbol	Count	P
hsa04550	Signaling pathways regulating pluripotency of stem cells	ACVR2A, AKT3, AXIN2, FGF2, FGFR1, SMAD3	6	8.61E-05
hsa04310	Wnt signaling pathway	AXIN2, BTRC, SIAH1, SMAD3, TBL1XR1, ZNRF3	6	.000175
hsa04115	p53 signaling pathway	CHEK1, COP1, SIAH1, ZMAT3	4	.000560
hsa04015	Rap1 signaling pathway	AKT3, FGF2, FGFR1, INSR, KDR, MAP2K3	6	.000725
hsa04010	MAPK signaling pathway	AKT3, FGF2, FGFR1, INSR, KDR, MAP2K3, TAOK1	7	.000859
hsa04014	Ras signaling pathway	AKT3, FGF2, FGFR1, INSR, KDR, SHOC2	6	.001346
hsa04218	Cellular senescence	AKT3, BTRC, CHEK1, MAP2K3, SMAD3	5	.001531

3.4. PPI network

PPI networks were constructed using the STRING online tool, with 35 target genes containing 35 nodes and 39 edges (Fig. 5A). Disconnected nodes were hidden in the network. The number of connected nodes of the top 15 genes are shown using a bar graph (Fig. 5B). The top 3 hub genes included beta-transducin repeat containing E3 ubiquitin protein ligase (BRTC), siah E3 ubiquitin protein ligase 1 (SIAH1), and ubiquitin conjugating enzyme E2 V1 (UBE2V1).

Figure 5.

PPI network analysis of the target mRNAs. An interaction score of 0.7 (high confidence) was regarded as the cutoff criterion, and the disconnected nodes in the network were hidden. (A) PPI network of target mRNAs. (B) The number of connected nodes for the top 15 genes. PPI = protein-protein interaction.

4. Discussion

CircRNAs are a group of widespread and conserved endogenous noncoding RNAs functioning as miRNA sponges and regulating the posttranscriptional activity of various genes in mammalian cells.[ It has been documented that circRNAs are involved in the pathological processes of many diseases, including hepatocellular carcinoma, heart disease, osteoarthritis, and IDD.[ Song et al reported that circRNA_104670 plays a vital role in IDD by functioning as a competing endogenous RNA (ceRNA).[ Cheng et al reported that the circRNA VMA21 protects against IDD by targeting miR-200c and XIAP.[ However, the roles of circRNAs and the mechanisms underlying their activity in the pathogenesis of IDD remain largely unclear. To explore the roles of circRNAs in human IDD-associated pathological processes, we identified 104 circRNAs that were differentially expressed between IDD NP and normal NP specimens from GSE67566 dataset. Among them, the most upregulated circRNA, hsa_circ_0040039, was selected for further analysis, and then, circRNA-miRNA-mRNA network construction, GO analysis, KEGG analysis, and PPI network construction were conducted. Our results suggested that hsa_circ_0040039 probably acted as a sponge of hsa-miR-424-5p and hsa-miR-15b-5p and upregulated the downstream target genes; it may be involved in IDD-associated pathological processes via the signaling pathways enriched in the aforementioned analyses.

However, the question, “what biological processes and regulatory mechanisms are mainly involved in the roles of hsa_circ_0040039?” remains unanswered. Our PPI network indicated that the top 3 hub genes are BRTC, SIAH1, and UBE2V1, which are mainly enriched in the ubiquitin proteasome system; further, in the GO analysis, processes such as β-catenin binding, ubiquitin protein ligase activity, ubiquitin-like protein ligase activity, and ubiquitin protein transferase activity were enriched (Table 2). The results of the KEGG analysis indicated that various pathways, including the Wnt signaling pathway, the p53 signaling pathway, and cellular senescence, especially that associated with the Wnt signaling pathway, were significantly enriched (Table 3). The Wnt signaling pathway is extremely important for developmental and homeostatic processes, including cell proliferation and differentiation and tissue formation.[ The Wnt signaling pathway comprises 2 different intracellular signaling pathways: the β-catenin-dependent signaling pathway (frequently referred to as the canonical Wnt/β-catenin signaling pathway) and the β-catenin-independent signaling pathway.[ Abnormal activation of the canonical Wnt/β-catenin signaling pathway is linked to many types of diseases, including pigmentary disorders, chronic obstructive pulmonary diseases, neurodegenerative diseases, and bone diseases.[ Furthermore, the expression levels of pivotal molecules involved in the canonical Wnt/β-catenin signaling pathway, such as β-catenin, are closely regulated by the ubiquitin proteasome system.[ Previous literature has documented that BRTC binds to a highly conserved consensus substrate recognition motif in β-catenin; this regulates the ubiquitylation and subsequent proteasomal degradation of β-catenin in a phosphorylation-dependent manner.[ On the contrary, SIAH1 represses the Wnt/β-catenin signaling pathway by targeting β-catenin for ubiquitylation and subsequent proteasomal degradation through a phosphorylation-independent manner.[ Some studies have documented that the Wnt/β-catenin signaling pathway involved in the occurrence and development of IDD. Xie et al reported that AQP3 exerted protective effects against IDD through the inhibition of the Wnt/β-catenin signaling pathway.[ Nilsson et al also reported that IDD was associated with the underactivation of the Wnt/β-catenin signaling pathway; further, the concomitant increase in β-catenin expression may lead to a regenerative response.[ In addition, Sun et al reported that miR-532 downregulated the Wnt/β-catenin signaling pathway via targeting Bcl-9 and induced the apoptosis of human intervertebral disc NP cells.[ To the best of our knowledge, the roles of circRNAs associated with the Wnt/β-catenin signaling pathway in NP cells under conditions of IDD are not yet well understood. The bioinformatics analyses performed in this study, together with the findings of previous reports, suggest that hsa_circ_0040039 may be a promising biomarker for IDD and may be involved in IDD-associated pathological processes through the Wnt/β-catenin signaling pathway. Future research is needed to validate the findings of this preliminary study and to elucidate the molecular mechanisms underlying the role of hsa_circ_0040039 in IDD.

5. Conclusions

A total of 104 differentially expressed circRNAs were identified between IDD NP tissues and normal NP tissues. Among these differentially expressed circRNAs, hsa_circ_0040039 may serve as a sponge of hsa-miR-424-5p and hsa-miR-15b-5p to regulate downstream genes (such as BRTC, SIAH1, and UBE2V1); it may be involved in IDD-associated pathological processes via the Wnt/β-catenin signaling pathway. Further studies are required to confirm the potential roles of hsa_circ_0040039 in IDD.

Author contributions

Conceptualization: Sheng Huang.

Data curation: Junlong Zhong, Qihua Qi.

Formal analysis: Qihua Qi, Gangan Liu.

Funding acquisition: Sheng Huang.

Investigation: Junlong Zhong, Qihua Qi.

Methodology: Junlong Zhong, Qihua Qi.

Project administration: Sheng Huang, Junlong Zhong.

Resources: Qihua Qi, Gangan Liu.

Software: Qihua Qi, Gangan Liu.

Supervision: Junlong Zhong, Ming Gong.

Validation: Sheng Huang, Junlong Zhong.

Visualization: Junlong Zhong, Ming Gong.

Writing – original draft: Junlong Zhong.

Writing – review & editing: Junlong Zhong, Ming Gong.

All authors read and approved this final manuscript.