Comprehensive Profile Analysis of Differentially Expressed circRNAs in Glucose Deprivation-Induced Human Nucleus Pulposus Cell Degeneration

Nucleus pulposus (NP) is the core substance to maintain the homeostasis of intervertebral disc and stability of biomechanics. The insufficient supply of nutrition (especially glucose) is an important factor that leads to the degeneration of NP cells. circRNAs play an important role in the process of intervertebral disc degeneration (IDD) by regulating the functions of NP cells. However, glucose deprivation-related circRNAs and their functions in IDD have not been reported. In this study, the differentially expressed circRNAs in NP cells after 0, 6, 12, and 24 h of glucose deprivation culture were detected by a microarray assay. Besides, time series clustering analysis by STEM software obtained the differentially up- and downregulated circRNAs during glucose deficiency. Then, the main functions and pathways of up- and downregulated circRNAs were predicted by the functional enrichment analysis. By constructing the circRNA-miRNA regulatory network, the potential mechanisms of the most differentially expressed circRNAs were predicted. In addition, according to in vitro validation, circ_0075062 was upregulated in degenerating NP tissues and glucose deprivation-induced NP cell degeneration. Based on Sanger sequencing and RNase tolerance assay, circ_0075062 was the circular transcript. Interfering with circ_0075062 expression could potentially alleviate the imbalance of extracellular matrix (ECM) synthesis and degradation in the NP cells induced by glucose deprivation. Together, these findings help us gain a comprehensive understanding of the underlying mechanisms of IDD, and circ_0075062 may be a promising therapeutic target of IDD.


Introduction
Intervertebral disc degeneration (IDD) is the main contributor to low back pain, which causes heavy economic burdens on society and families, as well as negative impacts on the quality of patients' work and life [1]. Multiple factors (such as heredity, environment, age, mechanics, and nutrition) and inflammatory factors (TNF, IL-1) are related to IDD [2][3][4][5]. However, the specific mechanism is still unknown. At present, medicine and operation are common ways of treating degenerative disc disease. But the success rate of relieving long-term pain is rather low, which may also induce abrasions among adjacent segments [6]. Hence, it is of high clinical significance to gain a deep understanding of the pathogenesis and pathological process of IDD, so as to seek therapeutic methods based on pathogeny.
Nucleus pulposus (NP) cells play an important role in maintaining the integrity of the intervertebral disc by secreting the extracellular matrix [7]. Abnormal NP cell function, such as reduced proliferation, cluster formation, and cellular senescence, can accelerate the development of IDD [8,9]. The intervertebral disc is the largest nonvascular tissue in the human body. NP cells gain nutrition from the diffusion of capillaries in the endplate; thus, their nutrition supply is limited for long time [10]. Glucose is one of the nutritional ingredients for NP cells. The concentration gradient of glucose in adults' NP cells declines from 5.5 mM in endplates to 0.5 mM and even lower in the center of NP [11]. According to previous studies, the decreasing supply of nutrition, especially glucose, has a significant impact on the proliferation or survival of NP cells [12][13][14]. For this reason, further research and studies on the potential mechanism of NP cell dysfunction in the shortage of glucose may provide new strategies for preventing and treating IDD.
circRNAs are a new type of RNA with covalent and closed-loop structures and no 5 ′ cap or 3 ′ tail. Due to the cyclic structures, circRNAs are not prone to the degradation of RNA exonuclease and more stable than linear RNA. Moreover, circRNAs demonstrate unique traits of illnesses and developmental stages in different pathological contexts, indicating that circRNAs have more advantages in serving as disease diagnosis and therapeutic targets [15,16]. More and more evidence shows that circRNAs play important roles in IDD development [17][18][19]. However, glucose deprivationrelated circRNAs and the functions in IDD have not been reported yet.
In this study, by microarray analysis and time series clustering analysis, we obtained up-and downregulated circRNAs in NP cells under glucose deprivation for 0, 6, 12, and 24 h. Then, the role of differentially expressed circRNAs in glucose deprivation-induced NP cell degeneration is investigated by bioinformatics prediction and in vitro validation. Our results contribute to the understanding of the underlying mechanisms of IDD and establish a novel research target for IDD.

