This article is mentioned in:

Introduction

Low back pain is mainly caused by degenerative disc disease (DDD) (1). A number of pathological factors may be responsible for DDD, such as apoptosis, inflammation and cartilage matrix degradation (2–4); metabolic exchange reduction, which results in nutrition insufficiency and waste accumulation, is considered to be the most important of these mechanisms (5). The intervertebral disc (IVD) is an avascular organ, and metabolic exchange primarily depends on the diffusion of cartilage endplate (CEP) in mature IVDs (6). The CEP is a horizontal, thin layer of hyaline cartilage, which separates the IVD from the vertebrae. Blood vessels in adjacent vertebrae end at the IVD-vertebra interface without reaching the inner component of the IVD (7). CEP is the most important channel of metabolic substance exchange, and its degeneration is considered to be the origination of DDD (8).

Our previous study reported the existence of CEP-derived stem cells (CESCs), which exhibited the ability for both chondrogenic and osteogenic differentiation (9). This differentiation potential enable CESCs to serve important roles in the restoration and regeneration of cartilage, and the differentiation status of CESCs may be responsible for CEP degeneration.

As an avascular organ, the IVD remains in a physiologically hypoxic microenvironment (10). Hypoxia has been reported to regulate the biological activities of mesenchymal stem cells (MSCs), such as apoptosis, proliferation and differentiation (11,12), which indicated that physiological hypoxia may serve vital roles in regulating the physiological functions of CESCs.

Data from our previous study on the effects of hypoxia on the biological activities of CESCs indicated that alternative splicing (AS) may serve as a ‘middle link’ between them, through which hypoxia may regulate the biological activities of CESCs (9). This assumption may be due to two reasons. First, several AS events may be hypoxia-inducible; for example, inhibitory PAS domain protein generated unique splicing variants of exon 6 and exon 3 exclusion in a hypoxic microenvironment, which subsequently formed a new negative-feedback modulation mode of adaptive responses to ischemia (13). In addition, a spliced variant of tyrosine kinase receptor A (TrkA), TrkA-III, was demonstrated to be induced by hypoxia (14), and the ability of this spliced variant to form stress-resistant and angiogenic isoforms suggested its potential roles in the protective effects against hypoxia (14). Second, the roles of AS in stem cell fate also have been reported previously. For example, the variant of transcription initiation factor TFIID subunit 4 (TAF4) that lacks the human TAF4-TAF homology domain has been demonstrated to lead to the promotion of chondrogenic differentiation in human MSCs (15). A conserved embryonic stem cell-specific AS event was previously identified that altered the DNA-binding activity of forkhead box P1, which indicated that the pluripotency and reprogramming of embryonic stem cells may be regulated (16).

High-throughput screening is a powerful technology for investigating the changes in the transcription profile and AS genome-wide. A previous study used exon microarrays to examine the gene expression profiles and AS events in human umbilical vein endothelial cells under hypoxic conditions (17). Another study used microarray-screening technology to identify a hypoxia-induced splice variant, laminin subunit α3, which was closely associated with poor prognosis in patients with head and neck cancers (18). Our previous studies revealed gene expression profiles and AS events during the process of chondrogenesis in CESCs under hypoxia and normoxia (19,20); however, high-throughput analysis of the general regulatory effects of hypoxia on the various aspects of CESCs biological function (without the induction of differentiation) has not yet been reported. In the present study, isolated CESCs were cultured in normal growth medium under normoxic and hypoxic conditions. RNA was extracted, purified and analyzed using the Affymetrix Human Transcriptome Array 2.0 System, and a comparative analysis of gene expressions and AS events between the normoxic and hypoxic groups was conducted. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis were used for the functional enrichment of genes of interest to explain the possible mechanisms of transcription and alternative splicing. Data from the present study may aid in identifying the roles of the physiological hypoxic microenvironment on the fate of CESCs, which may be beneficial to the understanding the pathological mechanisms of DDD and may lead to the development of new therapies.

Materials and methods

Patients and subjects

CEP samples were obtained from three patients that underwent discectomy operations at the Xinqiao Hospital, Third Military Medical University (Chongqing, China; Table I). All procedures in the present study were approved by the Ethics Committee of Xinqiao Hospital, Third Military Medical University and were conducted in accordance with the Declaration of Helsinki; written informed consent was obtained from each patient.

