This article is mentioned in:

Introduction

Intervertebral disc (IVD) disease is a common surgical disease, and IVD degeneration is hypothesized to be the first step in degenerative spinal changes (1). IVD degeneration is typically associated with disc herniation and back pain, which has a marked effect on the sufferer's life and causes severe chronic pain (2). Furthermore, lower back pain may limit the activity of individuals <45 years, which has a marked socio-economic impact (3). The etiology of lower back pain is unclear, but in 40% of cases it is associated with IVD degeneration (4). It is reported that ~40% of individuals <30 years of age, and >90% of individuals >55 years, exhibit moderate-to-severe levels of IVD degeneration (5,6). It is estimated that the costs associated with lumbar disc and lower back disorders exceed $100 billion/year in the USA (7). Traditionally, inflammation has mostly been observed as detrimental and is associated with disease progression. Although the etiology remains unknown, the osmolality has been identified to be associated with IVD disease (8–10). Katz (7) reported that the osmotic environment exerted an appreciable effect on gene expression and also affected responses to mechanical stimuli. Mavrogonatou and Kletsas (11) indicated that high osmolality was able to activate the G1 and G2 cell cycle checkpoints and affect the DNA integrity of nucleus pulposus IVD cells, triggering an enhanced DNA repair response. Although changes in the extracellular osmolality markedly influenced the behavior of IVD cells, their response to this condition has not yet been fully elucidated. In the present study, the associated biological processes and potential biomarkers in the response of IVD to osmotic stimuli were investigated using gene expression analysis.

Materials and methods

Gene expression profile