Ethical Statement.
This research program has been approved by the ethical review committee of Yangpu Hospital, Tongji University. We have acquired the informed consent form from all participants involved in the study.
2.2. Cell Culture and Treatment. Human NP cells were purchased from ScienCell Research Laboratories (Carlsbad, CA, USA), which had been extracted from nucleus pulposus. After being resuscitated, NP cells were cultured in DMEM/F12 with 10% fetal bovine serum and 100 mg/ml streptomycin at the temperature of 37°C and 5% CO 2 humidity. The culture medium was replaced once every three days. The secondgeneration cells were used in subsequent experiments. In the complete medium with no glucose, nucleus pulposus cells were incubated for 0 h, 6 h, 12 h, and 24 h. Total RNA was extracted from NP cells at each time node with three biological repetitions.

RNA Extraction.
Total RNA containing small RNA was extracted using the Trizol reagent (Invitrogen) and purified with the mirVana miRNA Isolation Kit (Ambion, Austin, TX, USA) according to the manufacturer's protocol. In addition, a spectrophotometer (NanoDrop ND-1000) was utilized to determine RNA purity and concentration from OD260/280 reading. RNA integrity was confirmed by 1% agarose gel electrophoresis.

circRNA Microarray
Imaging and Data Analysis. The extracted RNA was digested, dephosphorylated, transformed, extended, and marked with Cy3-dCTP according to the manufacturer's protocol. The marked RNA was mixed with the microarray (CapitalBio Technology Human Cir-cRNA Array v2) which also included circRNA probes from 170,304 persons. Then, the Agilent Scanner G2505C (Santa Clara, USA) was used in cleaning, fixating, and scanning the arrays to obtain images. Furthermore, GeneSpring software V13.0 (Agilent) summarized circRNA array data and analyzed standardization and quality control. In order to select differentially expressed circRNA, we adopted threshold values ≥ 2 and ≤−2. The Benjamini-Hochberg revised p value was 0.05. In addition, the Adjust Data function of the software CLUSTER 3.0 was used to conduct log 2 conversion of data and confirm the midvalue-based centered upon genes. Lastly, further analysis was carried out using hierarchical clustering with average linkage.

Time Series Clustering Analysis. Short Time-series
Expression Miner (STEM) software was used for data clustering analysis. This software program is specifically designed for analyzing short time series (3 to 8 temporal points) microarray gene expression data. The genes were allocated to a set of predefined model configuration files, so as to correctly detect temporal distribution with relevant associative functions [20,21].