Table I. Patient information.

Table I.

Patient information.

Case	Sex	Age (years)	Diagnosis	Degenerated disc level	Surgery
1	Male	50	Spondylolisthesis	L4-L5	TLIF
2	Male	55	Spondylolisthesis	L4-L5	TLIF
3	Female	52	Spondylolisthesis	L4-L5	TLIF

[i] L, lumbar; TLIF, transforaminal lumbar interbody fusion.

CESC isolation and culture

CEP samples were cut into pieces and digested in Dulbecco's modified Eagle's medium (DMEM)/Ham's F12 medium (Hyclone; GE Healthcare Life Sciences, Logan, UT, USA) containing 1% fetal calf serum (FCS; Gibco; Thermo Fisher Scientific, Inc., Waltham, MA, USA) and 0.2% collagenase II (Sigma-Aldrich; Merck KGaA, Darmstadt, Germany) for 12 h at 37°C. The suspension was passed through a cell filter (70 µm) and centrifuged at 110 × g for 5 min at room temperature. Following removal of the supernatant, the pellet was resuspended in 3 ml DMEM/F12 containing 10% FCS and 1% penicillin-streptomycin (Hyclone; GE Healthcare Life Sciences). The cells were then cultured in a 25-cm2 cell-culture flask at 37°C and 5% CO2 for 6 weeks. Cells were subjected to agarose selection following the first passage of expansion, as previously described (9). Briefly, Costar culture dishes (Corning Inc., Corning, NY, USA) were pre-coated with 1% low melting-point agarose, and a mixture containing DMEM/F12 (0.5 ml), 2% low melting-point agarose (0.5 ml) and 5×104 suspended CEP cells in culture medium (1 ml) was transferred to the culture dishes and incubated in humidified atmosphere in 5% CO2 and 37°C. The culture medium was changed twice per week and, 6 weeks later, cell aggregates were carefully aspirated with a sterile Pasteur pipette and cultured in a 25 cm2 cell culture flask in 5% CO2 and 37°C for 2 weeks to reach 80–100% confluence. CESCs which were sub-cultured to the third passage, after ~6 weeks of culture, were used.

Oxygen deprivation

For hypoxic culture conditions, CESCs at a density of 2.5×104/ml were cultured in a 1% O2 microenvironment at 37°C. For normoxic culture conditions, CESCs were cultured in a 21% O2 microenvironment at 37°C. The culture medium was changed twice per week over a 21-day period.

GeneChip human transcriptome array (HTA) 2.0

Total RNA was extracted from cells (~5×104) using TRIzol (Invitrogen; Thermo Fisher Scientific, Inc.) according to the manufacturer's protocol. RNA was purified with an RNeasy kit (Qiagen GmbH, Hilden, Germany) according to the manufacturer's protocol. The RNA purity was determined by the ratio of optical density (OD) value in 260 nm/OD value in 280 nm using a NanoDrop spectrophotometer (NanoDrop Technologies; Thermo Fisher Scientific, Inc., Wilmington, DE, USA). RNA quality was assessed by 1% agarose gel electrophoresis. The samples with a bright band of ribosomal 28S/18S RNA at a ratio >1.5:1 were used for the microarray analysis. For gene expression profiling, RNA was hybridized using the GeneChip HTA 2.0 system (Affymetrix; Thermo Fisher Scientific, Inc.). A total of 100 ng RNA was used to generate biotinylated and amplified sense-strand cDNA from the whole expressed genome using GeneChip WT PLUS Reagent kit (Applied Biosystems; Thermo Fisher Scientific, Inc.) according to the manufacturer's protocol. GeneChip HTA 2.0 cartridges were hybridized in a 45°C incubator (GeneChip Hybridization Oven 640; Thermo Fisher Scientific, Inc.) for 16 h under rotation at 60 rpm. After the hybridization, the microarrays were stained using the GeneChip Fluidics Station 450 (Thermo Fisher Scientific, Inc.) and screened using the GeneChip Scanner 3000 7G (Thermo Fisher Scientific, Inc.) according to the manufacturer's protocol. Gene probes of the HTA 2.0 target exons and junctions and are able to detect gene expression and AS events simultaneously. Microarray screening was performed by CapitalBio Corporation (Beijing, China) according to the manufacturer's protocol. Image signals of the probes were obtained and saved as DAT files. The GeneChip Command Console (AGCC) software version 3.2 (Affymetrix; Thermo Fisher Scientific, Inc.) was used to convert the image signals (DAT files) into digital signals (CEL files), which were converted to CHP files for probeset signal integration, quantile normalization and background correction through the robust multichip analysis (RMA) algorithm using the Affymetrix Expression Console software version 1.3 (Affymetrix; Thermo Fisher Scientific, Inc.). The CHP files were then analyzed by Affymetrix Transcriptome Analysis Console software version 3.0 (Affymetrix; Thermo Fisher Scientific, Inc.) to determine the differentially expressed genes (DEGs) and alternatively spliced genes (ASGs). GO terms and related signaling pathway enrichments were determined using web tools, such as DAVID (https://david.ncifcrf.gov), Molecule Annotation System (http://bioinfo.capitalbio.com/mas3) and KEGG (http://www.genome.jp/kegg. The workflow of the present study procedures was as previously described (20). The microarray data have been submitted to the Gene Expression Omnibus (accession no. GSE84484).

