Sequencing and bioinformatics analysis of the differentially expressed genes in herniated discs with or without calcification.

The purpose of this study was to detect the differentially expressed genes between ossified herniated discs and herniated discs without ossification. In addition, we sought to identify a few candidate genes and pathways by using bioinformatics analysis. We analyzed 6 samples each of ossified herniated discs (experimental group) and herniated discs without ossification (control group). Purified mRNA and cDNA extracted from the samples were subjected to sequencing. The NOISeq method was used to statistically identify the differentially expressed genes (DEGs) between the 2 groups. An in-depth analysis using bioinformatics tools based on the DEGs was performed using Gene Ontology (GO) enrichment, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment, and protein-protein interaction network analysis. The top 6 DEGs were verified using reverse transcription-quantitative polymerase chain reaction (RT-qPCR). A total of 132 DEGs was detected. A total of 129 genes in the ossified group were upregulated and 3 genes were found to be downregulated as compared to the control group. The top 3 cellular components in GO ontologies analysis were extracellular matrix components. GO functions were mainly related to the glycoprotein in the cell membrane and extracellular matrix. The GO process was related to completing response to stimulus, immune reflex and defense. The top 5 KEGG enrichment pathways were associated with infection and inflammation. Three of the top 20 DEGs [sclerostin (SOST), WNT inhibitory factor 1 (WIF1) and secreted frizzled related protein 4 (SFRP4)] were related to the inhibition of the Wnt pathway. The ossified discs exhibited a higher expression of the top 6 DEGs [SOST, joining chain of multimeric IgA and IgM (IGJ; also known as JCHAIN), defensin alpha 4 (DEFA4), SFRP4, proteinase 3 (PRTN3) and cathepsin G (CTSG)], with the associated P-values of 0.045, 0.000, 0.008, 0.010, 0.015 and 0.002, respectively, as calculated by the independent sample t-test. The gene expression profiling of the 2 groups revealed differential gene expression. Thus, our data suggest that Wnt pathway abnormality and local inflammation may be related to disc ossification.

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.

Figure 1

CT and micro CT imaging of ossified disc. (A) CT scan showing the ossified spot in the herniated disc; (B) micro CT screening confirmed the ossification. CT, computed tomography.

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).

Table I

Sequences of primers for RT-qPCR.

Gene	Upstream primer	Downstream primer
GAPDH	AGCCACATCGCTCAGACAC	GCCCAATACGACCAAATCC
SOST	GTGGCAGGCGTTCAAGAATG	CCCGGTTCATGGTCTTGTTG
IGJ	GGAGTCCTGGCGGTTTTTAT	GGATCTTCGGAAGAACGGAT
DEFA4	GCCCTCCTCGCTGCTATTCT	ATGTCCTGGTCTTCTGGCCC
SFRP4	ACTGCGAGCCCCTCATGAAG	CTTCAGGCGAGATGCACACG
PRTN3	CCTGCAGGAGCTCAATGTCA	GAGTCTCCGAAGCAGATGCC
CTSG	AACGGAAGGCTGCCTTCAAG	CTGGAGGAACCCCTGACGAC

PCR, polymerase chain reaction; GAPDH, glyceraldehyde 3-phosphate dehydrogenase.

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.

Table II

Sequencing data for each specimen.

Sample	Age (years)	Gender	Raw reads no.	Clean reads no.	Clean datarate (%)	Clean Read1 Q20 (%) ≥95	Clean Read1 Q30 (%) ≥90	Clean Reads ≥10 M	Gene mapping ratio (%) ≥80	Genome mapping ratio (%) ≥50
Ossified 1	62	F	13,127,663	13,113,920	99.89	98.9 (Y)	96.0 (Y)	13.11M (Y)	92.22 (Y)	86.64 (Y)
Ossified 2	33	M	13,127,599	13,100,270	99.79	98.8 (Y)	95.9 (Y)	13.10M (Y)	92.03 (Y)	86.20 (Y)
Ossified 3	62	F	13,127,628	13,106,545	99.83	98.9 (Y)	96.0 (Y)	13.11M (Y)	90.67 (Y)	88.20 (Y)
Ossified 4	63	F	13,127,584	12,524,724	95.4	98.9 (Y)	96.1 (Y)	12.52M (Y)	83.85 (Y)	88.96 (Y)
Ossified 5	50	M	13,127,257	13,110,977	99.87	99.1 (Y)	97.0 (Y)	13.11M (Y)	91.05 (Y)	88.56 (Y)
Ossified 6	26	M	13,127,394	12,597,800	95.96	99.0 (Y)	96.6 (Y)	12.60M (Y)	82.21 (Y)	88.91 (Y)
Control 1	58	F	13,127,248	13,107,152	99.84	98.8 (Y)	95.9 (Y)	13.11M (Y)	92.30 (Y)	86.93 (Y)
Control 2	12	F	13,127,788	12,907,909	98.32	98.8 (Y)	95.9 (Y)	12.91M (Y)	92.56 (Y)	85.80 (Y)
Control 3	53	F	13,127,674	13,059,401	99.47	98.8 (Y)	95.9 (Y)	13.06M (Y)	93.05 (Y)	86.58 (Y)
Control 4	25	M	13,127,590	13,079,633	99.63	99.1 (Y)	97.0 (Y)	13.08M (Y)	91.84 (Y)	86.71 (Y)
Control 5	24	F	13,127,758	13,114,894	99.9	99.0 (Y)	96.4 (Y)	13.11M (Y)	92.16 (Y)	86.40 (Y)
Control 6	38	M	13,127,963	13,109,683	99.86	99.2 (Y)	97.3 (Y)	13.11M (Y)	93.13 (Y)	82.50 (Y)

