ETM Experimental and Therapeutic Medicine 1792-0981 1792-1015 D.A. Spandidos 10.3892/etm.2020.9390 ETM-0-0-09390 Articles Comprehensive evaluation of differential long non-coding RNA and gene expression in patients with cartilaginous endplate degeneration of cervical vertebra YuanJinghong 1 2 3 JiaJingyu 1 2 3 WuTianlong 1 2 3 LiuXijuan 4 HuShen 1 ZhangJian 1 DingRui 1 PangChongzhi 1 ChengXigao 1 2 3 Department of Orthopedics, The Second Affiliated Hospital of Nanchang University, Nanchang, Jiangxi 330006, P.R. China Institute of Orthopedics of Jiangxi Province, The Second Affiliated Hospital of Nanchang University, Nanchang, Jiangxi 330006, P.R. China Institute of Minimally Invasive Orthopedics of Nanchang University, The Second Affiliated Hospital of Nanchang University, Nanchang, Jiangxi 330006, P.R. China Department of Pediatrics, The Second Affiliated Hospital of Nanchang University, Nanchang, Jiangxi 330006, P.R. China Correspondence to: Professor Xigao Cheng, Department of Orthopedics, The Second Affiliated Hospital of Nanchang University, 1 Minde Road, Donghu, Nanchang, Jiangxi 330006, P.R. China 228206846@qq.com 12 2020 27 10 2020 20 6 260 13 04 2020 31 07 2020 Copyright: © Yuan et al. 2020 This is an open access article distributed under the terms of the Creative Commons Attribution-NonCommercial-NoDerivs License, which permits use and distribution in any medium, provided the original work is properly cited, the use is non-commercial and no modifications or adaptations are made.

Long non-coding RNAs (lncRNAs) are emerging as key regulators in gene expression; however, little is currently known regarding their role in cartilaginous endplate (CE) degeneration (CED) of cervical vertebra. The present study aimed to investigate the expression levels of lncRNAs and analyze their potential functions in CED of cervical vertebra in patients with cervical fracture and cervical spondylosis. Human competitive endogenous RNA (ceRNA) array was used to analyze lncRNA and mRNA expression levels in CE samples from patients with cervical fracture and cervical spondylosis, who received anterior cervical discectomy and fusion. Differentially expressed lncRNAs (DELs) or differentially expressed genes (DEGs) were identified and functionally analyzed, using Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses. An lncRNA-microRNA(miRNA)-mRNA ceRNA regulatory network was constructed based on the DELs and DEGs, and the ceRNA network was visualized using Cytoscape 3.7.2 software. In total, one downregulated mRNA, one upregulated miRNA and five downstream regulated lncRNAs were identified using reverse transcription-quantitative PCR in CED and healthy CE samples. A total of 369 lncRNAs and 246 mRNAs were identified as differentially expressed in CE. The GO and KEGG analyses demonstrated that the majority of GO and KEGG enrichments were associated with CED. Furthermore, a ceRNA network was established, including 168 putative miRNA response elements, 189 upregulated and 37 downregulated lncRNAs and 47 upregulated and 10dow regulated DEGs. The present study analyzed the function of DEGs in the ceRNA network and filtered out the same items as in DEG-function enrichment analysis. These results provide a new perspective for an improved understanding of ceRNA-mediated gene regulation in cervical spondylosis, and provide a novel theoretical basis for further studies on the function of lncRNA in cervical spondylosis. However, further experiments are required to validate the results of the present study.

 endplate cartilage long non-coding RNAs competitive endogenous RNA mRNAs bioinformatics Introduction

Lower back pain and neck pain are two common symptoms caused by intervertebral disc degeneration (IDD) (1). The global lower back and neck pain disability rates increased by 59.5% from 1990 to 2015, and will likely increase further with an aging population (2). Currently, therapies for cervical IDD (CIDD) are intended to relieve a multitude of symptoms, offering only temporary benefits rather than a permanent recovery (3). This is mainly a result of poor understanding of the exact etiology and pathogenesis. Xia et al (4) previously reported that transplantation of gelatin microspheres loaded with nucleus pulposus-like cells (NP-1Cs) and growth and differentiation factor-5 (GDF-5) into the intervertebral discs of the acupuncture rat tail vertebrae partially regenerated damaged intervertebral discs in a rat model. However, cartilage endplate (CE) degeneration (CED) is an obstacle for effective stem cell therapy (5). Therefore, it is necessary to determine the mechanisms underlying CIDD and CED to develop novel effective treatment methods for these diseases.

The outer layer of intervertebral discs is annulus fibrosus tissue, which surrounds the nucleus pulposus tissue and connects the superior and inferior vertebral bodies with CE. The supply of nutrition to intervertebral discs depends on the pores of the superior and inferior CE of the vertebral body (6,7). Previous studies have reported that CE serves a significant role in maintaining the basic function of the intervertebral disc (8). Therefore, CED is the primary factor leading to IDD, and maintaining the physiological function of CE is essential to prevent and treat IDD. Previous studies have identified factors associated with cartilage (9). For example, one study reported that cartilage degeneration is associated with the downregulation of Cbp/p300 interacting transactivator with Glu/Asp rich carboxy-terminal domain 2 (CITED2) expression (9).

A recent study revealed that long non-coding RNAs (lncRNAs), a subset of non-coding transcripts with a length of >200 nucleotides that have a low protein-coding potential, are involved in a series of biological processes such as glucose and lipid metabolism (10). Some lncRNAs can serve as competitive endogenous RNAs (ceRNAs) of microRNA (miRNA/miR) to modulate the downregulation of miRNA targets (11). The functions of lncRNAs in cancer and several other diseases have been investigated (12-15). For instance, Chen et al (12) reported that LINC00173.v1 upregulated the expression of VEGFA via sponging miR-511-5p to promote the proliferation and migration of vascular endothelial cells. Wang et al (13) also revealed that CARL lncRNA inhibited mitochondrial fission and apoptosis by upregulating prohibitin 2 expression via decreasing the expression of miR-539. However, to the best of our knowledge, the expression profiles and potential functions of lncRNAs in CED of cervical vertebra remain unknown.

In the present study, lncRNA and mRNA microarrays were used to identify differentially expressed lncRNAs (DELs) and differentially expressed genes (DEGs) between degenerate and healthy CE. Furthermore, a ceRNA regulatory network was constructed based on DELs and DEGs. The expression levels of selected DELs and DEGs were further assessed using reverse transcription-quantitative PCR (RT-qPCR). The present study aimed to identify novel lncRNAs and genes that are relevant to CED, which may further guide investigations and contribute to the development of novel therapeutic strategies for the treatment of CIDD.

 Materials and methods <sec> <title>Tissue samples

A total of 38 CE specimens from 19 patients with CIDD and 19 healthy subjects from February 2017 to February 2020 were enrolled in the study. The CIDD and healthy samples were collected from The Second Affiliated Hospital of Nanchang University (Nanchang, China). Details of the samples are presented in Table SI. Degenerated cervical vertebral CE specimens were collected from patients with CIDD that suffered from cervical spondylosis myelopathy, who had received anterior cervical discectomy and fusion (ACDF). CE specimens of healthy subjects were obtained from patients with a cervical fracture who received ACDF. None of the enrolled subjects had undergone radiotherapy or chemotherapy, and none had a history of surgery. Sample weight was evaluated before the study was conducted. Of the lesions collected, six legions weighing >50 mg were used for microarray detection. The remaining 32 lesions, which weighed <50 mg, were used for RT-qPCR and western blot analyses. The present study was approved by the Ethics Committee of The Second Affiliated Hospital of Nanchang University. Written informed consent was obtained from all participants. All tissue samples were preserved in liquid nitrogen until RNA and protein extraction.

 lncRNA and mRNA microarray analysis

Total RNAs from three degenerative CE and three healthy CE (HCE) tissues were extracted and purified using the RNeasy micro kit (Qiagen GmBH) and the RNase-Free DNase set (Qiagen GmBH). The RNA integrity coefficient number, as a measure of RNA integrity, was measured using an Agilent Bioanalyzer 2100 (Agilent Technologies, Inc.). The purity and quantity of initial total RNA chip samples were determined using a NanoDrop ND-2000 spectrophotometer (NanoDrop Technologies; Thermo Fisher Scientific, Inc.) and an Agilent Bioanalyzer 2100 to exclude genomic DNA contamination. RNA integrity was also determined on an ethidium bromide-stained 2% agarose gel. Qualified RNAs, of which the absorbance ratio at 260 and 280 nm is 2.0 was used for the following chip experiments.

 RNA amplification and labeling