Criteria for the detection of DEGs and ASGs

Data from the cell samples cultured under normoxia were used as a control to calculate the fold change in gene expression. The default filter criteria were set as fold change ≤-2 or ≥2 (linear) for significant DEGs. For ASG identification, the splicing index (SI) calculation mode was used to analyze the level of exon inclusion/exclusion. SI indicates the ratio of the signal intensity of an exon normalized to that of the target gene between two groups. SI values were calculated using the following formulas:

Where NI(i,j)A represents the normalized intensity (NI) of the ith exon signal normalized to that of the jth gene signal under condition A, N represents the normoxic culture condition and H represents the hypoxic culture condition. The default filter criteria were defined as SI (linear) ≤-2 or ≥2.

DEG validation by reverse transcription-quantitative polymerase chain reaction (RT-qPCR)

Total RNA was extracted from cells (~5×104 cells/ml) using TRIzol reagent (Invitrogen; Thermo Fisher Scientific, Inc.) according to the manufacturer's protocol. RNA was purified with an RNeasy kit (Qiagen GmbH, Hilden, Germany) according to the manufacturer's protocol. The RNA purity was determined by the ratio of optical density (OD) value in 260 nm/OD value in 280 nm using a NanoDrop spectrophotometer (NanoDrop Technologies; Thermo Fisher Scientific, Inc.). RNA quality was assessed by 1% agarose gel electrophoresis. The samples with a bright band of ribosomal 28S/18S RNA at a ratio >1.5:1 were used. RNA was reverse transcribed into cDNA with the PrimeScript RT Master Mix kit (cat no. RR047A; Takara Bio, Inc., Otsu Japan). qPCR was performed with the ABI 7900H Real-Time PCR system and the SYBR Premix Ex Taq II kit (cat no. RR820A; Takara Bio, Inc.) with the following conditions: At 95°C for 30 sec, then 40 cycles at 95°C for 5 sec and at 60°C for 34 sec; a dissociation curve analysis was performed in a temperature range of 60–95°C. The mRNA expression levels of each gene were normalized to the level of β-actin expression and analyzed using the 2−∆∆Cq method (21). Data were analyzed using StepOne software version 2.3 (Applied Biosystems; Thermo Fisher Scientific, Inc.). The primers used for qPCR are listed in Table II-A.

Table II. List of primers used for quantitative and semi-quantitative PCR.

Table II.

List of primers used for quantitative and semi-quantitative PCR.