Function and Pathway Enrichment Analysis. Gene
Ontology (GO) analysis [22] of predicting target genes and pathway analysis of Kyoto Encyclopedia of Genes and Genomes (KEGG) [23] were conducted through online tools. The database was used in annotation, visualization, and comprehensive discovery (DAVID, http://david.ncifcrf.gov/) [24]. GO analysis depicted our understanding of biology from three aspects: biological process (BP), cellular component (CC), and molecular function (MF). KEGG was the integrated database resource designed for explaining genomic sequence and other high-throughput data. p < 0:05 was regarded as the threshold value.

Construct a circRNA-miRNA Regulatory Network.
According to the latest research, circRNA could be used as the miRNA sponge to regulate gene expressions. In order to assess the potential function of circRNA, we used the miRanda algorithm to predict the target miRNA (miRanda score ≥ 140) of differentially expressed circRNA [25]. Target miRNA capable of combining more than three candidate circRNAs was identified as reliable interaction pairs of circRNA-miRNA. Cytoscape software was used to visualize the circRNA-miRNA regulatory network.  Normalized intensity  3 BioMed Research International hypertension or diabetes, etc. The gathered NP tissue was used in qRT-PCR and structure verification for the candidate circRNA.
2.9. Sanger Sequencing and RNase R Processing. NP tissue RNA was extracted. According to the circRNA circularization site, specific divergent primers were designed. Then, TA cloning and sequencing were carried out on the amplified products of divergent primers.
Total RNA of NP tissues was extracted and processed at 37°C with or without RNase R to obtain cDNA through inverse transcription of arbitrary primers. Then, cDNA was amplified by convergent primers and divergent primers.
Electrophoresis of the PCR product was carried out with 2% agarose gel electrophoresis.
2.10. Quantitative Real-Time PCR (qRT-PCR). The Trizol reagent was used to extract total RNA of NP cells and NP tissues. RNA concentration, purity, and integrity were detected by ND-1000 spectrophotometry and 1% agarose gel electrophoresis. Based on the instructions of PrimeScript™ RT Reagent Kit with Genomic DNA (gDNA) Eraser Kit, complementary DNA (cDNA) was synthesized through reverse transcription. qRT-PCR was carried out based on instructions of SYBR Premix Ex Taq II kits. 18s was the internal reference of circRNA. The relative expression level of data was  Table 1.
2.11. Western Blotting. Cells were transfected with the circ_ 0075062 interference plasmid using Lipofectamine 2000 (Invitrogen, Carlsbad, CA, USA), according to the manufacturer's protocol. Then, transfected cells were collected for western blotting analysis after 24 h of glucose deprivation incubation. Briefly, treated NP cells were completely lysed in ice-cold RIPA Lysis Buffer containing 1% Protease Inhibitor Cocktail (Sigma-Aldrich) and then centrifuged at 12000 r/min for 15 min. The bicinchoninic acid protein assay was performed to monitor the protein concentration. Then, proteins were separated in a 10% SDS-PAGE gel and transferred to PVDF membranes. The membranes were subsequently blocked with 3% bovine serum albumin/Tris-buffered saline containing Tween 20 for 2 hours and then incubated with anti-MMP-13 (Affinity, AF5355, 1 : 1000), anti-aggrecan (Affinity, DF7561, 1 : 1000), anti-collagen II (Affinity, AF0135, 1 : 1000), and anti-GAPDH (Sangon, D110016, 1 : 1000) at 4°C overnight and then incubated with a 1 : 5000 secondary antibody for 2 hours. The protein signals were detected by the enhanced chemiluminescence method. GAPDH was used as an endogenous normalization of protein. The intensities of bands were analyzed using ImageJ.
2.12. Statistical Analysis. Statistical analysis was conducted with GraphPad Prism 8.0. All data were represented by means ± SEM. Each group of the experiment was repeated three times. Student's t-test or one-way analysis of variance was used to analyze the significance of differences. p < 0:05 was regarded as having statistical significance.

circRNA Expression Profile Analysis in Glucose
Deprivation-Induced NP Cell Degeneration. In order to find out the potential role of circRNAs in glucose deprivationinduced NP cell degeneration, the NP cells were incubated in the complete medium with no glucose for 0 h, 6 h, 12 h, and 24 h. Then, high-throughput circRNA microarray assays were performed. The expression value of microarray data before pretreatment is shown in Figure 1(a). Standardized medians were on the same horizontal line with a favorable effect (Figure 1(b)). Principal component analysis (PCA) is used to visualize clusters in 12 samples (Figure s1). p ≤ 0:05 and |log 2 ðfold changeÞ | ≥2 were used as a significant cutoff criterion to conduct differential expression analyses with respect to each time point (0, 6, 12, and 24 h), and 1135 (12h_vs_6h), 7438 (24h_vs_12h), 12521 (24h_vs_6h), 2953 (6h_vs_0h), 3246 (12h_vs_0h), and 13578 (24h_vs_0h) differential expression circRNAs were obtained, respectively (Figures 1(c) and 1(d), Figure s2).

Time Series Clustering
Analysis. The STEM software was used to conduct time series cluster analysis and obtain 30     Files 1 and 2).