The gene expression dataset GSE1648 (8) was downloaded from the Gene Expression Omnibus (http://www.ncbi.nlm.nih.gov/geo/) database. There was a total of 11 IVD cell samples in GSE1648, including 4 hyperosmotic stimuli samples, 3 hypoosmotic stimuli samples and 4 isosmotic stimuli samples.

Data pre-processing and identification of differentially expressed genes (DEGs)

For the expression profile, the original data were converted into a recognizable format using R software version 3.1.1 (from bioconductor.org/packages/release/bioc/html/biomaRt.html), and the Affy package (12) version 1.48.0 (from bioconductor.org/packages/release/bioc/html/affy.html) was used for the background correction and normalization. The DEGs in hyperosmotic (designated DEGs-hyper) or hypoosmotic (DEGs-hypo) IVD cells, compared with isosmotic cells, were identified using the limma package (13) version 3.18.13 (from bioconductor.org/packages/release/bioc/html/limma.html) based on gene expression differences of P<0.05 and |log2(fold-change)|>0.05.

Gene Ontology (GO) term enrichment analysis

Functional annotation of DEGs is a necessary and critical step in the analysis of microarray data. The Database for Annotation, Visualization and Integrated Discovery (DAVID, http://david.abcc.ncifcrf.gov/) (14) is a tool providing a comprehensive set of functional annotation. Enriched GO terms from DEGs-hyper were identified using DAVID with a threshold of P<0.05.

Construction of the protein-protein interaction (PPI) network

The Human Protein Reference Database (HPRD) was used as a centralized platform to visually depict and integrate information pertaining to domain architecture, post-translational modifications, interaction networks and disease association for each protein in the human proteome (15). PPI network data for DEGs-hyper was obtained from HPRD and visualized using Cytoscape (16) software.

Establishment and analysis of the microRNA (miRNA)-regulated-gene network

The software TargetScan (17) was used to predict biological targets of miRNAs by searching for the presence of conserved 8-mer and 7-mer sites that matched the seed region of each miRNA. TargetScan screened for miRNAs associated with the regulation of DEGs-hyper, and miRNA-gene connections were obtained. An miRNA-regulated-gene network was established and visualized using Cytoscape software, and every node was analyzed according to the number of connections with other nodes (degree). miRNAs for which the degree was ≥1 were removed from the analysis.

Results

DEGs

A total of 43 DEGs (34 up- and 9 downregulated) were identified in DEGs-hyper, and the cluster graph is presented in Fig. 1. However, only 9 DEGs were obtained in DEGs-hypo, which were not studied further.

Figure 1. Cluster graph of differentially expressed genes in intervertebral disc cells of hyperosmotic stimuli samples compared with the isosmotic stimuli.

GO terms of DEGs-hyper

In total, 41 GO terms were enriched in DEGs-hyper, and the 10 most significantly enriched GO terms (e.g., regulation of apoptosis, regulation of programmed cell death and regulation of cell death) are presented in Table I.

Table I. Top 10 most significantly enriched GO terms of differentially expressed genes in intervertebral disc cells of hyperosmotic stimuli samples compared with the isosmotic stimuli.

Table I.

Top 10 most significantly enriched GO terms of differentially expressed genes in intervertebral disc cells of hyperosmotic stimuli samples compared with the isosmotic stimuli.

Category	GO ID	GO name	Number of genes	P-value
BP	GO:0042981	Regulation of apoptosis	9	3.48×10−4
BP	GO:0043067	Regulation of programmed cell death	9	3.73×10−4
BP	GO:0010941	Regulation of cell death	9	3.82×10−4
BP	GO:0006913	Nucleocytoplasmic transport	5	4.20×10−4
BP	GO:0051169	Nuclear transport	5	4.41×10−4
BP	GO:0012501	Programmed cell death	7	0.00232
BP	GO:0008219	Cell death	7	0.00519
BP	GO:0030183	B cell differentiation	3	0.00537
BP	GO:0016265	Death	7	0.00537
BP	GO:0051168	Nuclear export	3	0.00828

[i] GO, Gene Ontology; ID, identifier; BP, biological process.

PPI network of DEGs-hyper

The PPI network was constructed and is presented in Fig. 2, and included 376 pairs and 382 nodes. The top 20 nodes (e.g., NCOA3, SOS1 and XPO1), as defined by highest degree, are presented in Table II.

Figure 2. Protein-protein interaction network of differentially expressed genes in intervertebral disc cells of hyperosmotic stimuli samples compared with the isosmotic stimuli.

Table II. Top 20 nodes with a higher degree in the protein-protein interaction network.

Table II.

Top 20 nodes with a higher degree in the protein-protein interaction network.

Gene	Degree
NCOA3	49
SOS1	45
XPO1	45
RRAS2	27
TP53BP2	23
BCL6	22
NXF1	22
EFNB2	21
SGK1	16
ARNTL	13
LSM6	10
PTPRE	10
VCAM1	10
BDNF	9
GPSM2	7
FNBP4	6
GLUL	5
MALT1	4
MAPK1	4
PPM1D	4

[i] Degree indicates the number of connections of each node with other nodes in the network.

miRNA-regulated-gene network of DEGs-hyper

In total, 1,314 miRNA-regulated-gene connections were identified with TargetScan. A miRNA-gene-regulated network was established, including 1,314 connections and 422 nodes (Fig. 3). The top 20 nodes (e.g., ZBTB18, EFNB2 and SOBP) in the regulated network, as determined by degree, are presented in Table III.

Figure 3. microRNA-gene-regulated network of differentially expressed genes in intervertebral disc cells of hyperosmotic stimuli samples compared with the isosmotic stimuli.

Table III. Top 20 nodes with higher degree in the microRNA-gene-regulated network.

Table III.

Top 20 nodes with higher degree in the microRNA-gene-regulated network.

Gene	Degree
ZBTB18	81
EFNB2	73
SOBP	72
NXF1	68
BDNF	66
PTBP3	64
PPM1D	59
BAMBI	58
NCOA3	55
KLHL24	47
SGK1	45
GAREM	42
NUAK1	42
SOS1	42
BCL6	38
XPO1	37
FAM110B	34
TNFAIP8	34
KLHL9	29

[i] Degree indicates the number of connections of each node with other nodes in the network.

Discussion

The normal upper limit of water in the human IVD is ~500 mOsm/kg, which is higher than that routinely encountered in the majority of other parts of the body (~300 mOsm/kg) (18). It was previously reported that ~18% of the fluid of IVD was lost and re-imbibed during a diurnal cycle, with consequent changes in osmolality (19). Changes in osmolality are an important component of the physicochemical environment of IVD, as variations in disc loading leads to alterations of disc hydration. A study of IVD cells in a three-dimensional alginate culture system confirmed that the biological response to altered osmolarity is mediated, in part, by changes at the transcriptional level (20). In the present study, 43 DEGs were identified in hyperosmotic cells, whereas 9 DEGs were identified in hypoosmotic cells, compared with isosmotic cells. These results suggest that IVD is more sensitive to hyperosmotic stimuli than to hypoosmotic stimuli, and that the effects of hyperosmotic stimuli may be far greater at the transcriptional level. Hyperosmotic stimuli were previously demonstrated to elicit calcium transience in IVD cells that were modulated by the stability of the actin cytoskeleton (21,22). In the present study, DEGs-hyper were the basis of further research, whereas DEGs-hypo were not further studied.

The enriched GO terms for DEGs-hyper were predominantly associated with the biological processes of apoptosis and cell death, and the regulation of these processes. Three subfamilies were identified: Extracellular-signal-regulated kinase (ERK1/2), p38 mitogen-activated protein kinase (p38 MAPK, p38) and c-Jun N-terminal kinase (JNK1/2). Nucleus pulposus cell apoptosis is one of the changes associated with IVD degeneration (23). p38 and JNK1/2 are serine/threonine protein kinases activated by various stress stimuli, including osmotic shock, toxic compounds and pro-inflammatory cytokines (24). ERK1/2 was previously demonstrated to be activated by high osmolality in nucleus pulposus cells in vitro (25). Dong et al (26) reported that high osmolality activated p38 MAPK, JNK1/2 and ERK1/2 in rabbit nucleus pulposus cells. The activated p38 MAPK and JNK1/2 induced cell apoptosis; by contrast, the activation of ERK1/2 promoted cell survival. A recent study (27) indicated that the effects of osmolality on nucleus pulposus cell apoptosis depended on the osmolality level (hypo-, iso- or hyper-) and osmolality mode (constant or cyclic). Furthermore, inhibition of the ERK1/2 pathway promoted nucleus pulposus cell apoptosis in this process (27). Therefore, it was suspected that the biological processes of apoptosis and cell death may be associated with the effect of hyperosmolality on IVD diseases.

The PPI network was analyzed, and nuclear receptor coactivator 3 (NCOA3), SOS Ras/Rac guanine nucleotide exchange factor 1 (SOS1) and exportin 1 (XPO1) were the top three nodes as defined by degree. NCOA3 was identified to be associated with the pathophysiology of osteoarthritis (OA) (28,29). A recent study concluded that NCOA3 was subject to a cis-acting expression quantitative trait locus in articular cartilage, which was associated with the OA association signal and associated with the OA-associated allele of the functional single-nucleotide polymorphism rs116855380 (30). Although there is no direct evidence that NCOA3 is associated with IVD degeneration, IVD degeneration and OA may exhibit similarities in their occurrence and development, as they are joint degeneration diseases. SOS1 primarily encoded membrane-bound guanine nucleotide-binding proteins in humans, which function in the transduction of signals that control cell growth and differentiation (31).

It was previously reported that SOS1 participates in the process of apoptosis (32,33). It was identified that the biological processes of apoptosis and cell death may be associated with the effect of high osmolality on IVD diseases. The mutation of SOS1 serves an important role in the salt-tolerance and osmotic stimuli of Arabidopsis and tobacco (34–37). XPO1 is a specific receptor for leucine-rich nuclear export sequences (38) and was identified to mediate the nuclear export of proteins in a range of species (39,40). Ferrigno et al (41) indicated that XPO1 mediated the export of high osmolality glycerol response MAPK to the cytoplasm of cells adapted to hyper-osmotic stimuli. Thus, NCOA3, SOS1 and XPO1 may be genes directly affected by osmotic stimuli in IVD cells, although more trials and clinical validation are required. Similarly, ZBTB18, EFNB2 and SOBP were the top 3 nodes in the miRNA-gene-regulated network, and they may exhibit an intimate association with the effects of osmotic stimuli on IVD.

The dataset GSE1648 was created by Boyd et al (8); however, the present study adopted different methods and obtained novel results compared with the Boyd et al study (8). First, differential expression analysis was performed for cells in hyperosmotic and hypoosmotic conditions in the present study, whereas Boyd et al (8) only conducted differential expression analysis of cells in hyperosmotic conditions. Secondly, the PPI network and miRNA-gene-regulated network were constructed in the present study; however, this was not carried out in the study by Boyd et al (8). Finally, potential target genes (NCOA3, SOS1, XPO1, ZBTB18, EFNB2 and SOBP) were identified from the PPI network and the miRNA-gene-regulated network, whereas they were not identified in the study by Boyd et al (8). These genes exhibited more interactions compared with those identified by Boyd et al, and therefore they were considered more reliable.

References