A, Quantitative PCR primers for DEG validation
Gene	Primer sequence (5′→3′)
VCAM1	F: GGACCACATCTACGCTGACA
VCAM1	R: TTGACTGTGATCGGCTTCCC
COMP	F: CTTCAGGGCCTTCCAGACAG
COMP	R: CGCATAGTCGTCATCCGTGA
FGF9	F: CGGCTACAACGCTCCGC
FGF9	R: TGGCGGCGACAAATCTCC
IRAK3	F: AGGATTTCCGCGGTTGTGTA
IRAK3	R: ACTCAACACTGCTCCCAGGC
CFD	F: GCGGTGAGGAGGCCTGG
CFD	R: GAACCTGCACCTTCCCGTC
TFRC	F: CCAGGCTATAAACCGCCG
TFRC	R: CCAGGCTGAACCGGGTATATG
PGAM2	F: TAGCAAGGAGCGTCGGTACG
PGAM2	R: GCCTGGTCTGACATCCCTTC
CAPN6	F: CTTGGGGCAGGTTACTCTGG
CAPN6	R: CTAAGACCTGGCACCTGACG
β-actin	F: CAACCGGGAAGGAAATGAATGG
β-actin	R: GCCCAATACGACCAAATCAGAG
B, Semi-quantitative PCR primers for ASGs validation
Gene	Primer sequence (5′→3′)
GFRA2	F: GGACTCGGGATCTTCATCG
GFRA2	R: AGGTCACCGGCTCATAGGG
CADM1	F: AGTCGATGATGAAATGCC
CADM1	R: CAGAATGATGAGCAAGCA
NR4A3	F: CCGCTCCTCCTACACTCT
NR4A3	R: TCCGTGGTGTATTCCGAG
RXFP1	F: AGTGCTCCCTTGGCTATT
RXFP1	R: AGAAACCGATGGAACAGC
NDUFC1	F: AACGAGAACAACGGAGGC
NDUFC1	R: AGTTGGCAACAGAACCAG
TMEM106B	F: GTTATCCTACCCCTCCCC
TMEM106B	R: CCAGTCCATTCCTCATGT
FN1	F: TTATAGAATTACCACAACCC
FN1	R: TTCACAGGTGAGTAACGC
β-actin	F: CAACCGGGAAGGAAATGAATGG
β-actin	R: GCCCAATACGACCAAATCAGAG

[i] ASG, alternatively spliced gene; CADM1, cell adhesion molecule 1; CAPN6, calpain 6; CFD, complement factor D; COMP, cartilage oligomeric matrix protein; DEG, differentially expressed gene; F, forward; FGF9, fibroblast growth factor 9; FN1, fibronectin 1; GFRA2, GDNF family receptor α2; IRAK3, interleukin 1 receptor-associated kinase 3; NDUFC1, NADH:ubiquinione oxidoreductase subunit C1; NR4A3, nuclear receptor subfamily 4 group A member 3; PCR, polymerase chain reaction; PGAM2, phosphoglycerate mutase 2; R, reverse; RXFP1, relaxin/insulin-like family peptide receptor 1; TFRC, transferrin receptor; TMEM106B, transmembrane protein 106B; VCAM1, vascular cell adhesion molecule 1.

ASG validation by semi-quantitative RT-PCR

Total RNA extraction and reverse transcription of cDNA were carried out as aforementioned. For semi-quantitative PCR, cDNA (2 µl) was combined with Premix Taq (cat no. RR901A; Takara Bio, Inc.) and primers that were designed flanking the constitutively expressed exons; primers are listed in Table II-B. PCR thermocycling conditions were as follows: At 94°C for 30 sec, then 30 cycles at 94°C for 30 sec, at 55°C for 30 sec and at 72°C for 60 sec. PCR products were resolved on a 1.5% agarose gel, using a DNA ladder (Thermo Fisher Scientific, Inc.) and gold view nucleic acid gel stain regent (Beijing Solarbio Science & Technology Co., Ltd., Beijing, China). Gels were imaged using a fluorescence imaging analysis system (Bio-Rad Laboratories, Inc., Hercules, CA, USA). The mRNA expression levels of each ASG were normalized to β-actin. Differences of ASG between two groups were reflected in the ratio of exon-inclusion fragment expression level vs. β-actin expression level/exon-exclusion fragment expression level vs. β-actin expression level. The ASGs used for validation were in accordance with the following selection criteria: i) The whole exon exclusion/inclusion was included; ii) relatively higher SI absolute value was included; iii) the first or the last AS exon was excluded.

Statistical analysis

The data of 3 independent experiments were expressed as the mean ± standard deviation. Independent-samples t-test was used to determine the significance between groups. P<0.05 was considered to indicate a statistically significant difference. Statistical analysis was performed using SPSS software version 19.0 (IBM Corp., Armonk, NY, USA).