Total RNA was amplified and labeled using the Low Input Quick Amp Labeling kit, One-Color (Agilent Technologies, Inc.) according to the manufacturer's instructions. Cy3-labeled cRNA was purified using the RNeasy mini kit according to the manufacturer's instructions (Qiagen GmBH).

 Microarray hybridization

The SBC Human (4x180 K) ceRNA array (Shanghai Shibei Biotechnology Co., Ltd.) was used to detect lncRNA and mRNA expression profiles. Each slide was hybridized with 1.65 µg Cy3-labeled cRNA using a Gene Expression Hybridization kit (Agilent Technologies, Inc.) in a Hybridization Oven (Agilent Technologies, Inc.) at 65˚C, according to the manufacturer's protocols. After 17 h of hybridization, the slides were washed in staining dishes (Thermo Fisher Scientific, Inc.) using a Gene Expression Wash Buffer kit (Agilent Technologies, Inc.), according to the manufacturer's protocol. All raw data have been uploaded to Gene Expression Omnibus (GSE153761).

 Data acquisition

Slides were scanned using an Agilent Microarray scanner (Agilent Technologies, Inc.) with the default settings: Dye channel=Green, Scan resolution=3 µm, photomultiplier tube 100% and 20 bit. Data were extracted using Feature Extraction 12.0 software (Agilent Technologies, Inc.). Raw data were normalized using the Quantile algorithm and 3.11 version of limma packages (https://www.bioconductor.org/packages/release/bioc/html/limma.html) in R 3.6.0 software (https://www.R-project.org/). Statistically significant DELs and DEGs between the two groups were defined as fold-change ≥2 and P<0.05. Furthermore, heat maps were constructed to present expression profiles of DELs and DEGs using hierarchical clustering, which was performed using the 1.0.12 version of pheatmap package (https://cran.r-project.org/web/packages/pheatmap/index.html) in R. In addition to these, the repeatability and accuracy of ceRNA array in the current study was also assessed using the Coefficient of Variation (CV). The value of CV was calculated according to the 10 repeated probe signals. If the CV is <15%, the outcomes of ceRNA array can be considered to be of excellent stability. Statistical analyses were performed using R. A two-tailed P<0.05 was considered indicate a statistically significant difference for all tests.

 Function and pathway analysis of DEGs

Gene Ontology (GO; www.geneontology.org) is a widely used ontology resource in the bioinformatics field. The present study analyzed the association of DEGs in the Biological Processes (BP), Molecular Function (MF) and Cellular Component (CC) sets in the GO database. The P-value denotes the significance of the enriched GO terms among DEGs.

Kyoto Encyclopedia of Genes and Genomes (KEGG; https://www.genome.jp/kegg/) pathway analysis is a systematic analysis tool and database of gene function and genome information, which helps researchers to study and express information in the context of entire gene networks (16). KEGG enrichment analysis of DEGs can identify differentially enriched pathways, which is useful for identifying biological regulation pathways that are significantly altered under experimental conditions (16). The enrichment analysis was performed using Fisher's exact test in the PANTHER classification system (http://www.pantherdb.org) and the 3.11 version of Cluster Profiler package from R/Bioconductor (https://www.bioconductor.org/) (17,18). The standard of selection was the number of genes attaining a P<0.05 threshold.

 Construction of ceRNA regulatory network

To determine the interaction between lncRNAs and mRNAs, data regarding lncRNAs and mRNAs were combined with miRNA data to construct a lncRNA-miRNA-mRNA ceRNA regulatory network. Using the miRanda tool (version 3.3a; http://www.microrna.org/), the present study predicted the lncRNA-miRNA regulatory combinations. miRNAs that had regulatory associations with mRNAs in the ceRNA network were queried in three databases, miRDB (version 6.0; http://www.mirdb.org), TargetScan (version 7.2; http://www.ta-rgetscan.org/vert_72/) and miRTarBase (version 8.0; http://mirtarbase.mbc.nctu.edu.t-w), which included experimentally validated miRNA-target interactions (19-21). The overlapping genes between targets of the identified miRNAs and the total of DEGs were retained for the ceRNA network construction, with associated DELs and miRNAs. Based on the identified miRNAs and miRNA targets in DELs and DEGs, a ceRNA regulatory network was constructed and visualized using Cytoscape software version 3.7.2(22).

 Functional analysis of genes in the ceRNA regulatory network

In the ceRNA regulatory network, a GO enrichment analysis was performed using Fisher's exact test in the PANTHER classification system to assess the potential functions of DEGs in the categories of BP, MF and CC. Over-represented enriched KEGG pathways were searched for using the 3.11 version of ClusterProfiler package of R/Bioconductor to identify pathways that were strongly associated with DEGs in the ceRNA network. P<0.05 was set as the threshold for significance.

 RT-qPCR

Total RNA was extracted from 32 samples (16 patients and 16 healthy control) of CE using TRIzol® reagent (Thermo Fisher Scientific, Inc.). Expression levels of one DEG [integrin subunit β8 (ITGB8)], miR-20a-5p and five DELs (ENST00000548900, lnc-MYBPC1-1:1, lnc-ARL13A-1:1, lnc-C2orf67-1:1 and lnc-DNAJB6-3:1) were detected using a PrimeScript™ RT reagent kit (Perfect Real-Time) and SYBR® Fast qPCR Mix (both from Takara Bio, Inc.). The temperature protocol for reverse transcription was: 16˚C for 30 min, 42˚C for 30 min and 85˚C for 5 min. cDNA was subjected to initial denaturation at 94˚C for 3 min, followed by 40 cycles at 94˚C for 10 sec and 60˚C for 40 sec, followed by extension at 72˚C for 10 min, using the specific primers. All experiments were repeated three times. β-actin was used as an internal reference for these lncRNAs and ITGB8 whereas U6 was used as an internal reference for miR-20a-5p. In the present study, 2-ΔΔCq method was used for relative quantification (23). Primer sequences are presented in Table SII.

 Western blotting

Western blot analyses were performed on the proteins extracted from 22 samples (11 patients and 11 healthy control) of CE tissue samples by using the RIPA lysis buffer (cat. no. P0013B; Beyotime Institute of Biotechnology). Protein concentrations were measured using a BCA protein assay reagent (Beyotime Institute of Biotechnology). Each lane was loaded with 15 µg protein. Samples were then separated via 10% SDS-PAGE and separated proteins were transferred to PVDF membranes. After blocking in 5% non-fat dry milk in TBS-T (0.1% Tween 20) for 2 h at 37˚C, the membranes were incubated with anti-ITGB8 (4˚C; 8 h; 1:500; cat. no. ab80673) and anti-GAPDH (4˚C; 8 h; 1:5,000; cat. no. ab9485) primary antibodies (both from Abcam), and then with the secondary goat anti-rabbit IgG antibodies (37˚C; 2 h; 1:5,000; cat. no. ab205718; Abcam), detection was performed using chemiluminescence (Beyotime Institute of Biotechnology). ImageJ software (1.52v; National Institutes of Health) was used to quantify protein band intensity.

 Statistical analysis

Correlation was assessed by using Pearson's correlation coefficient. Except for the Microarray, each experiment was repeated ≥ three times. Relative expression levels of ITGB8, miR-20a-5p and the five DELs were calculated using the 1.1.2 version of pcr package (https://www. rdocumentation.org/packages/pcr/versions/1.1.2). Statistical significance was determined using an unpaired Student's t-test between two groups. P<0.05 was considered to indicate a statistically significant difference. All statistical analyses were performed using SPSS (version 20.0; IBM, Corp.) and GraphPad Prism (version 8.0; GraphPad Software, Inc.) software.

 Results <sec> <title>Identifying DELs and DEGs via microarray

In order to investigate the expression pattern of lncRNAs in degenerative CE and HCE tissues, the present study used ceRNA microarray to identify DELs and DEGs. A box plot was used to demonstrate the distribution of the hybridization data and the degree of dispersion. After log2 normalization, no abnormal distributions of data were observed in the six healthy and CED samples (Fig. 1A). Analysis of correlations between expression and CE disease status showed that the expression level of genes was highly similar among Dg (>0.98) and less so for Hg (0.93-0.96; Fig. 1B). The variation between CED and HCE was demonstrated in a scatter plot of the lncRNA expression profile (Fig. 1C). A volcano plot identified the DELs at different P-values and fold-changes between the two groups (Fig. 1D). CED samples were compared with healthy samples, where a total of 369 DELs, including 316 upregulated and 53 downregulated lncRNAs (Fig. 1C), and a total of 246 DEGs, including 171 upregulated and 75 downregulated mRNAs (Fig. 1D), were identified from the microarray results.

Hierarchical clustering identified that lncRNA and gene expression levels were distinguishable between CED and healthy samples, as presented by the associated heat map (Fig. 2). The top 20 fold-change upregulated or downregulated DELs are presented in Fig. 3A and Table I. The top 20 fold-change upregulated and downregulated DEGs are presented in Fig. 3B and Table II.

 GO and KEGG pathway enrichment analysis on basis of total DEGs

In order to identify the potential etiological factors and key genes associated with the pathogenesis of CIDD, BP, MF and CC analyses in GO and KEGG enrichment pathways were performed for DEGs. The top 20 results are presented in Fig. 4 and Tables SIII-VI. BP, MF, CC and KEGG terms that were significantly associated with DEGs included ‘integrin binding (GO: 0005178)’ in Fig. 4B, ‘ossification (GO: 0001503)’ in Fig. 4A, ‘extracellular matrix (ECM)-receptor interaction (hsa04512)’ in the Fig. 4D and ‘skeletal system development (GO: 0001501)’ in the Fig. 4A.

 Construction of a ceRNA network

To evaluate the potential functions of DELs in CIDD, the same putative miRNA response element (MREs) in DELs and DEGs were identified. Using the miRanda tool, the present study determined the regulatory associations between miRNAs and lncRNAs. Data from the miRDB, TargetScan and miRTarBase databases were integrated to determine the regulatory associations between miRNAs and mRNAs (Table III). A ceRNA regulatory network was then constructed by combining the identified miRNA-mRNA and lncRNA-miRNA associations (Fig. S1). The ceRNA network was composed of 451 nodes [168 miRNAs, 226 DELs (including 189 upregulated and 37 downregulated DELs) and 57 DEGs (including 47 upregulated and 10 downregulated DEGs)] and 1,087 lines (including 852 lncRNA-miRNA connections and 235 mRNA-miRNA connections).

 GO and KEGG pathway enrichment analysis on basis of DEGs in the ceRNA network

To determine the potential functions of identified DEGs and DELs in the ceRNA network, GO and KEGG pathway enrichment analyses were performed, and the results are presented in Fig. 5 and Tables SVII-X, which included ‘Ossification’, ‘Blood vessel development’, ‘Integrin binding’ and ‘SMAD binding’. The overlapping items in both of the functional analysis results of the ceRNA network and functional analysis results of all DEGs were considered to represent significant enrichments from BP and MF analyses in GO and KEGG. The results are presented in Fig. 6 and Table IV, which included ‘Ossification’, ‘Integrin binding’, ‘Skeletal system development’ and ECM.

 Validation of lncRNAs, miRNA and gene expression in tissue samples using RT-qPCR and western blot analyses

To assess the results from the ceRNA microarray, the present study selected one DEG, ITGB8(24), miR-20a-5p (25) and its related five DELs (ENST00000548900, lnc-MYBPC1-1:1, lnc-ARL13A-1:1, lnc-C2orf67-1:1 and lnc-DNAJB6-3:1), based on a literature review of the miRNAs and genes identified in the functional analysis in the ceRNA network by searching ‘geneID’ and ‘cartilage’ on PubMed (https://pubmed.ncbi.nlm.nih.gov/) and the selected five DELs were associated with miR-20a-5p in the ceRNA network. The expression patterns of seven RNAs from the RT-qPCR results demonstrated that ITGB8 (0.45-fold change) and lnc-DNAJB6-3:1 (0.34-fold change) expression levels were significantly downregulated in CED samples compared with HCE samples. The expression levels of miR-20a-5p (2.10-fold change) was upregulated in CED samples compared with those in HCE samples (Fig. 7A). The RT-qPCR results of ITGB8 and lnc-DNAJB6-3:1 expression levels are consistent with the microarray data (Fig. 7B). Moreover, the western blotting results shown that ITGB8 is downregulated in CED samples compared with that in HCE samples (Fig. 7C and D).

 Discussion

CED is the main cause of IDD (26). Although the etiology of CED remains unknown, it is a complex, multifactorial disorder that is influenced by cartilage degeneration due to angiogenesis, inflammation, ECM degeneration, oxidative stress, cell hyperproliferation, chondrocyte apoptosis and autophagy (27,28). The study into the pathogenesis of CED is ongoing.

lncRNAs have been reported to serve essential roles in angiogenesis, inflammation and ECM degeneration (29). However, to the best of our knowledge, the expression profile of lncRNAs and their potential biological functions in CED have not yet been reported. The present study utilized a lncRNA and mRNA microarray to provide comprehensive lncRNA and mRNA expression profiles in CED and HCE, with the aim of investigating the potential involvement of dysregulated lncRNAs in the development of CED. The present study identified a large number of DELs and DEGs from diverse genomic locations, with 246 DEGs (171 upregulated and 75 downregulated) and 369 lncRNAs (316 upregulated and 53 downregulated) found to be differentially expressed between cervical vertebra CED samples and healthy controls. A ceRNA network was established, which included 168 miRNAs, 226 DELs (including 189 upregulated and 37 downregulated DELs) and 57 DEGs (including 47 upregulated and 10 downregulated DEGs). To the best of our knowledge, the present study was the first to report the comprehensive lncRNA expression profile in CED. Therefore, the results of the present study provide a novel theoretical basis for further investigations into the function of lncRNA in CED.

lncRNAs are rich in miRNA binding sites and can act as ceRNAs in cells, which leads to relief of the inhibitory effect of miRNAs on their target genes and, ultimately, increased expression of the target genes (30). To determine the potential function of lncRNAs as miRNAs sponges in cervical CED, GO BP, GO MF, GO CC and KEGG pathway enrichment analyses were performed on total DEGs and ceRNA network DEGs in the present study. The common significant enrichments between DEGs in the ceRNA network and all DEGs from GO BP, GO MF, GO CC and KEGG pathway analysis were associated with cartilage development, degeneration and regeneration.

Integrins are well-known as cell adhesion molecules and act via ECM-receptor interaction (31). Moreover, integrins are suggested to control intracellular signaling pathways both physically and chemically as mechanoreceptors (32). Previous studies have reported that if the fibrillar network of the ECM is dysfunctional or assembled from fibers that cannot bind and activate α5β1 integrins, chondrocytes increase the MMP expression, which inhibits cartilage regeneration (33). The increased expression of α2β1 and activation of integrin signaling may also induce a pathological increase in MMP expression following exposure to static compression (33,34).

The significant items in the results of enrichment analyses in the present study included ‘integrin binding’, ‘ossification’, ‘ECM’ and ‘skeletal system development’ Li et al (35) previously reported that the loss of OPG leads to the occurrence of IDD by promoting the ‘ossification’ biological process of the cartilage endplate in the osteoprotegerin (OPG)-knockout mice. Qu et al (36) reported that DEGs were primarily involved in ‘skeletal system development’ in a co-expression network in IDD and serve an important role the pathobiology of the disease.

Signaling pathways and the term ‘signal receptor binding’ are similar. Xue et al (37) previously reported that the decreased expression of NF-κB interacting lncRNA, a miR-145 sponge, inhibited SP1 expression and regulated the NF-κB signaling pathway to inhibit proliferation and promote apoptosis of chondrocytes. Furthermore, Gu et al (38) revealed that the majority of hub genes associated with osteoarthritic cartilage screened from GSE51588 were primarily enriched in ‘protein digestion and absorption’.

The present study identified a group of DEGs that serve important roles in cartilage development and degeneration from a literature review (Table V). Genes whose microarray results were consistent with those reported in the literature include ITGB8, lactotransferrin (LTF), msh homeobox 2 (MSX2), latent transforming growth factor β binding protein 2 (LTBP2), neuroblastoma suppressor of tumorigenicity 1 (NBL1) and CITED2. Genes whose microarray results are inconsistent with those reported in the literature include insulin like growth factor binding protein 5 (IGFBP5), collagen type III α 1 chain (COL3A1), COL5A1, COL5A2, fibroblast growth factor 9 (FGF9), ITGBL1, versican (VCAN), RAB23, frizzled class receptor 7 (FZD7) and ITGA11. LaPointe et al (24) reported that the downregulation of ITGB8 increased the expression of ITGA11 and caused ablation of COL2A1 expression, as detected in the chondrogenic differentiation of mesenchymal stem cells. Xue et al (39) observed that LTF can inhibit IL-1β-induced apoptosis via AKT1-induced cAMP responsive element binding protein 1 phosphorylation in human articular chondrocytes. Furthermore, Nishimura et al (40) suggested that MSX2, a homeobox Runx2 family member, may promote endochondral ossification via upregulation of Osterix (Sp7) expression. He et al (9) also reported that the transcriptional regulator CITED2 promoted chondroprotection by inhibiting the expression of MMP13, and that the decreased expression of CITED2 induced cartilage degradation via increased expression of MMP13. Moreover, Sideek et al (41) reported that exogenous LTBP2 inhibits elastinogenesis in ear cartilage chondrocytes in culture, while Goessler et al (42) reported that LTBP2 was activated during chondrogenic dedifferentiation to decrease the repairability of cartilage defects in vitro. In addition, Wei et al (43) used microarray technology to demonstrate that the attenuation of cartilage anabolism by increased expression of NBL1, a secreted BMP antagonist, is a molecular mechanism underlying cartilage degeneration.

Genes from the present study that were not consistent with those in the literature may decrease CED. For instance, Weimer et al (44) revealed that IGFBP5, a potentiator of IGF-I, could increase the expression of IGF-1 to maintain the proliferation and survival of chondrocytes. Type III collagen, encoded by COL3A1, together with type I collagen, provides a structural framework for the synovium surrounding the synovial joint to protect articular cartilage (45). Type V collagen, encoded by COL5A1 and COL5A2, is present in tissues containing type I collagen and appears to regulate the assembly of shaped fibers composed of type I and type V collagen, with the type V collagen being associated with articular cartilage maturity (46,47). Exogenous FGF9 attenuates cartilage degradation in mice (48), and its increased expression may represent a protective mechanism that remains intact during cartilage degeneration. Song et al (49) reported that overexpression of ITGBL1 inhibits integrin signaling in developing chondrocytes to promote cartilage formation. Moreover, Choocheep et al (50) observed that VCAN is required for chondrocyte differentiation during cartilage formation. Yang et al (51) also demonstrated that too much or too little RAB23 protein induced the expression of SOX9 and caused the failure of chondrogenic differentiation. In addition, this differential expression profile of RAB23 may not have been specific to chondrocyte differentiation in their study by using gain-of-function and loss-of-function experiments (51). Randall et al (52) reported that activation of FZD7, a component of the Wnt planar cell polarity pathway, was able to promote growth plate chondrocyte column formation.

Of the genes identified in the present study that were consistent with those in the literature, only four miRNAs with MREs within ITGB8 were associated with cartilage degeneration and calcification, including miR-20a-5p, miR-145-5p, miR-93-5p and miR-17-5p. Guo et al (25) reported that miR-20a-5p in exosomes from breast cancer cells promoted the proliferation and differentiation of osteoclasts. Furthermore, Chen et al (53) revealed that decreased miR-145-5p expression promoted the proliferation of osteoclasts to aggravate cartilage erosion by targeting OPG. Xiao et al (54) also found that increased osteoclasts in the CE of ovariectomized mice resulted in cartilage remodeling to stimulate IDD. It has been demonstrated that increased HOX transcript antisense RNA expression, decreased the expression of miR-17-5p to promote cartilage injury by mediating the Wnt/β-catenin pathway via fucosyltransferase 2(55). Moreover, Xue et al (56) suggested that increased expression of miR-93-5p promoted chondrocyte viability and suppressed chondrocyte apoptosis in vitro.

Previous studies have reported that the ceRNA mechanism can be used for the analysis of multiple diseases (36,57). Qu et al (36) constructed a miRNA-lncRNA-mRNA ceRNA regulatory network on the basis of DEGs and lncRNAs in patients with IDD. In addition, Wang et al (57) identified a novel ceRNA network that provides a novel insight for investigating the underlying molecular mechanism of diabetic peripheral neuropathy. In the present study, ITGB8 was inhibited by increased expression of MREs resulting from decreased expression of lncRNAs in the ceRNA network. After screening, the expression levels of IGTB8, miR-20a-5p and its five regulated lncRNAs within ceRNA network were validated in 16 CIDD and 16 healthy subjects using RT-qPCR. It was suggested that downregulated lnc-DNAJB6-3:1 may have decreased the expression of ITGB8 to promote cartilage degeneration via upregulation of miR-20a-5p expression. To the best of our knowledge, the present study was the first to report that lnc-DNAJB6-3:1 is downregulated in CED, as well as demonstrating that lnc-DNAJB6-3:1 regulated miR-20a-5p/ITGB8 axis in the CED.

In the present research, the box plots (Fig. 1A) were used to visualize the distributions of lncRNAs for two groups. Although the normalized data were skewed, the distributions of the log2 ratios among the six samples presented an excellent similarity after normalization, which indicated that the data could be further analyzed. The phenomenon that the normalized data were skewed (long tail above mean/median) was also demonstrated in previously published studies (58-60). It was suggested that the sampling error and selective bias are responsible for this phenomenon, and enlarging the sample sizes is an effective strategy to decrease the sampling error and selective bias. In the present study, gene expression was highly similar among Dg (>0.98) and less so for Hg (0.93-0.96), which can also be due to sampling error and selective bias.

Technical replicates were performed in order to further determine the reliability of detection platform for identifying DEGs using microarray. The SBC Human ceRNA array (Agilent-078298 human ceRNA array V1.0 4X180K) has demonstrated excellent repeatability in previously published papers (61-63). Moreover, the Shanghai Shibei Biotechnology Co., Ltd. has obtained official certification of Microarray repeatability and accuracy from Agilent Technologies, Inc. (R2>0.95, MAQC Plan) (61-63). The CV of the all samples in present study were <15% [Dg1, 4.23209%; Dg2, 4.88349%; Dg3, 5.57233%; Hg1, 4.19178%; Hg2, 3.85391%; Hg3, 5.49377%). Therefore, technical replicates were not conducted in the current study.

Previous studies (61-63) have improved the reliability of microarray results using biological repetition. Biological repetition for six different individuals (three samples from healthy control and three samples from patients with CIDD) was performed in the present study. Subsequently, the authors further verified the outcomes of ceRNA array using RT-qPCR in 32 samples (16 patients and 16 control individuals), and demonstrated that the ceRNA array and RT-qPCR results had an excellent reliability.

Therefore, a total of 32 samples (16 patients and 16 normal control) were used in the current study to further assess the expression levels of lnc-DNAJB6-3:1, ITGB8 and miR-20a-5p using RT-qPCR or western blot analysis. Previous reports have aimed at identifying the association between disease and potential lncRNA based on the semi-supervised learning framework and database (64-66). In future research, with the accumulation of cases and the collection of follow-up information, it will be possible to construct a clinical prognosis model of lnc-DNAJB6-3:1 in CED.

In conclusion, to the best of our knowledge, the present study was the first to use a microarray to investigate the lncRNA expression profile in CED. lnc-DNAJB6-3:1 was significantly downregulated in CED, indicating that lnc-DNAJB6-3:1 may be a candidate biomarker for CED. Collectively, the present results provide a new perspective toward an improved understanding of ceRNA-mediated gene regulation in CED and suggest a novel theoretical basis for further studies on the function of lncRNA in CED-associated IDD.

 Supplementary Material Figure S1. Construction of ceRNA network. Triangular green nodes denote miRNAs. Diamond-shaped red nodes indicate lncRNAs and oval blue nodes indicate mRNAs. ceRNA, competing endogenous RNA; miRNA/miR, microRNA; lnc, long non.coding RNA. Table SI. Details of the Dg and Hg. Table SII. Specific primer sequences for reverse transcription-quantiative PCR. Table SIII. Top 20 of GO biological process enrichment analysis for differentially expressed genes (P<0.05). Table SIV. Top 20 of GO molecular function enrichment analysis for differentially expressed genes (P<0.05). Table SV. Top 20 of GO cellular component enrichment analysis for differentially expressed genes (P<0.05). Table SVI. Top 20 of Kyoto Encyclopedia of Genes and Genomes pathway enrichment analysis for differentially expressed genes (P<0.05). Table SVII. Top 20 of GO biological process enrichment analysis for differentially expressed genes in competing endogenous RNA network (P<0.05). Table SVIII. Top 20 of GO molecular function enrichment analysis for differentially expressed genes in competing endogenous RNA network (P<0.05). Table SIX. Top 20 of GO cellular component enrichment analysis for differentially expressed genes in competing endogenous RNA network (P<0.05). Acknowledgements

The authors would like to thank Dr Wenzhou Huang and Dr Xinxin Miao (Department of Orthopedics, Second Affiliated Hospital of Nanchang University, Nanchang, China) for providing technical help, including assistance with the follow-up and advice regarding manuscript preparation.

 Funding

The present study was supported by the National Natural Science Foundation of China (grant nos. 8166090137 and 8186090165).

 Availability of data and materials

The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.

 Authors' contributions

JY conceived the study and XC designed the study. JY, RD, JZ and CP performed the experiments. JJ, TW, SH and XL collected the clinical samples and analyzed the data. JY wrote the study. All authors read and approved the final manuscript.

 Ethics approval and consent to participate

The present study was approved by the Ethics Committee of The Second Affiliated Hospital of Nanchang University. Written informed consent was obtained from all participants.

 Patient consent for publication

Not applicable.

 Competing interests

The authors declare that they have no competing interests.

 References LivshitsGPophamMMalkinISambrookPNMacgregorAJSpectorTWilliamsFMLumbar disc degeneration and genetic factors are the main risk factors for low back pain in women: The UK Twin Spine StudyAnn Rheum Dis701740174520112164641610.1136/ard.2010.137836 HurwitzELRandhawaKYuHCotePHaldemanSThe Global Spine care initiative: A summary of the global burden of low back and neck pain studiesEur Spine J2779680120182948040910.1007/s00586-017-5432-9 GantenbeinBIllien-JüngerSChanSCWalserJHaglundLFergusonSJIatridisJCGradSOrgan culture bioreactors-platforms to study human intervertebral disc degeneration and regenerative therapyCurr Stem Cell Res Ther1033935220152576419610.2174/1574888x10666150312102948 XiaKZhuJHuaJGongZYuCZhouXWangJHuangXYuWLiLIntradiscal injection of induced pluripotent stem cell-derived nucleus pulposus-like cell-seeded polymeric microspheres promotes rat disc regenerationStem Cells Int2019680654020193119167910.1155/2019/6806540 BendtsenMBüngerCEZouXFoldagerCJørgensenHSAutologous stem cell therapy maintains vertebral blood flow and contrast diffusion through the endplate in experimental intervertebral disc degenerationSpine (Phila Pa 1976)36E373E37920112137264910.1097/BRS.0b013e3181dce34c RodriguezAGSlichterCKAcostaFLRodriguez-SotoAEBurghardtAJMajumdarSLotzJCHuman disc nucleus properties and vertebral endplate permeabilitySpine (Phila Pa 1976)3651252020112124004410.1097/BRS.0b013e3181f72b94 YuanWCheWJiangYQYuanFLWangHRZhengGLLiXLDongJEstablishment of intervertebral disc degeneration model induced by ischemic sub-endplate in rat tailSpine J151050105920152563746610.1016/j.spinee.2015.01.026 ColombierPClouetJHamelOLescaudronLGuicheuxJThe lumbar intervertebral disc: From embryonic development to degenerationJoint Bone Spine8112512920142393272410.1016/j.jbspin.2013.07.012 HeZLeongDJXuLHardinJAMajeskaRJSchafflerMBThiMMYangLGoldringMBCobelliNJSunHBCITED2 mediates the cross-talk between mechanical loading and IL-4 to promote chondroprotectionAnn N Y Acad Sci144212813720193089176610.1111/nyas.14021 LuWCaoFWangSShengXMaJLncRNAs: The regulator of glucose and lipid metabolism in tumor cellsFront Oncol9109920193185018910.3389/fonc.2019.01099 WangHHuoXYangXRHeJChengLWangNDengXJinHWangNWangCSTAT3-mediated upregulation of lncRNA HOXD-AS1 as a ceRNA facilitates liver cancer metastasis by regulating SOX4Mol Cancer1613620172881092710.1186/s12943-017-0680-1 ChenJLiuAWangZWangBChaiXLuWCaoTLiRWuMLuZLINC00173.v1 promotes angiogenesis and progression of lung squamous cell carcinoma by sponging miR-511-5p to regulate VEGFA expressionMol Cancer199820203247364510.1186/s12943-020-01217-2 WangKLongBZhouLYLiuFZhouQYLiuCYFanYYLiPFCARL lncRNA inhibits anoxia-induced mitochondrial fission and apoptosis in cardiomyocytes by impairing miR-539-dependent PHB2 downregulationNat Commun5359620142471010510.1038/ncomms4596 ChenXYanCCZhangXYouZHLong non-coding RNAs and complex diseases: From experimental results to computational modelsBrief Bionform1855857620172734552410.1093/bib/bbw060 ChenXSunYZGuanNNQuJHuangZAZhuZXLiJQComputational models for lncRNA function prediction and functional similarity calculationBrief Funct Genomics18588220193024750110.1093/bfgp/ely031 KanehisaMFurumichiMTanabeMSatoYMorishimaKKEGG: New perspectives on genomes, pathways, diseases and drugsNucleic Acids Res45D353D36120172789966210.1093/nar/gkw1092 MiHMuruganujanAHuangXEbertDMillsCGuoXThomasPDProtocol update for large-scale genome and gene function analysis with the PANTHER classification system (v.14.0)Nat Protoc1470372120193080456910.1038/s41596-019-0128-8 YuGWangLGHanYHeQYclusterProfiler: An R package for comparing biological themes among gene clustersOMICS1628428720122245546310.1089/omi.2011.0118 ChenYWangXmiRDB: An online database for prediction of functional microRNA targetsNucleic Acids Res48D127D13120203150478010.1093/nar/gkz757 ChouCHShresthaSYangCDChangNWLinYLLiaoKWHuangWCSunTHTuSJLeeWHmiRTarBase update 2018: A resource for experimentally validated microRNA-target interactionsNucleic Acids Res46D296D30220182912617410.1093/nar/gkx1067 AgarwalVBellGWNamJWBartelDPPredicting effective microRNA target sites in mammalian mRNAsElife4e0500520152626721610.7554/eLife.05005 ShannonPMarkielAOzierOBaligaNSWangJTRamageDAminNSchwikowskiBIdekerTCytoscape: A software environment for integrated models of biomolecular interaction networksGenome Res132498250420031459765810.1101/gr.1239303 LivakKJSchmittgenTDAnalysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) methodMethods2540240820011184660910.1006/meth.2001.1262 LaPointeVLVerpoorteAStevensMMThe changing integrin expression and a role for integrin β8 in the chondrogenic differentiation of mesenchymal stem cellsPLoS One8e8203520132431240010.1371/journal.pone.0082035 GuoLZhuYLiLZhouSYinGYuGCuiHBreast cancer cell-derived exosomal miR-20a-5p promotes the proliferation and differentiation of osteoclasts by targeting SRCIN1Cancer Med85687570120193138546410.1002/cam4.2454 Berg-JohansenBHanMFieldsAJLiebenbergECLimBJLarsonPEGunduz-DemirCKazakiaGJKrugRLotzJCCartilage endplate thickness variation measured by ultrashort Echo-time MRI is associated with adjacent disc degenerationSpine (Phila Pa 1976)43E592E60020182898473310.1097/BRS.0000000000002432 ZhangJFWangGLZhouZJFangXQChenSFanSWExpression of matrix metalloproteinases, tissue inhibitors of metalloproteinases, and interleukins in vertebral cartilage endplateOrthop Surg1030631120183047432410.1111/os.12409 ZhangYHeFChenZSuQYanMZhangQTanJQianLHanYMelatonin modulates IL-1β-induced extracellular matrix remodeling in human nucleus pulposus cells and attenuates rat intervertebral disc degeneration and inflammationAging (Albany NY)11104991051220193177214510.18632/aging.102472 XuTWuKZhangLZhengSWangXZuoHWuXTaoGJiangBZhangLLong non-coding RNA LINC00858 exerts a tumor-promoting role in colon cancer via HNF4alpha and WNK2 regulationCell Oncol (Dordr)4329731020193188457710.1007/s13402-019-00490-8 LiaoCLongZZhangXChengJQiFWuSHuangTLncARSR sponges miR-129-5p to promote proliferation and metastasis of bladder cancer cells through increasing SOX4 expressionInt J Biol Sci1611120203189284110.7150/ijbs.39461 AlankoJIvaskaJEndosomes: Emerging platforms for integrin-mediated FAK signallingTrends Cell Biol2639139820162694477310.1016/j.tcb.2016.02.001 HiroseNOkamotoYYanoshitaMAsakawaYSumiCTakanoMNishiyamaSSuSCMitsuyoshiTKunimatsuRProtective effects of cilengitide on inflammation in chondrocytes under excessive mechanical stressCell Biol Int4496697420203187632310.1002/cbin.11293 YanZPanYWangSChengMKongHSunCHuKChenTDongQChenJStatic Compression induces ECM remodeling and integrin α2β1 expression and signaling in a rat tail caudal intervertebral disc degeneration modelSpine (Phila Pa 1976)42E448E45820172754857910.1097/BRS.0000000000001856 Almonte-BecerrilMGimenoLIVillarroyaOBenito-JardonMKouriJBCostellMGenetic abrogation of the fibronectin-α5β1 integrin interaction in articular cartilage aggravates osteoarthritis in micePLoS One13e019855920182987055210.1371/journal.pone.0198559 LiXFXueCCZhaoYJChengSDZhaoDFLiangQQChenLWangQLuSShiQDeletion of Opg leads to increased neovascularization and expression of inflammatory cytokines in the lumbar intervertebral disc of miceSpine (Phila Pa 1976)42E8E1420172719601610.1097/BRS.0000000000001701 QuZQuanZZhangQWangZSongQZhuangXFuCXuFLiuYWangYComprehensive evaluation of differential lncRNA and gene expression in patients with intervertebral disc degenerationMol Med Rep181504151220182990113610.3892/mmr.2018.9128 XueHYuPWangWZNiuYYLiXThe reduced lncRNA NKILA inhibited proliferation and promoted apoptosis of chondrocytes via miR-145/SP1/NF-κB signaling in human osteoarthritisEur Rev Med Pharmacol Sci2453554820203201695510.26355/eurrev_202001_20030 GuHYYangMGuoJZhangCLinLLLiuYWeiRXIdentification of the biomarkers and pathological process of osteoarthritis: Weighted Gene Co-expression network analysisFront Physiol1027520193094105910.3389/fphys.2019.00275 XueHTuYMaTLiuXWenTCaiMXiaZMeiJLactoferrin inhibits IL-1β-induced chondrocyte apoptosis through AKT1-induced CREB1 activationCell Physiol Biochem362456246520152627944710.1159/000430206 NishimuraRHataKMatsubaraTWakabayashiMYonedaTRegulation of bone and cartilage development by network between BMP signalling and transcription factorsJ Biochem15124725420122225344910.1093/jb/mvs004 SideekMAMenzCParsiMKGibsonMALTBP-2 competes with tropoelastin for binding to fibulin-5 and heparin, and is a negative modulator of elastinogenesisMatrix3411412320142414880310.1016/j.matbio.2013.10.007 GoesslerURBugertPBiebackKDemlMSadickHHormannKRiedelFIn-vitro analysis of the expression of TGFbeta-superfamily-members during chondrogenic differentiation of mesenchymal stem cells and chondrocytes during dedifferentiation in cell cultureCell Mol Biol Lett10345362200516010298 WeiTKulkarniNHZengQQHelveringLMLinXLawrenceFHaleLChambersMGLinCHarveyAAnalysis of early changes in the articular cartilage transcriptisome in the rat meniscal tear model of osteoarthritis: Pathway comparisons with the rat anterior cruciate transection model and with human osteoarthritic cartilageOsteoarthr Cartil18992100020102043457410.1016/j.joca.2010.04.012 WeimerAMadryHVenkatesanJKSchmittGFrischJWezelAJungJKohnDTerwilligerEFTrippelSBCucchiariniMBenefits of recombinant adeno-associated virus (rAAV)-mediated insulinlike growth factor I (IGF-I) overexpression for the long-term reconstruction of human osteoarthritic cartilage by modulation of the IGF-I axisMol Med1834635820122216039210.2119/molmed.2011.00371 HuYWuRLiHGuYWeiWExpression and significance of metalloproteinase and collagen in vaginal wall tissues of patients with pelvic organ prolapseAnn Clin Lab Sci47698705201729263043 WillardKMannionSSaundersCJCollinsMSeptemberAVThe interaction of polymorphisms in extracellular matrix genes and underlying miRNA motifs that modulate susceptibility to anterior cruciate ligament ruptureJ Sci Med Sport21222820182892797110.1016/j.jsams.2017.08.017 WuJJWeisMAKimLSCarterBGEyreDRDifferences in chain usage and cross-linking specificities of cartilage type V/XI collagen isoforms with age and tissueJ Biol Chem2845539554520091910359010.1074/jbc.M806369200 ZhouSWangZTangJLiWHuangJXuWLuoFXuMWangJWenXExogenous fibroblast growth factor 9 attenuates cartilage degradation and aggravates osteophyte formation in post-traumatic osteoarthritisOsteoarthritis Cartilage242181219220162747355810.1016/j.joca.2016.07.005 SongEKJeonJJangDGKimHESimHJKwonKYMedina-RuizSJangHJLeeARRhoJGITGBL1 modulates integrin activity to promote cartilage formation and protect against arthritisSci Transl Med10eaam748620183030545410.1126/scitranslmed.aam7486 ChoocheepKHatanoSTakagiHWatanabeHKimataKKongtawelertPWatanabeHVersican facilitates chondrocyte differentiation and regulates joint morphogenesisJ Biol Chem285211142112520102040434310.1074/jbc.M109.096479 YangLClintonJMBlackburnMLZhangQZouJZielinska-KwiatkowskaATangBLChanskyHARab23 regulates differentiation of ATDC5 chondroprogenitor cellsJ Biol Chem283106491065720081821862010.1074/jbc.M706795200 RandallRMShaoYYWangLBallockRTActivation of Wnt planar cell polarity (PCP) signaling promotes growth plate column formation in vitroJ Orthop Res301906191420122267435110.1002/jor.22152 ChenYWangXYangMRuanWWeiWGuDWangJGuoXGuoLYuanYmiR-145-5p Increases osteoclast numbers in vitro and aggravates bone erosion in collagen-induced arthritis by targeting osteoprotegerinMed Sci Monit245292530020183005949110.12659/MSM.908219 XiaoZFHeJBSuGYChenMHHouYChenSDLinDKOsteoporosis of the vertebra and osteochondral remodeling of the endplate causes intervertebral disc degeneration in ovariectomized miceArthritis Res2020720183020105210.1186/s13075-018-1701-1 HuJWangZShanYPanYMaJJiaLLong non-coding RNA HOTAIR promotes osteoarthritis progression via miR-17-5p/FUT2/β-catenin axisCell Death Dis971120182990776410.1038/s41419-018-0746-z XueHTuYMaTWenTYangTXueLCaiMWangFGuanMmiR-93-5p attenuates IL-1β-induced chondrocyte apoptosis and cartilage degradation in osteoarthritis partially by targeting TCF4Bone12312913620193093029410.1016/j.bone.2019.03.035 WangCXuXChenJKangYGuoJDuscherDYangXGuoGRenSXiongHThe construction and analysis of lncRNA-miRNA-mRNA competing endogenous RNA network of Schwann cells in diabetic peripheral neuropathyFront Bioeng Biotechnol849020203252394310.3389/fbioe.2020.00490 ZhangSChangYYGongYWGaoYJGuoQWangYHZhaoYLWangZPComprehensive analysis of microRNA-messenger RNA regulatory network in gemcitabine-resistant bladder cancer cellsJ Cell Biochem1206347636020193030454910.1002/jcb.27922 CaoMZhangLWangJHZengHPengYZouJShiJZhangLLiYYoshidaSIdentifying circRNA-associated-ceRNA networks in retinal neovascularization in miceInt J Med Sci161356136520193169291710.7150/ijms.35149 Wang S, Zhan J, Lin X, Wang Y, Wang Y and Liu Y: CircRNA-0077930 from hyperglycaemia-stimulated vascular endothelial cell exosomes regulates senescence in vascular smooth muscle cells. Cell Biochem Funct: Apr 19 doi: 10.1002/cbf.3543.2020 (Epub ahead of print). WangRZhangSChenXLiNLiJJiaRPanYLiangHEIF4A3-induced circular RNA MMP9 (circMMP9) acts as a sponge of miR-124 and promotes glioblastoma multiforme cell tumorigenesisMol Cancer1716620183047026210.1186/s12943-018-0911-0 GaoAGongYZhuCYangWLiQZhaoMMaSLiJHaoSChengHChengTBone marrow endothelial cell-derived interleukin-4 contributes to thrombocytopenia in acute myeloid leukemiaHaematologica1041950196120193079220010.3324/haematol.2018.214593 ZhaoXLiDYangFLianHWangJWangXFangESongHHuAGuoYLong noncoding RNA NHEG1 drives β-catenin transactivation and neuroblastoma progression through interacting with DDX5Mol Ther2894696220203198203710.1016/j.ymthe.2019.12.013 ChenXYanGNovel human lncRNA-disease association inference based on lncRNA expression profilesBioinformatics292617262420132400210910.1093/bioinformatics/btt426 HerszageJCensorNModulation of learning and memory: A shared framework for interference and generalizationNeuroscience39227028020183012568610.1016/j.neuroscience.2018.08.006 GershmanSJDawNDReinforcement learning and episodic memory in humans and animals: An integrative frameworkAnnu Rev Psychol6810112820172761894410.1146/annurev-psych-122414-033625

Expression profile of lncRNAs detected using a microarray in degeneration CE compared with healthy CE. (A) Box plots were used to visualize the distributions of lncRNAs for the two groups. After normalization, the distributions of the log2 ratios among six samples were nearly the same. (B) A block diagram for the correlation coefficient r. The blue or red of the color represents the degree of correlation between samples, and the deeper the red, the higher the correlation between the two samples. (C) Volcano plot visualizing differential lncRNA expression between the two groups. The vertical lines correspond to a 2.0-FC (log2 scaled), showing upregulation and downregulation. The horizontal line represents a P<0.05 (log10 scaled). The red and blue points in the plot present the statistically significantly upregulated and downregulated lncRNAs, respectively. Red points mark upregulated expression of lncRNAs in CED vs. healthy CE (P<0.05, FC ≥2). Blue points mark downregulated expression of lncRNAs in CED vs. healthy CE (P<0.05, FC ≤0.5). (D) Scatter plot demonstrating the variation of lncRNAs expression in Dg (y-axis) vs. Hg (x-axis). The values of the x- and y-axes are the averaged normalized signal values of each group (log2 scaled). The middle line represents no difference between the two groups. In the figure, representing the probe point in the two groups, the signal difference was FC ≥2; red represents upregulation and blue represents downregulation. FC, fold-change; Dg, degeneration group; Hg, healthy group; lnc, long non-coding RNA.

Heat map of DEGs and DELs. Heat map presenting a distinguishable expression profile of (A) long non-coding RNAs and (B) mRNAs between the Dg and Hg. The values correspond to the different colors representing the fold change (log2 transformed) of each sample. Black stands for 0, indicating no change in gene expression; red represents upregulation; and green represents downregulation. The brightness of the color represents the degree of increased or decreased gene expression. Dg, degeneration group; Hg, healthy group.

The top 20 upregulated and downregulated DEGs and DELs. Heat map of the top 20 upregulated and downregulated (A) lncRNAs and (B) genes. Fold-change ≥2 and P<0.05. Dg, degeneration group; Hg, healthy group; lnc, long non-coding RNA.

Top 20 of GO (BP, MF and CC) and KEGG pathway enrichment analysis based on total DEGs. (A) GO BP. (B) GO MF. (C) GO CC. (D) KEGG pathway. GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; BP, Biological Processes; MF, Molecular Function; CC, Cellular Component.

Functional enrichment analysis based on total DEGs. (A) GO BP. (B) GO MF. (C) GO CC. (D) KEGG pathway. The overlapping terms between Figs. 4 and 5 were marked in red. GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; BP, Biological Processes; MF, Molecular Function; CC, Cellular Component.

Chord plot displaying the significant enrichments of GO and KEGG pathway enrichments associated with cartilage. GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; ECM, extracellular matrix; FC, fold change.

Expression levels of ITGB8, miR-20a-5p, ENST00000548900, lnc-MYBPC1-1:1, lnc-ARL13A-1:1, lnc-KANSL1L-1:1 and lnc-DNAJB6-3:1 in CED samples compared with healthy CE. (A) Results of RT-qPCR were evaluated using the 2-ΔΔCq method. (B) RT-qPCR results of miR-20a-5p was evaluated using the 2-ΔΔCq method. (C) Similarity of RT-qPCR and microarray analysis results of ITGB8 and the five DEGs. (D) Western blotting results of ITGB8 expression, which were (E) semi-quantified using ImageJ software. Data are presented as the mean ± standard deviation. *P<0.05 and **P<0.01 vs. HG.

Top 20 upregulated and downregulated lncRNAs in degenerative CE vs. healthy CE were sorted according to FC.

lncRNA-ID Log2FC P-value Regulation
ENST00000460164 -4.0654447321111 0.010925989 Down
ENST00000444114 3.15320239767769 0.005654941 Up
lnc-ARHGAP26-1:3 -2.97917502630228 0.035859204 Down
ENST00000518605 -2.80565656425262 0.02692457 Down
lnc-AC068473.1-1:1 2.78060506466701 0.026317694 Up
NR_033909 -2.77339818775925 0.01209491 Down
NONHSAT013575 2.72406891252417 0.006270757 Up
ENST00000606068 2.71093666603456 0.007273469 Up
lnc-OTOP1-1:1 2.6826499083849 0.004281053 Up
lnc-AC073416.2-3:1 -2.55407843103034 0.022005983 Down
lnc-AC233263.1-12:1 -2.54967640076035 0.03761393 Down
NR_024469 2.52303431159003 0.000477517 Up
lnc-AC073416.2-5:2 -2.51306945648041 0.017447387 Down
lnc-FZD6-1:2 2.50613392195141 0.014306711 Up
lnc-BNIP1-3:3 -2.47300950984058 0.041505299 Down
lnc-BNIP1-4:3 -2.4537777630073 0.017295824 Down
lnc-LONRF2-1:3 2.30742206138951 1.99172E-06 Up
ENST00000517816 -2.30653135357461 0.048332914 Down
NR_120454 2.30548698577468 0.030404793 Up
NR_030732 2.29451834582489 0.019540145 Up

In the current analysis, FC referred to the absolute ratio (not log scale) of normalized intensities between two groups. FC ≥2 and P<0.05. lncRNA, long non-coding RNA; CE, cartilage endplate; FC, fold-change.

Top 20 upregulated and downregulated genes in degenerative CE vs. healthy CE were sorted according to FC.

Gene symbol log2FC P-value Regulation
SCGB2A2 -6.38069654020843 0.016239341 Down
HLA-DRB5 -5.60888034944606 0.000344626 Down
SYT7 -4.12512677137271 0.010160235 Down
DEFA4 -3.84312421870298 0.03700143 Down
KLK11 -3.73246320843984 0.003181221 Down
HBD -3.56787972602964 0.030431781 Down
PLA2G2A 3.30813047877424 0.015774342 Up
MMP8 -3.288405759072 0.004677465 Down
RNASE2 -3.22555107120361 0.00318142 Down
HLA-DQB1 -3.18445779878179 0.019791381 Down
CEACAM8 -3.16284892181939 0.015625014 Down
CCNA1 -3.05409089387395 0.040614958 Down
GINS3 3.03459748011711 0.034104828 Up
BIRC7 2.99886168508531 0.026568151 Up
CLC -2.8225612164875 0.030603049 Down
GPR20 2.8001431232719 0.039633617 Up
LCN2 -2.7633057628854 0.015383988 Down
CEBPE -2.76040800943709 0.008924767 Down
HP -2.73937831942498 0.029404319 Down
LTF -2.7084261374062 0.031266564 Down

In the present analysis, FC referred to the absolute ratio (not log scale) of normalized intensities between two groups. FC ≥2 and P<0.05. CE, cartilage endplate; FC, fold-change.

Relationship between mRNAs and miRNA predicted by all three databases (miRTarBase, TargetScan and miRDB) in competitive endogenous RNA network.

Gene miR
ABHD13 miR-6507-5p
ADAMTS14 miR-4667-5p, miR-423-5p, miR-6762-5p, miR-8089, miR-6845-5p, miR-4700-5p
ANKRD50 miR-93-5p, miR-20b-5p, miR-20a-5p, miR-3934-3p, miR-548t-5p, miR-153-5p, miR-6736-3p, miR-106b-5p, miR-17-5p, miR-6884-3p, miR-106a-5p
BNC2 miR-567, miR-513c-5p, miR-514b-5p, miR-8485, miR-6867-5p
CCDC80 miR-3916
CDCP1 miR-1233-5p, miR-6778-5p, miR-665, miR-6878-5p, miR-1827
CHL1 miR-182-5p
CITED2 miR-3121-3p, miR-548t-3p, miR-548t-5p, miR-548aa, miR-182-5p
COL3A1 let-7b-5p, miR-29b-3p, miR-767-5p
COL5A1 miR-6785-5p, miR-4728-5p, miR-6878-5p, miR-4516, miR-6883-5p, miR-149-3p
COL5A2 miR-29b-3p, miR-143-3p, miR-767-5p
COLEC12 miR-6770-5p, miR-7-5p
COMMD3-BMI1 miR-195-3p, miR-3120-3p
DYRK2 miR-6767-3p, miR-193a-3p, miR-484, miR-193b-3p, miR-8077, miR-24-3p, miR-3613-3p, miR-8485, miR-4447, miR-1275
EBF1 miR-4776-3p
FGF9 miR-942-5p, miR-182-5p, miR-622, miR-140-5p, miR-8485
FOXN2 miR-6124, miR-6809-3p, miR-367-3p, miR-25-3p, miR-363-3p, miR-3613-3p, miR-92b-3p, miR-5692a, miR-92a-3p, miR-7109-3p, miR-188-5p, miR-4724-5p
FZD1 miR-7109-3p
FZD7 miR-145-5p
GPR173 miR-8485, miR-6785-5p, miR-4524a-3p, miR-4516, miR-6883-5p, miR-92a-2-5p, miR-1249-5p, miR-4728-5p, miR-6825-5p, miR-6797-5p, miR-149-3p, miR-3202, miR-6867-5p
HES6 miR-4739, miR-6129, miR-6127, miR-4510, miR-6133, miR-6130
IGFBP5 miR-548aa, miR-185-5p, miR-3613-3p, miR-3925-5p, miR-143-3p, miR-548t-3p, miR-4755-5p, miR-4446-5p
ITGA11 miR-1827, miR-6867-5p, miR-650, miR-3612
ITGB8 miR-3613-3p, miR-93-5p, miR-20a-5p, miR-145-5p, miR-6507-5p, miR-17-5p
ITGBL1 miR-8485
KCNJ12 miR-1915-3p, miR-6764-5p, miR-8485
KCNQ5 miR-543, miR-6507-5p, miR-6867-5p
LTBP2 miR-6880-5p, miR-4708-3p, miR-3165
LTF miR-214-3p
MSX2 miR-5196-5p, miR-4747-5p
NBL1 miR-4270, miR-6754-5p, miR-4441
NKX2-5 miR-1538, miR-4745-3p
NOVA1 miR-338-3p
NUAK2 miR-589-5p, miR-3925-5p
OLFM4 miR-486-5p
PDZRN4 miR-7844-5p, miR-3613-3p
PHYHIP miR-6785-5p, miR-1296-3p, miR-6883-5p, miR-4728-5p, miR-6825-5p, miR-3616-3p, miR-149-3p, miR-6081, miR-6731-5p, miR-8085
PLSCR4 miR-15b-5p
PLXDC1 miR-4732-3p, miR-670-5p, miR-150-5p
PTPRF miR-24-3p, miR-298
RAB23 miR-545-3p, miR-1468-3p, miR-518a-5p, miR-424-5p, miR-527, miR-195-5p, miR-497-5p, miR-548p
S1PR1 miR-363-3p
SLC4A1 miR-504-3p, miR-4430, miR-5698, miR-3652
SNTB2 miR-20a-5p, miR-17-5p, miR-20b-5p, miR-5582-5p, miR-106b-5p, miR-93-5p
SPATS2L miR-103a-3p
SPOCK2 miR-8485
SPTLC3 miR-4534
SVOP miR-532-3p, miR-3612, miR-4667-3p, miR-4727-5p, miR-3972, miR-6852-5p, miR-1202, miR- 6749-3p, miR-7162-5p, miR-650, miR-1827, hsa-miR-4768-3p
SYT7 miR-4632-5p, miR-3202, miR-4436b-3p, miR-6735-5p, miR-6760-5p, miR-4270, miR-6879-5p, miR-1207-5p, miR-7843-5p, miR-4516, miR-4441, hsa-miR-6887-3p, hsa-miR-4763-3p
SYTL4 miR-8485, miR-4789-3p
TBXA2R miR-31-5p, miR-1275, miR-6779-5p, miR-149-3p, miR-6785-5p, miR-508-5p, miR-1273h-5p, miR-6883-5p, miR-7160-5p, miR-4478, miR-7106-5p, hsa-miR-4728-5p
THSD7A miR-3617-5p, miR-153-5p, miR-5581-5p, miR-4297, miR-6867-5p, miR-641
THY1 miR-6778-3p, miR-6825-5p
TMEM119 miR-2114-3p
TSPAN12 miR-140-5p, miR-196a-5p
VCAN miR-103a-3p, miR-107, miR-578
VWA1 miR-765, miR-6825-5p

miR/miRNA, microRNA.

Majority of significant enrichments and genes of Gene Ontology and Kyoto Encyclopedia of Genes and Genomes pathway enrichments.

Enrichment terms Gene
Extracellular matrix organization ADAMTS14, CCDC80, COL3A1, COL5A1, COL5A2, ITGA11, ITGB8, SPOCK2, VCAN, VWA1
Ossification COL5A2, FGF9, FZD1, IGFBP5, ITGA11, LTF, MSX2, S1PR1, TMEM119, VCAN
Skeletal system development FGF9, COL3A1, VWA1, VCAN, COL5A2, TMEM119, MSX2, RAB23, CITED2, ITGB8, LTF
Extracellular matrix structural constituent COL3A1, COL5A1, COL5A2, LTBP2, VCAN, VWA1
Integrin binding COL5A1, COL3A1, THY1, ITGBL1, ITGB8
Signaling receptor binding FZD7, NBL1, COL5A1, FGF9, FZD1, COL3A1, THY1, ITGBL1, PLSCR4, ITGB8, SYTL4, S1PR1
Extracellular matrix VCAN, FGF9, CCDC80, SPOCK2, LTBP2, VWA1, COLEC12, COL3A1, COL5A2, COL5A1
Extracellular region VCAN, SLC4A1, FGF9, CDCP1, CHL1, PLXDC1, OLFM4, CCDC80, SPOCK2, PTPRF, ITGB8, LTBP2, NBL1, VWA1, THSD7A, LTF, ADAMTS14, COLEC12, SYT7, COL3A1, ITGBL1, THY1, COL5A2, IGFBP5, COL5A1, CHL1
Extracellular space VCAN, SLC4A1, FGF9, CHL1, PLXDC1, OLFM4, SPOCK2, PTPRF, ITGB8, LTBP2, NBL1, VWA1, LTF, COLEC12, SYT7, COL3A1, THY1, COL5A2, IGFBP5, COL5A1, CHL1
Protein digestion and absorption COL3A1, COL5A1, COL5A2
Cell adhesion molecules ITGB8, PTPRF, VCAN
Extracellular matrix-receptor interaction ITGA11, ITGB8

Identifying crucial genes using a literature review.

Gene group Upregulated genes in cartilage degeneration Downregulated genes in cartilage degeneration
A MSX2, LTBP2, NBL1 ITGB8, LTF, CITED2
B IGFBP5, COL3A1, COL5A1, COL5A2, FGF9, ITGBL1, VCAN, RAB23, FZD7, ITGA11 -

A, Genes whose microarray results are consistent with those reported in the literature; B, Genes whose microarray results are inconsistent with those reported in the literature.