Functional Enrichment Analysis of DECs.
For host genes of the selected up-and downregulated circRNAs, GO and KEGG pathway enrichment analyses were performed to predict their potential functions. The top 30 GO and pathway items are shown in the bubble diagram ( Figure 3). We found that the gene symbol of downregulated circRNAs is mainly enriched in biological processes such as cytoskeleton organization and cell migration, as well as pathways such as cytokine-cytokine receptor interaction and metabolic pathways (Figures 3(a) and 3(b)). Functions of the gene symbol in upregulated circRNAs are enriched in biological processes like ion transmembrane transport and regulation of ion transport (Figure 3(c)) and pathways like focal adhesion, ECM-receptor interaction, and calcium signaling pathway (Figure 3(d)).

circRNA-miRNA Regulatory Network Construction.
Based on the miRanda algorithm, the target miRNAs of the top 5 up-and downregulated circRNAs ( Table 2) were predicted. As we all know, one type of miRNA could combine with one or more circRNAs. We selected miRNA that could combine with more than three circRNAs to construct the circRNA-miRNA regulatory network, as shown in Figure 4. circ_0075062 (degree, 12) and circ_0075063 (degree, 12) were two upregulated circRNAs with the highest connectivity. Besides, miR-5008-5p (degree, 5), miR-1915-3p (degree, 4), and miR-4726-3p (degree, 4) were the target miRNAs with high connectivity. These circRNAs and miRNAs may be the important regulator in glucose deprivation-induced NP cell degeneration.

Expression and Characterization of circ_0075062.
In vitro validation was performed on the selected key circRNAs. qRT-PCR results showed that circ_0075062 apparently increased in the NP cells under glucose deprivation condition and degenerative NP tissues, which were in accordance with microarray data (Figures 5(a) and 5(b)). Sanger sequencing was used to verify PCR amplified product, which indicated that circ_0075062 was head-to-tail splicing ( Figure 5(c)). circ_0075062 originated from the fourth exon of STC2 genes. Due to the special topological structures, circRNA was quite resistant to RNase R. After being digested by RNase R, the treatment group could still detect circ_0075062, but linear RNA could not be detected ( Figure 5(d)). The above results showed that circ_0075062 was a circular transcript. To evaluate the specificity of the circ_0075062 interference plasmid targeting the circular structure of circ_0075062, we determined the relative expression of circ_0075062 and linear 0075062 in NP cells ( Figure 5(e)). The results showed that the circ_0075062 interference plasmid downregulated the expression of circ_0075062, while the expression of linear 0075062 was not affected, indicating successful circ_ 0075062 interference. Downregulation of circ_0075062 promotes ECM synthesis (collagen II and aggrecan) and inhibits ECM degradation (MMP-13) under glucose-restricted condition. Collectively, these results suggest that circ_0075062 is a novel target and downregulation of circ_0075062 can inhibit glucose deprivation-induced NP cell degeneration.