Results

Detection, validation and functional analysis of DEGs under normoxic and hypoxic conditions in CESCs

HTA 2.0 data analysis identified a total of 448 DEGs, of which 330 (73%) were upregulated and 118 (27%) were downregulated. Based on these results, 10 DEGs were chosen for validation by RT-qPCR, and 8 of these were validated successfully. Since the GeneChip signal increases linearly as target gene expression increases, but the qPCR signal increases exponentially as target gene expression increases, the ratio calculated by GeneChip and qPCR was different. Thus, a similar tendency for increase or decrease in the results from both methods was considered to indicate successful validation of the DEGs. In the present study, the expression of 5 genes [vascular cell adhesion molecule 1, cartilage oligomeric matrix protein, fibroblast growth factor 9, interleukin 1 receptor-associated kinase 3 (IRAK3) and complement factor D (CFD)] was upregulated in the hypoxia group compared with the normoxia group. Conversely, the expression of 3 genes (transferrin receptor, phosphoglycerate mutase 2 and calpain 6) was downregulated in the hypoxia group compared with the normoxia group (Fig. 1).

Figure 1. Validation of the DEGs identified for CESCs under normoxic and hypoxic conditions. CESCs were cultured under normoxia or hypoxia for 21 days and 10 of the identified DEGs were selected for further validation. A total of 8 of the 10 DEGs were successfully validated by reverse transcription-quantitative polymerase chain reaction as either (A) upregulated or (B) downregulated in hypoxic conditions. β-actin was used as internal control. Data are presented as the mean ± standard deviation of three independent experiments. *P<0.05 vs. normoxia. CAPN6, calpain 6; CESC, cartilage endplate-derived stem-cells; CFD, complement factor D; COMP, cartilage oligomeric matrix protein; FGF9, fibroblast growth factor 9; IRAK3, interleukin 1 receptor-associated kinase 3; PGAM2, phosphoglycerate mutase 2; TFRC, transferrin receptor; VCAM1, vascular cell adhesion molecule 1.

Enriched biological processes, molecular functions and cellular components in CESCs were identified by GO analysis of the DEGs under normoxia and hypoxia. GO terms were highly enriched in CESCs under hypoxia, such as multicellular organismal process, cargo receptor activity and extracellular region part; the top 10 GO terms from each category are indicated in Fig. 2A-C.

Figure 2. GO term and KEGG pathway enrichment analysis in CESCs under normoxic and hypoxic conditions. The top ten GO functional terms are presented that were identified as enriched in (A) biological process, (B) molecular function and (C) cellular component. (D) Top ten KEGG enriched pathway. CESC, cartilage endplate-derived stem cells; GO, gene ontology; JAK-STAT, Janus kinase-signal transducer and activator of transcription; KEGG, Kyoto Encyclopedia of Genes and Genomes; MAP, mitogen-activated protein kinase.

The KEGG mapping tool was used to identify the enriched functional signaling pathways for the identified DEGs under hypoxia in CESCs. The results revealed that several signaling pathways were significantly enriched, such as drug metabolism-cytochrome P450, hematopoietic cell lineage and rheumatoid arthritis; the top 10 KEGG pathways are in Fig. 2D.

Detection, validation and functional analysis of ASGs under normoxic and hypoxic conditions in CESCs

Genome-wide AS analysis between the normoxic group and hypoxic group identified 9,475 alternatively spliced exons, which belonged to 1,808 ASGs. In addition, 3,677 (39%) of the alternatively spliced exons with an SI value ≥2 were defined as ‘exon inclusion’ events, whereas the remaining 5,798 (61%) exons were considered as ‘exon exclusion’ events (with an SI value ≤-2). Under hypoxic conditions, each ASG had an average of 5.2 (9,475/1,808) alternatively spliced exons, which confirmed the occurrence of multiple alternative splicing events in the same gene. For example, nine alternatively spliced exons were identified for nuclear receptor subfamily 4 group A member 3 (NR4A3), which suggested that AS regulation was complex. In addition, 114 of the 1,808 ASGs were also identified as DEGs, indicating an underlying link between AS and differential gene expression. A total of 10 ASGs were chosen for validation by semi-quantitative PCR analysis, of which 7 were successfully validated. Successful validation was confirmed when the HTA 2.0 array results and semi-quantitative PCR demonstrated that an exon was more likely to be spliced or preserved under hypoxic conditions compared with under normoxic conditions. In the present study, AS events were observed and validated in GDNF family receptor α2, cell adhesion molecule 1, NR4A3, relaxin/insulin-like family peptide receptor 1, NADH:ubiquinione oxidoreductase subunit C1, transmembrane protein 106B and fibronectin 1 under hypoxic conditions compared with under normoxic conditions (Fig. 3).

