Introduction
Spinal disc herniation is known as degenerated discs protruding into the spinal canal or the foraminal canal (1). The associated symptoms, such as lower back pain and radiating leg pain (2) tend to severely affect the quality of life of patients and impose a massive economic burden on society (3,4). Herniated discs may exhibit severe ossification which may adhere to the epidural membrane, and can lead to ossification of the epidural membrane with a high risk of cerebrospinal fluid leakage and this complicates surgical treatment (5). However, the mechanisms of disc ossification remain unclear.
Herniated disc ossification is related to disc degeneration. Patients with a prolonged history of spinal disc herniation are susceptible to disc ossification (6). Degenerated disc cells have been shown to express bone morphogenetic proteins (BMPs) and their receptors (7). Although some studies have documented the BMP-induced stimulation of chondrogenesis and the regeneration of the extracellular matrix of the intervertebral disc (8,9), another study demonstrated no effect of BMPs in preventing the degeneration of the intervertebral disc, but demonsrated that it led to intervertebral disc calcification (10). The expression of runt-related transcription factor 2 (Runx2), a key transcription factor for chondrocyte hypertrophic differentiation (11), has been shown to be upregulated in the degenerative discs of dogs (12) and humans (13). Unlike dystrophic calcification, disc ossification is a positive osteogenic process associated with osteoblasts (14); however, the origin of osteogenic cells remains unclear. According to Risbud et al, these osteogenic cells differentiate from the initial cells residing in intervertebral disc tissue (15). Studies have demonstrated the osteogenic potential of annulus fibrosus cells which leads to ossification by osteogenic inductive stimulation (16,17). Owing to the linkage between ossification and degeneration, the local environmental factors appear to be instrumental in causing disc ossification.
The association of disc degeneration (or herniation) and inflammation is well documented, with several inflammatory factors now known to be implicated in the degeneration of intervertebral discs (18). However, whether these inflammatory factors are the cause of degeneration or the result of degeneration is yet to be completely understood. The presence of local inflammation and angiogenesis in herniated and ossified discs indicates that the ossification process may be related to these two processes (19). However, whether the inflammatory and angiogenic factors involved in disc herniation are akin to those in ossified discs has yet to be determined.
High-throughput sequencing is a useful method to detect the underlying molecular mechanisms of any disease and can provide valuable pointers to guide future research. There are very few studies which have applied bioinformatics analysis and identified group of genes and single nucleotide polymorphisms (SNPs) associated with intervertebral disc degeneration and disc-related disorders (20). These studies have mainly centered on the genes involved in degeneration of the intervertebral discs, and do not dwell on the osteogenic process in the degenerated discs. Screening for differentially expressed genes (DEGs) focuses on the mRNA expression in a tissue, and can allow for the detection of different genes between 2 groups, and which in turn may help to unravel the probable mechanisms of a certain pathogenesis at the gene level.
The purpose of this study was to detect the DEGs between ossified and non-ossified herniated discs. In addition, bioinformatics tools were employed to identify potential candidate genes and pathways for further research.
Materials and methods
Case selection and grouping
Patients diagnosed with lumbar intervertebral disc herniation or herniation combined with spinal canal stenosis at the Peking University Third hospital from February, 2015 to May, 2015 were included in the study. All cases were scheduled for lumbar spinal surgery for the first time. We sampled 6 discs from 6 patients with ossified intervertebral discs as the experiment group. After sampling, these cases were confirmed by computed tomography (CT) scan and verified by micro CT scan as shown in Fig. 1. Another 6 iscs from 6 patients suffering from intervertebral disc herniation without disc ossification served as the control group. Patients with gout, hypercalcemia, congenital bone metabolism disorders or bone malformation were excluded from the study. This study was approved by Medical Scientific Research Ethics Committee of Peking University Third Hospital, Beijing, China with the certification number 2015006. Written informed consent was obtained from all patients prior to obtaining the samples.
Reagents
TRIzol reagent (15596-026) for total RNA extraction was obtained from Invitrogen, Life Technologies (Carlsbad, CA, USA). The RevertAid first-strand cDNA synthesis kit (K1633) from ThermoScientific (Waltham, MA, USA) was used for reverse transcription. FastStart Universal SYBR-Green Master from Roche (Basel, Switzerland) was used for carrying out quantitative polymerase chain reaction (qPCR). Primers were synthesized by Invitrogen Biotechnology Co., Ltd. (Shanghai, China).
Sampling
After the laminas were removed by the posterior lumbar spinal procedure, the nerve root and the posterior longitudinal ligament of the herniated side of the intervertebral disc was retracted, leaving space for a circular incision to be made in the annular fibrosus. A nucleus pulposus clamp was used to extract the first piece of the herniated disc, and used as a specimen. For the purpose of uniformity, totally ossified discs which could not be sampled by this method were excluded. Following their excision, the specimens were immediately kept in frozen storage tubes and placed in liquid nitrogen. As the pre-procedure CT scan images had already been obtained, the specimens were directly transported for the micro CT scan and later stored at −80°C until further analyses.
Micro CT scan
A Siemens Inveon Micro CT scanner (Siemens Medical Solutions, Knoxville, TN, USA) was used to analyze the specimens using the following parameters: X-ray beam voltage, 80 kV; current, 500 μA; and effective resolution 13.6 μm. Referring to the semi-quantitative grading criteria of Rutges et al (14), the following criteria were adopted for this analysis: the absence of calcification was indicated as -; the presence of a single area of calcification as ±; the presence of 2 clear areas of calcification as +; and the clear presence of multiple areas of calcification as ++. We designated - or ± for a negative CT scan as the control group; ++ with positive CT scan as experiment group. At least 2 of the authors collaborated to assess the ossification from the CT radiograph and the ossification grade according to the micro CT analysis.
mRNA extraction
Following micro-CT analysis, the ossified disc group was considered as the experiment group, and the degenerated herniated disc group without ossification as the control group. For the mRNA extraction, the specimens were treated with TRIzol reagent and grinded sufficiently. The specimens were then centrifuged (8000 × g, 4°C) and reconstituted in methenyltrichloride and propyl alcohol. The total RNA was stored at −80°C for further sequencing and verification.
Sequencing and bioinformatics analysis
Sequencing was performed at the Beijing Genomics Institute (BGI). The total RNA samples were treated with DNase I to avoid DNA contamination. The enriched mRNA was mixed and fragmented into short fragments using fragmentation buffer. After the double-strand cDNA fragments were synthesized and purified, end reparation and 3′-end single nucleotide A (adenine) addition was performed. Finally, the sequencing adaptors to the fragments were ligated. Following enrichment by PCR amplification, the fragments were sequenced using a Illumina HiSeq™ 2000 sequencer (Illumina Inc., Santiago, CA, USA). Primary sequencing data generated by Illumina HiSeq™ 2000 was referred to as raw reads. The raw reads are filtered into clean reads which were aligned to the reference sequences subsequently by using the Burrows-Wheeler Alignment BWA (21)/Bowtie2 (22) tool. The NOISeq (23) method was used to screen DEGs between 2 groups. Furthermore, an in depth analysis using bioinformatics tools based on the DEGs was performed, including GO enrichment analysis, KEGG pathway enrichment analysis, and protein-protein interaction network analysis. After mapping all the DEGs to GO terms according to the database in the website, http://www.geneontology.org/, the numbers for each GO term were calculated; the significantly enriched GO terms were found by using 'GO::TermFinder' tool on the website, http://www.yeastgenome.org/help/analyze/go-term-finder. All the DEGs annotated in the GO database were used to perform GO functional classification using WEGO (24) software for understanding the distribution of gene functions from the macro level. DEGs for KEGG enrichment analysis were mapped to the KEGG database. After the GO and KEGG data were analyzed, a P-value was obtained. The protein-protein interaction network of the top 20 DEGs was completed based on the local database of BGI which integrated the Biomolecular Interaction Network Database (BIND), Biological General Repository for Interaction Datasets (BioGRID) and the Human Protein Reference Database (HPRD).
Verification
ELK2P was tagged for the member of the ETS oncogene family, and the expression of bone gamma-carboxyglutamate protein (BGLAP) had been tested in our previous study (25). Hence, we verified the top 6 DEGs excluding ELK2P and BGLAP by RT-qPCR after the screening of the DEGs. Reverse transcription for cDNA synthesis was performed as per the manufacturer's instructions (RevertAid First Strand cDNA synthesis kit). Glyceraldehyde 3-phosphate dehydrogenase (GAPDH) was used as the reference gene, and qPCR was performed as per the protocol recommended by the manufacturer (FastStart Universal SYBR-Green Master) with the primers listed in Table I. The PCR program was set as follows: pre-denaturation at 95°C for 10 min, cooling to 60°C for 15 sec, held for 60 sec for 40 cycles, elongation at 60°C for 5 min, and then an increase from 75°C to 95°C at 1°C/20 sec to obtain the melting curve. Relative gene expression was normalized to the control group by referring to the Livak and Schmittgen method (26).
Statistical analysis
The NOISeq method was used to assess the DEGs in the 2 groups using the following criteria as default: fold-change ≥2.0 and diverge probability ≥0.8. The associated P-value for the GO and KEGG enrichment analysis was subjected to Bonferroni Correction and a corrected P-value obtained. A Corrected P-value <0.05 was considered indicative of a statistically significant difference. Statistical analysis of the verification was performed using SPSS 20.0 (International Business Machines Corporation, Armonk, NY, USA). An independent sample t-test was used to analyze the RT-qPCR verification results between the ossified and non-ossified degeneration group. Inter-group differences with an associated P-value <0.05 were considered as statistically significant.
Results
Sequencing and raw data filtering
We sequenced 12 samples by RNA-Seq technology and generated an average of 13,127,595.5 raw reads. After the low quantity reads were filtered, an average of 12,994,408 clean reads was acquired. The data pertaining to each sample are listed in Table II. The numbers of expressed genes and the expression ratio in each specimen are shown in Fig. 2.
Differentially expressed genes
A total of 132 DEGs was detected in this study. A total of 129 in the ossification group were upregulated as compared to those in the control group, while 3 genes were found to be downregulated. The top 20 DEGs according to the probability value are listed in Table III.
GO enrichment analysis
The GO functional classification used in this study is illustrated in Fig. 3. Referring to the corrected P-values, the top 3 enrichment terms in all the 3 GO ontologies were listed: cellular component, molecular function and biological process (Table IV).
KEGG enrichment analysis
The most enriched pathways based on KEGG were asthma, malaria, African trypanosomiasis, staphylococcus aureus infection, collecting duct acid secretion and primary immunodeficiency. The top 10 KEGG enriched pathways are listed in Table V.
Protein-protein interaction network analysis
Protein-protein interaction network analysis was performed using the coding proteins of the top 20 DEGs based on the interaction network database BIND; the results are illustrated using Medusa software (27) (Fig. 4).
Verification by RT-qPCR
The results of RT-PCR are shown in Fig. 5. The ossified discs exhibierd a higher expression of the top 6 DEGs: SOST, IGJ, DEFA4, SFRP4, PRTN3 and CTSG, with associated P-values of 0.045, 0.000, 0.008, 0.010, 0.015 and 0.002, respectively, and F-values of 4.868, 159.719, 10.201, 8.401, 8.266 and 17.783, respectively, calculated using an independent sample t-test. The results were consistent with those of the DEG sequencing.
Discussion
Disc ossification is not a rare phenomenon (28). The ossification rate is known to correlate with the degeneration grade of intervertebral discs and the Modic type endplate (19). The underlying mechanisms responsible for the ossification of intervertebral discs are unclear. To the best of our knowledge, the present study is the first to determine the DEGs between ossified discs and non-ossified discs by employing high-throughput sequencing technology. After sequencing, the NOISeq method was applied to determine the DEGs. These results were verified and were found to be consistent with those of DEGs obtained by performing RT-qPCR.
We selected the top 20 DEGs and focused on some of these genes as our candidate genes for further research. According to the diverge probability calculated by the NOISeq method, SOST was the topmost DEG. The coding product of SOST is sclerostin, a glycoprotein secreted by osteocytes, which is an important marker of late-stage osteoblasts/osteocytes. An earlier study suggested sclerostin to act as a BMP-inhibitor by competitively binding to BMP receptors, thereby inhibiting osteoblast differentiation and osteogenesis induced by BMPs (29). By contrast, other studies have suggested that the inhibitory effect of sclerostin was mediated through the Wnt signaling pathway and not by BMP (30,31). Since there are contradicting results related to the effect of sclerostin on the ossification of discs (32), further studies are required to delineate its role.
We were interested in another DEG, BGLAP, which codes for osteocalcin protein. Osteocalcin is an important protein marker of late stage osteogenesis secreted by osteoblasts (33). It is abundant in γ-carboxyglutamic acid (Gla) residue and can combine with the calcium ions in hydroxyapatite to form osteoid element (34,35). The upregulated expression of SOST and BGLAP indicated that the calcification of the intervertebral disc was not a dystrophic calcification, but an osteogenic process involving the participation of osteoblasts. This observation is consistent with some earlier reports (14,15), and it also serves to verify the rationale behind our sampling process using CT and micro CT scan.
All the top 3 cellular components of GO ontologies analysis based on DEGs were extracellular components which indicate that these proteins are of the secretory type. GO functions are mainly related to the interaction of the glycoprotein in the cell membrane with the extracellular matrix. The GO processes are associated with completing response to stimulus, immune reflex and defense.
The top 5 KEGG enrichment pathways are associated with infection and inflammation. On scrutinizing the top 20 DEGs, we found that 3 genes: SOST, WIF1 and SFRP4 were related to the inhibition of the Wnt signaling pathway (36–38). SOST, the topmost DEG, is not termed by KEGG. The KEGG enrichment pathway analysis was performed using all the DEGs, which, to some extent, could have served to obscure the identification of the Wnt pathway. We believe that disc ossification is associated with the inhibition of the Wnt pathway due to the identification of several gene products related to Wnt inhibition. Whether the ossification is caused by Wnt pathway inhibition or Wnt was inhibited by the ossification is still not clear from the findings of our study.
Degenerative discs do exhibit the expression of BMP-2 (39), which is a widely known inducer of bone formation. Furthermore, the existence of crosstalk between BMP-2 and Wnt pathways in the degenerated discs suggests that BMP-2 pathway can inhibit the Wnt pathway (40). However in a recent study, Wnt/β-catenin signaling was found to activate BMP-2 expression in osteoblasts (41). Taking into consideration the association between the Wnt/β-catenin signaling and BMP-2 pathway, we hypothesized that the BMP-2 pathway may play a role in the the overexpression of Wnt inhibitors. This needs to be further investigated in in vitro cell experiments.
Although inflammation is thought to play an important role in intervertebral disc degeneration and herniation (18), definitive evidence of this association has yet to be obtained. Several inflammatory factors, such as tumor necrosis factor (TNF)-α, interleukin (IL)-1β, IL-6 and IL-8, have been implicated in the causation of disc degeneration and herniation (42–45). By contrast, this study found that none of the DEGs encoded for these inflammatory factors. We found a comprehensive expression of inflammation-related factors, such as defencin α4, proteinase 3, cathepsin G, elastase, azurocidin 1 and defensin α1B, which are encoded by the genes, DEFA4, PRTN3, CTSG, ELANE, AZU1 and DEFA1B. These genes were among the top 20 differentially expressed genes between the experimental and the control groups. As regards the link between the intervertebral disc calcification and inflammation around the herniated discs (19), as well as the association between inflammation and ossification (46), we believe that the inflammatory factors coding DEGs may be specifically related to disc ossification.
The main limitation of the present study is its descriptive nature. The abnormal Wnt pathway and the overexpression of some inflammation-related factors documented in this study appear to be two crucial elements involved in causing disc ossification. However, in vivo and in vitro studies are required to thoroughly assess this observed association. We performed an elementary study of the molecular mechanisms underlying disc ossification. We found a potential mechanism, which, however, needs to be further investigated as a step towards the development of gene therapy targeting disc ossification in the future. We believe that our study is helpful to finally elucidate out the mechanisms responsible for disc ossification, and patients may benefit from the prevention and treatment of herniated disc with ossification at the gene level.