Mean age of ossified group, 49.3 years; mean age of control group, 35 years. (Y) indicates that the number passed the quality control criteria.

Figure 2

Gene expression in each specimen. With a total gene number of 24,048 in the database, the numbers of expressed genes in our study ranged from 17,026 in 'ossified 6' to 18,232 in 'ossified 4'; the expression ratios were from 70.8 to 75.82%.

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.

Table III

Top 20 DEGs.

Gene ID	Symbol	log2Ratio(ossified/control)	Probability	Up/downregulation	Description
50964	SOST	9.658888588	0.980396	Up	Sclerostin
3512	IGJ	6.673912182	0.967813	Up	Immunoglobulin J polypeptide
632	BGLAP	5.492969178	0.962151	Up	Bone γ-carboxyglutamate (gla) protein
2003	ELK2AP	5.437789846	0.960819	Up	ELK2A
1669	DEFA4	5.317314168	0.956617	Up	Defensin, α4
6424	SFRP4	5.055694219	0.95017	Up	Secreted frizzled-related protein 4
5657	PRTN3	4.933061142	0.947346	Up	Proteinase 3
1511	CTSG	4.657946704	0.941332	Up	Cathepsin G
1991	ELANE	4.591945637	0.941047	Up	Elastase
566	AZU1	4.584099727	0.940888	Up	Azurocidin 1
100423062	IGLL5	4.310197783	0.933947	Up	Immunoglobulin λ-like polypeptide 5
2812	GP1BB	10.73075316	0.931232	Up	Glycoprotein Ib (platelet), β polypeptide
63924	CIDEC	6.76165846	0.930184	Up	Cell death-inducing DFFA-like effector c
3543	IGLL1	5.1749338	0.928683	Up	Immunoglobulin λ-like polypeptide 1
728358	DEFA1B	3.988886669	0.92734	Up	Defensin, α1B
4353	MPO	4.209517648	0.927108	Up	Myeloperoxidase
55363	HEMGN	4.89235857	0.927006	Up	Hemogen
6037	RNASE3	4.64090047	0.926718	Up	Ribonuclease, RNase A family, 3
6036	RNASE2	4.154822975	0.919141	Up	Ribonuclease, RNase A family, 2
11197	WIF1	5.058808822	0.917018	Up	WNT inhibitory factor 1

DEGs, differentially expressed genes.

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).

Figure 3

GO functional classification on DEGs. The x-axis represents the number of DEGs available at the top of the column. The y-axis represents GO terms. All GO terms are grouped into 3 ontologies: biological process, cellular component and molecular function; blue indicates biological process, brown indicates cellular component and orange indicates molecular function. DEGs, differentially expressed genes; GO, Gene Ontology.

Table IV

Top 3 enrichment terms in all three GO ontologies.

GO ontologies	GO terms	Cluster frequency	Genome frequency of use	Corrected P-value
Cellular component	GO:0044421 - extracellular region part	32 (101), 31.7%	1055 (16090), 6.6%	1.70E-12
Cellular component	GO:0005576 - extracellular region	32 (101), 31.7%	1077 (16090), 6.7%	3.02E-12
Cellular component	GO:0016023 - cytoplasmic membrane- bounded vesicle	12 (101), 11.9%	609 (16090), 3.8%	0.03246
Molecular function	GO:0030247 - polysaccharide binding	11 (95), 11.6%	159 (15165), 1.0%	3.48E-07
Molecular function	GO:0001871 - pattern binding	11 (95), 11.6%	169 (15165), 1.1%	6.60E-07
Molecular function	GO:0005539 - glycosaminoglycan binding	10 (95), 10.5%	142 (15165), 0.9%	1.54E-06
Biological process	GO:0006952 - defense response	25 (108), 23.1%	601 (14596), 4.1%	6.26E-10
Biological process	GO:0002376 - immune system process	33 (108), 30.6%	1106 (14596), 7.6%	8.06E-10
Biological process	GO:0050896 - response to stimulus	54 (108), 50.0%	3380 (14596), 23.2%	4.94E-07

GO, Gene Ontology.

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.

Table V

The top 10 KEGG enriched pathways.

Pathway ID	Pathway	DEGs with pathway annotation(108)	All genes with pathway annotation(17252)	P-value
ko05310	Asthma	7 (6.48%)	47 (0.27%)	1.58E-08
ko05144	Malaria	7 (6.48%)	76 (0.44%)	4.75E-07
ko05143	African trypanosomiasis	6 (5.56%)	54 (0.31%)	1.06E-06
ko05150	Staphylococcus aureus infection	6 (5.56%)	139 (0.81%)	0.00024
ko04966	Collecting duct acid secretion	4 (3.7%)	51 (0.3%)	0.00029
ko05340	Primary immunodeficiency	4 (3.7%)	61 (0.35%)	0.00058
ko04610	Complement and coagulation cascades	6 (5.56%)	181 (1.05%)	0.00097
ko05322	Systemic lupus erythematosus	6 (5.56%)	188 (1.09%)	0.0012
ko04640	Hematopoietic cell lineage	5 (4.63%)	152 (0.88%)	0.0027
ko03320	PPAR signaling pathway	4 (3.7%)	135 (0.78%)	0.01

KEGG, Kyoto Encyclopedia of Genes and Genomes; DEGs, differentially expressed genes.

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).

Figure 4

Protein-protein interaction network of the coding proteins of the top 20 DEGs. DEGs, differentially expressed genes.

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.

Figure 5

Verification of top 6 DEGs. The x-axis represents the top 6 DEGs; each DEG contains the comparison of the ossified and control groups. The y-axis represents the relative mRNA expression of the DEG in each group normalized to the ossified group. *P<0.05 and **P<0.01 as compared to the control group as shown by an independent sample t-test. DEGs, differentially expressed genes.

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.