Figure 3. Validation of the ASGs identified for CESCs under normoxic and hypoxic conditions. CESCs were exposed to normoxic or hypoxic conditions for 21 days. A total of 10 ASGs were selected for further validation by RT-PCR, of which 7 were successfully validated. β-actin was used as an internal control. ASG, alternatively spliced gene; CADM1, cell adhesion molecule 1; CESC, cartilage endplate-derived stem-cells; excl, exclusion; FN1, fibronectin 1; GFRA2, GDNF family receptor α2; H, hypoxia; incl, inclusion; N, normoxia; NDUFC1, NADH:ubiquinione oxidoreductase subunit C1; NR4A3, nuclear receptor subfamily 4 group A member 3; RT-PCR, reverse transcription-polymerase chain reaction; RXFP1, relaxin/insulin-like family peptide receptor 1; TMEM106B, transmembrane protein 106B.

GO term enrichment analysis of ASGs in CESCs under normoxic and hypoxic conditions was conducted to determine the enriched biological processes, molecular functions and cellular components. The results revealed enrichment in GO terms, such as antigen processing and presentation of peptide antigen via major histocompatibility complex (MHC) class I, protein binding and cytoplasm; the top 10 GO terms for ASGs under hypoxic conditions in CESCs are presented in Table III.

Table III. List of the top 10 GO terms identified for alternatively spliced genes under hypoxic conditions in cartilage endplate stem cells.

Table III.

List of the top 10 GO terms identified for alternatively spliced genes under hypoxic conditions in cartilage endplate stem cells.

A, Biological process
GO term	Count	P-value
GO:0002474 antigen processing and presentation of peptide antigen via MHC class I	60	1.97×10−119
GO:0044419 interspecies interaction between organisms	82	2.15×10−116
GO:0006355 regulation of transcription, DNA-dependent	145	7.89×10−110
GO:0006955 immune response	101	1.72×10−103
GO:0019882 antigen processing and presentation	61	5.82×10−100
GO:0006350 transcription	122	4.81×10−79
GO:0007155 cell adhesion	74	2.23×10−70
GO:0007165 signal transduction	128	3.81×10−68
GO:0055114 oxidation reduction	57	7.07×10−54
GO:0006468 protein amino acid phosphorylation	52	1.62×10−44
B, Molecular function
GO term	Count	P-value
GO:0005515 protein binding	496	0
GO:0008270 zinc ion binding	222	2.30×10−211
GO:0046872 metal ion binding	237	1.14×10−174
GO:0000166 nucleotide binding	173	3.72×10−155
GO:0005524 ATP binding	136	9.80×10−133
GO:0005509 calcium ion binding	103	9.50×10−105
GO:0016740 transferase activity	114	1.62×10−91
GO:0003700 transcription factor activity	79	1.40×10−70
GO:0008233 peptidase activity	64	2.37×10−65
GO:0016787 hydrolase activity	93	1.10×10−55
C, Cellular component
GO term	Count	P-value
GO:0005737 cytoplasm	430	0
GO:0005634 nucleus	445	0
GO:0016020 membrane	411	1.05×10−308
GO:0016021 integral to membrane	355	2.29×10−286
GO:0005886 plasma membrane	217	1.30×10−162
GO:0005576 extracellular region	176	1.46×10−157
GO:0042612 MHC class I protein complex	60	1.39×10−115
GO:0005887 integral to plasma membrane	107	5.12×10−99
GO:0005829 cytosol	99	8.65×10−98
GO:0005794 Golgi apparatus	77	5.77×10−72

[i] GO, gene ontology; MHC, major histocompatibility complex.