Discussions
IDD has always been a health issue around the world, which has brought heavy burdens to the healthcare system and the economy [26]. NP, as an important part of intervertebral disc, plays a critical role in maintaining the microenvironmental stability and mechanical property of the intervertebral disc. A comprehensive understanding of the NP dysfunction mechanism will lay the foundation for the recovery, regeneration, and tissue engineering strategies of IDD-related disease [27].  I  I  I  II  I I  I  I  I  I  I   Recent studies showed that circRNAs can regulate the development and progression of IDD, including cell proliferation, apoptosis, ECM synthesis/degradation, and the generation of proinflammatory cytokine [28]. Inadequate supply of nutrition (especially glucose) is a primary obstacle in the regeneration of NP cells [12][13][14]. However, the expression and function of glucose deficiency-associated circRNAs in IDD remain unknown. In this study, we obtained circRNA expression profiles in response to glucose deficiency using microarray technology and time series clustering analysis in an experimental nutrient deprivation cell model. By using the tools of bioinformatics, biological functions and regulatory mechanisms of DECs were predicted. Further in vitro experiments suggest that circ_0075062 may be a new target for IDD intervention.
With the development of biotechnology, microarray analysis has become a powerful tool for characterizing many pathophysiological processes. Microarray technology has also identified many diagnostic markers and therapeutic targets in IDD [29,30]. Exposure of NP cells to serum deprivation [31], oxidative stress [32], and pressure response [33] is an effective in vitro model to study the molecular mechanism of IDD. In rat NP cells with serum starvation, microarray analysis indicated that significant differences in the genetic expressions of cell cycle DNA damage checkpoints [31]. In the compressed human NP cells, circRNA microarray analysis and qRT-PCR revealed that circRNA-CIDN was significantly downregulated and overexpression of circRNA-CIDN inhibited compression-induced apoptosis and NP ECM degradation [33]. Our previous studies found that glucose deprivation induced time-dependent morphological changes, proteoglycan degradation, and apoptosis in NP cells [34]. More and more research demonstrates that the analysis of microarray time series data plays an important role in bioinformatics studies and biomedical applications. The STEM software can obtain models with statistical significance from short time series microarray experiments. Besides, it can also compare the datasets among various experiments, presenting its analysis of the data in a highly visual and interactive manner [21]. Hence, our study can provide a more reliable and novel research target for IDD.
According to the GO analysis, downregulated circRNAs were mainly related to biological processes including cytoskeleton organization and cell migration. Upregulated cir-cRNAs were related to biological processes including ion transmembrane transport and regulation of ion transport. IVD involves three cytoskeleton elements (F-actin, β-microtubulin, and vimentin), which are of significance in cell division, movement, protein transport, and secretion. F-actin can regulate the mechanical conduction between ECM and IVD. Cytoskeleton network disorder may contribute to the loss of IVD stability and catabolism phenotype [35]. Recently, Zhang et al. found that SM/J early degeneration of the intervertebral disc was related to changes in genetic expressions of the ion transport system [36]. Transient receptor potential (TRP) channels, a superfamily of multimodal ion channels, were perhaps the potential trigger for intervertebral disc disease. Sadowska et al. confirmed the significance of TRP channels in IVD stability and pathology, which may consider them pharmacological targets for the treatment of IDD and low back pain [37]. These evidences indicate that DECs may be involved in the pathophysiological processes of IDD by regulating the above biological processes. KEGG pathway analysis indicated that downregulated circRNAs were involved in cytokine-cytokine receptor interaction and metabolic pathways, while upregulated circRNAs were involved in focal adhesion and ECM-receptor interaction signal channels, which are closely related to NP cell regeneration, extracellular matrix synthesis, and catabolism [5,38,39]. Nutrition deprivation could lead to NP cell apoptosis and reduce the expression of collagen II and aggrecan, which are two basic and important elements in the stability of the intervertebral disc [31,34]. Hence, we predicted that DECs may play important roles in IDD progression by regulating ECM metabolism. circRNA can act as a sponge for miRNA, thus influencing its activity and increasing expression of miRNA target genes [40]. Cui and Zhang showed that circ_001653 combines with miR-486-3p to regulate biological characteristics of NP cells and ECM synthesis [18]. Guo et al. reported that circ-FAM169A enhanced ECM catabolism through the miR-583/BTRC axis, thereby promoting IDD progression [19]. An increasing number of circRNAs regulating the progression of IDD have been discovered. However, circRNAmediated regulatory mechanisms are still poorly understood. In this study, by constructing the circRNA-miRNA regulatory network, we screened key circRNAs as well as important target miRNAs during glucose deprivation-induced NP cell degeneration. Further in vitro validation showed that circ_ 0075062 is involved in the ECM metabolic process of the NP cells under glucose deprivation condition, which is consistent with our bioinformatics predictions.

Conclusions
We explored the changes in the expression of circRNAs in NP cells during glucose deprivation. Through bioinformatics analysis and in vitro validation, circ_0075062 was identified as the key circRNA and predicted the potential mechanism of its regulation of IDD progression. This study enhanced the understanding of the molecular mechanism of IDD. The novel circ_0075062 may be the potential therapeutic target for IDD.