The 1,808 identified ASGs were further examined by KEGG pathway analysis to identify the enriched functional pathways. These results revealed that several vital pathways were significantly enriched, such as focal adhesion, cell adhesion molecules and axon guidance; a list of the top 10 KEGG pathways identified for ASGs under hypoxic conditions in CESCs is presented in Table IV.

Table IV. List of the top 10 enriched signaling pathways identified by Kyoto Encyclopedia of Genes and Genomes for alternatively spliced genes under hypoxic conditions in cartilage endplate stem cells.

Table IV.

List of the top 10 enriched signaling pathways identified by Kyoto Encyclopedia of Genes and Genomes for alternatively spliced genes under hypoxic conditions in cartilage endplate stem cells.

Pathway	Count	P-value
Focal adhesion	37	1.01×10−22
Cell adhesion molecules	28	2.60×10−19
Axon guidance	26	1.56×10−17
MAPK signaling pathway	35	1.44×10−16
ECM-receptor interaction	21	1.57×10−16
Type I diabetes mellitus	15	2.00×10−14
Small cell lung cancer	19	5.00×10−14
Wnt signaling pathway	23	5.78×10−13
Calcium signaling pathway	25	6.98×10−13
Allograft rejection	13	1.06×10−12

[i] ECM, extracellular matrix; MAPK, mitogen-activated protein kinase.

Discussion

The CEP is an avascular organ, and the oxygen tension within it is as low as 1% (22). A previous study reported that in healthy juveniles and adolescents with no history of low back pain, the CEP was free of blood vessels. Conversely, in adults and seniors with a history of DDD, defects in the CEP were noted through which blood vessels could be observed (23). Blood vessel invasion was also noted in painful IVD and osteoarthritis model (24,25). In degenerated CEP, the destruction of the avascular microenvironment may lead to the destruction of the physiological hypoxic microenvironment of CESCs; therefore, an improved understanding of the effects of the hypoxic microenvironment on CESCs may aid in our understanding of CEP degeneration.

MSC differentiation under hypoxic conditions has previously been investigated: Bone marrow cells at an oxygen concentration of 4–7%, adipocytes at 3.8–9.6% O2 and muscle tissue at a concentration of 1–10% (26–29). In addition, several studies have focused on the regulatory effects of hypoxia on various biological processes in MSCs, including cell cycle, survival, proliferation, differentiation, morphology, immunophenotype and cytogenetics (11,12,30,31). However, the effects of hypoxia on stem cells have not been examined on a genome-wide scale. In the present study, gene expression profiling and AS events in CESCs under hypoxia were examined with the Affymetrix HTA 2.0 system, and the GO terms and KEGG pathways were analyzed with relevant web-based bioinformatics tools.

A total of 448 DEGs were identified between CESCs in the normoxic group and CESCs in the hypoxic group. GO term and KEGG pathway analyses of the DEGs revealed that hypoxia may influence CESC fate in a number of aspects of biological processes, molecular functions, cellular components and signaling pathways, which was consistent with a previous study (20). Several of the identified genes were chosen for validation by PCR, including IRAK3 and CFD. Members of the IRAK protein family have been reported to be necessary components of toll/interleukin-receptor immune signaling pathways (32,33). IRAK3 is mainly expressed in macrophages and monocytes, and negatively regulates the toll-like receptor signaling (34). The CFD protein is an essential component of the alternative complement pathway, and has been reported to serve a role in humoral immunity against infectious agents (35,36). CFD cleaves CFB, following CFB binding to complement component 3b, thereby forming the active complement component 3 convertase (37,38). IRAK3 and CFD are vital components of immune signaling pathways, which indicated their potential necessary roles in the hypoxia-regulated fate of CESCs.

Due to the avascular properties and low immunogenicity of IVD cells, the roles of the immune signaling pathway in IVD cells had not previously been investigated. However, a number of immune signals have been recently demonstrated to exhibit non-immune effects. For example, although the complement pathway is known to serve essential roles in immune surveillance, it has been reported to be involved in apoptotic waste clearance, stem cell differentiation and recruitment of stem cells (39). Several markers of the complement cascades have been demonstrated to be expressed in both MSCs and osteoblasts, such as C5a anaphylatoxin chemotactic receptor (C5aR), C5, C3aR and C3, in addition to cell-surface markers CD59, CD55 and CD46 (40,41). As hypoxia is an important activator of the immune response (42,43), immune signaling may serve an essential role in the CESC response to the hypoxic microenvironment.

The present study aimed to obtain a global view of AS events that occur in CESCs under hypoxic conditions. A total of 9,475 alternatively spliced exons belonging to 1,808 ASGs were detected and analyzed for GO term enrichment. Consistent with the DEG analysis results, both antigen processing/presentation of peptide antigen via MHC class I and immune response were also significantly enriched and ranked within the top 10 GO terms of biological processes. MHC class I protein complex was also significantly enriched and ranked within the top 10 GO terms of molecular function, indicating an abundant occurrence of AS events in the immune signaling pathways. A previous study demonstrated that AS may regulate the immune response by generating different protein isoforms (44); therefore, it is reasonable to consider that AS is highly involved in the immune signaling pathway in CESCs under hypoxia.

The validated AS events from the present study may be considered as examples to highlight the function of AS on immune pathways. For example, glial cell line-derived neurotrophic factor (GDNF) family receptor α2 (GFRA2) is a glycosylphosphatidylinositol-linked surface receptor for both GDNF and neurturin (NTN), and mediates RET receptor tyrosine kinase activation. GFRA2 is a highly-expressed ligand-binding chain in B cells, T cells and monocytes. The RET/GFRA2 signals have been reported to be essential for immune cell development and the protective immunity of central nervous system (45,46). A GFRA2 isoform that lacks exon 2 was revealed to be expressed at low levels in lipopolysaccharide (LPS) + interferon γ (IFN-γ)-activated peripheral blood mononuclear cells (PBMCs), resting monocytes and neuroblastoma cells. Expression of the full-length GFRA2 isoform was observed in resting and (LPS + IFN-γ)-activated PBMCs, resting and (LPS + IFN-γ)-activated monocytes and B cells, particularly following activation with Staphylococcus aureus Cowan I (SAC) (46). The lack of exon 2 resulted in the loss of N-glycosylation sites, which may have influenced the LPS, IFN-γ and SAC-mediated RET/GFRA2 immune signals (46). In the present study, the ratio of variant lacking exon 2/full-length variant was lower in the hypoxia exposed group compared with that in normoxia exposed group, which indicated that the RET/GFRA2 immune signals may be involved in CESCs under hypoxia.

In addition, GO term analysis revealed enrichment in the biological process of DNA-dependent regulation of transcription, the molecular function of nucleotide binding and transcription factor activity. A number of vital regulatory molecules are transcription factor, such as SRY box 9 and runt-related transcription factor 2, the functions of which are to bind the promoter region to facilitate or suppress transcription (47,48). Additionally, the cellular component terms of cytoplasm and nucleus were significantly enriched, which indicated an abundant and frequent occurrence of cytoplasm/nucleus translocation by which many important molecules, such as nuclear factor κB (49) and AKT (50), serve their regulatory roles. These results indicated a relationship between AS events and nuclear signal transduction in CESCs under hypoxia.

KEGG analysis of ASGs identified several AS-related signaling pathways in CESCs under hypoxic conditions, such as focal adhesion, cell adhesion molecules and extracellular matrix (ECM)-receptor interaction. Focal adhesion with the cytoskeleton may influence stem cell morphology, mechanical nature and the direction of differentiation (51). In the process of focal adhesion, signal transduction of matrix-cell adhesion is based on the interaction between stem cells and ECM-receptors (52). The results of KEGG analysis indicated that AS may serve regulatory roles through these mentioned signaling pathways in CESCs under hypoxic conditions.

In conclusion, the present study demonstrated the differentially expressed genes and AS events in CESCs under hypoxia on a genome-wide scale, followed by GO and KEGG analysis. These results provided preliminary insight into the mechanisms of hypoxia-regulated CESC differentiation at the level of gene expression and alternative splicing, which may aid in our understanding of CEP degeneration processes.

Acknowledgements

We sincerely thank Mr Yi Zha for assistance with microarray data analysis. This work was supported by the National Natural Science Foundation of China (grant nos. 81472076, 81271982, 81401801 and 81272028).

Glossary

Abbreviations

Abbreviations:

References