Identification of genes and signaling pathways associated with the pathogenesis of juvenile spondyloarthritis.

The aim of the present study was to identify key genes and signaling pathways associated with the pathogenesis of juvenile spondyloarthritis (JSA). The gene expression profile dataset GSE58667, including data from 15 human whole blood samples collected from 11 patients with JSA and four healthy controls, was analyzed to identify differentially expressed genes (DEGs) associated with disease characteristics. Additionally, Gene Ontology term and Kyoto Encyclopedia of Genes and Genomes pathway enrichment analyses of the DEGs were performed. Protein‑protein, microRNA‑transcription factor and chemical‑gene interaction networks were constructed. A total of 326 DEGs, 196 upregulated and 130 downregulated, were identified. DEGs, including C‑X‑C motif chemokine ligand 5 (CXCL5), BCL2 interacting protein 3 like (BNIP3L), dual specificity phosphatase 5 (DUSP5) and tumor protein p53 (TP53) were enriched in functions associated with apoptosis, the cell cycle and immune responses. KEGG pathway enrichment analysis revealed that pathways associated with inflammation and the mitogen‑activated protein kinase 1 (MAPK) signaling pathway were the most enriched by DEGs. The results of the present study indicated that the MAPK signaling pathway and four genes, including CXCL5, BNIP3L, DUSP5 and TP53, may be implicated in the pathogenesis of JSA.

Introduction

Spondyloarthritis (SpA) includes a group of immune-mediated inflammatory diseases with similar genetic and clinical manifestations, including ankylosing spondylitis (AS), psoriatic arthritis (PsA) and juvenile SpA (JSA) (1,2). In particular, JSA is a group of chronic inflammatory diseases associated with human leukocyte antigen B27, affecting children at ≤16 years of age. The prevalence of spondyloarthropathy is ~1.9% worldwide (3,4). The primary clinical manifestations of JSA include enthesitis and peripheral arthritis, although symptoms including acute anterior uveitis, bowel inflammation, psoriasis and cardiac disease may affect certain patients with JSA (5). An increased incidence is observed among males compared with females (3). Patients with JSA are likely to develop AS during their disease course (6). Therefore, research into the mechanisms underlying JSA is necessary for development of improved diagnostic and treatment approaches.

A previous study suggested that certain genes, including Toll-like receptor 4, NLR family pyrin domain containing 3, C-X-C motif chemokine receptor 4 (CXCR4), dual specificity phosphatase 6, mitogen-activated protein kinase kinase 2 and protein tyrosine phosphatase, non-receptor type 12 are involved in the development of JSA (2). CXCR4 was additionally hypothesized to be involved in the pathogenesis of SpA (7). In one study, CXCR4 and its ligand CXCL12 were demonstrated to serve a role in retaining inflammatory cells within the joint (8). Previous studies demonstrated that CXCR4 and CXCL12 are involved in the development of inflammatory bowel disease (IBD) (8,9). Furthermore, overexpression of CXCR4 in patients with JSA may lead to mucosal disregulation and IBD (2). The HOX signaling pathway and the granule cell survival pathway are additionally associated with the development of JSA (2). However, there has been a limited number of genetic studies in patients with JSA and these studies were limited by small numbers of patients (6). At present, there are no effective drugs or therapies available for the treatment of patients with JSA. As a consequence, it is necessary to study the mechanism underlying the development of JSA and identify targets for more effective therapies using bioinformatics analysis.

The present study aimed to screen target genes and pathways implicated in development of JSA. Differentially expressed genes (DEGs) in patients with JSA were identified by screening the dataset GSE58667 (2). Furthermore, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway and Gene Ontology (GO) term analyses were performed to verify the pathways associated with JSA. Subsequently, protein-protein interactions (PPI), transcription factors (TF), microRNA (miRNA) and small molecule interaction network analyses were performed to screen for genes and molecules associated with JSA. The results of the present study may provide information for further investigation of the mechanisms underlying JSA and for the development of potential treatment methods.

Materials and methods

Microarray data

mRNA expression profile dataset GSE58667 was obtained from the National Center of Biotechnology Information (NCBI) Gene Expression Omnibus (GEO) database (www.ncbi.nlm.nih.gov/geo). The data were based on the GPL570 (HG-U133_Plus_2) Affymetrix Human Genome U133 Plus 2.0 Array platform (Affymetrix; Thermo Fisher Scientific, Inc., Waltham, MA, USA). A total of 15 human whole blood samples were collected from 11 patients with JSA and four healthy controls. In the present study, data from the 11 whole blood samples from patients with JSA were used as the experimental group. The expression data obtained from four whole blood samples from healthy controls were used as the control group.

Data preprocessing and analysis of DEGs

Raw data were downloaded in CEL format and the robust multiarray average method (10) in the Affy package (version 1.540) (11) was used to pre-process the mRNA expression data by performing background adjustment, quantile normalization and summarization. The differentially expressed probe in the experimental group compared with the control group was screened using the limma (ver.3.32.3) package (12). The cutoff thresholds for identification of differentially expressed probes were P<0.05 and |log2 fold change (FC) |>0.585, and the probes were annotated. Ultimately, a heatmap was generated using the pheatmap package in R (version 3.25; http://www.R-project.org/).

GO and pathway enrichment analyses

In the present study, GO and KEGG pathway analyses for DEGs were performed using the Database for Annotation, Visualization and Integrated Discovery (version 6.8; david.ncifcrf.gov) (13). The number of enriched genes >2 and P<0.05 were selected as cut-off criteria.

PPI network construction and analysis

The Search Tool for the Retrieval of Interacting Genes (STRING; version 10.0; string-db.org) database was used for the construction of the PPI network between DEGs and other genes (14). Required confidence (combined score) >0.7 was selected as the threshold value for PPIs. The PPI networks were constructed using Cytoscape software (version 3.3.0; www.cytoscape.org) (15) and the network was analyzed to identify hubs according to the ranking of network connectivity.

In the PPI network, proteins with similar functions tend to cluster together, and the node function was correlated with the distance between one node and another node in the network. Therefore, the bioinformatics analysis in the present study allowed for the identification of unknown functions of proteins by studying protein complexes or clusters of functional subnetworks in the PPI network. The subnetworks were evaluated with the Molecular Complex Detection (MCODE; version 1.4.2; apps.cytoscape.org/apps/mcode) (16) plugin in Cytoscape. Subsequently, the functions of subnetworks were evaluated by GO and pathway enrichment analyses.

Construction and analysis of a DEG-miRNA-TF regulatory network

TFs and miRNAs associated with DEGs were predicted based on the Encyclopedia of DNA Elements (http://amp.pharm.mssm.edu/Enrichr/#stats) and ChIP-X Enrichment Analysis consensus TFs (17) and TargetScan microRNA library enrichment analysis in the Enrichr database (amp.pharm.mssm.edu/Enrichr) (18). The screening threshold was set as adjusted P<0.01. Within the predicted TFs, the differentially expressed TFs in the present study were chosen for the downstream analysis. Finally, the target gene network of TFs and miRNAs was established.

Small molecule target network analysis

The Comparative Toxicogenomics Database (CTD, ctdbase.org) (19) includes chemical-gene interactions associated with spondylitis. The present study identified the chemical small molecules associated with the DEGs and constructed the chemical-gene interaction network.

Results

Analysis of DEGs

To determine DEGs between patients with JSA and controls, the publicly available microarray dataset GSE58667 was accessed from the GEO database and differences in expression were determined using the limma package. A total of 326 DEGs were identified using the threshold of P<0.05 and |logFC| >0.585, including 196 upregulated DEGs and 130 downregulated DEGs (Fig. 1).

Figure 1.

Volcano plot of differentially expressed genes. The blue dots represent differentially expressed genes (TRUE) and the red dots represent non-differentially expressed genes (FALSE).

Pathway enrichment analysis

The results demonstrated that a total of 128 GO terms or KEGG pathways were enriched. Based on the most significant P-values, the top five GO terms and KEGG pathways were selected. The DEGs were significantly enriched in five KEGG pathways, including pathways associated with IBD, the cell cycle and the mitogen-activated protein kinase (MAPK) signaling pathway (Fig. 2). GO analysis indicated that the DEGs were significantly enriched in apoptosis-, cell cycle- and immune response-associated terms. In particular, CXCL5 was revealed to be enriched in functions associated with immune response as well as pathways associated with rheumatoid arthritis.

Figure 2.

GO and KEGG pathway enrichment analyses of DEGs. Count represents the number of DEGs enriched in each term. The black trend line represents the value of-log10 (P-value). GO, Gene Ontology; BP, biological process; CC, cellular component; MF, molecular function; DEGs, differentially expressed gene; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Interaction network analysis of DEGs

STRING was used to analyze PPIs between DEGs. Subsequently, the PPI network was constructed, which included a total of 164 nodes and 762 interaction pairs (Fig. 3). Furthermore, nodes with highest degrees in the network were analyzed. Table I presents the top 30 target genes with the highest degrees, all downregulated, including BUB1 mitotic checkpoint serine/threonine kinase (BUB1), tumor protein p53 (TP53), aurora kinase B (AURKB) and cell division cycle 20 (CDC20).

Figure 3.

Protein-protein interaction network analysis of DEGs. The pink nodes represent upregulated DEGs. The green nodes represent downregulated DEGs. DEGs, differentially expressed genes.

Table I.

Top 30 nodes with the highest degrees.

Gene	Degree	Differential expression
BUB1	39.0	Downregulated
AURKB	38.0	Downregulated
CDC20	38.0	Downregulated
UBE2C	38.0	Downregulated
CCNB1	38.0	Downregulated
TPX2	38.0	Downregulated
TOP2A	38.0	Downregulated
NCAPG	37.0	Downregulated
RRM2	36.0	Downregulated
BIRC5	35.0	Downregulated
MELK	35.0	Downregulated
CCNB2	34.0	Downregulated
CHEK1	33.0	Downregulated
HMMR	33.0	Downregulated
DTL	33.0	Downregulated
CDKN3	33.0	Downregulated
DLGAP5	33.0	Downregulated
BUB1B	32.0	Downregulated
TYMS	32.0	Downregulated
NUSAP1	32.0	Downregulated
ASPM	32.0	Downregulated
GINS2	31.0	Downregulated
ZWINT	31.0	Downregulated
CEP55	31.0	Downregulated
MCM4	30.0	Downregulated
FOXM1	30.0	Downregulated
TP53	30.0	Downregulated

BUB1, BUB1 mitotic checkpoint serine/threonine kinase; AURKB, aurora kinase B; CDC20, cell division cycle 20; UBE2C, ubiquitin conjugating enzyme E2 C; CCNB1, cyclin B1; TPX2, TPX2, microtubule nucleation factor; TOP2A, DNA topoisomerase IIα; NCAPG, non-SMC condensin I complex subunit G; RRM2, ribonucleotide reductase regulatory subunit M2; BIRC5, baculoviral IAP repeat containing 5; MELK, maternal embryonic leucine zipper kinase; CCNB2, cyclin B2; CHEK1, checkpoint kinase 1; HMMR, hyaluronan mediated motility receptor; DTL, denticleless E3 ubiquitin protein ligase homolog; CDKN3, cyclin dependent kinase inhibitor 3; DLGAP5, DLG associated protein 5; BUB1B, BUB1 mitotic checkpoint serine/threonine kinase B; TYMS, thymidylate synthetase; NUSAP1, nucleolar and spindle associated protein 1; ASPM, abnormal spindle microtubule assembly; GINS2, GINS complex subunit 2; ZWINT, ZW10 interacting kinetochore protein; CEP55, centrosomal protein 55; MCM4, minichromosome maintenance complex component 4; FOXM1, forkhead box M1; TP53, tumor protein p53.

Furthermore, the highest scoring module identified in the PPI network, which scored equal to 29.806, was acquired through the analysis of PPI network via the MCODE plug-in. A total of 32 nodes and 462 interaction pairs were included in this module (Fig. 4). Enrichment analysis was conducted for the nodes in the module. A total of 69 GO functions and KEGG pathways were enriched (data not shown).

Figure 4.

A module including a sub-network in the protein-protein interaction network. The green nodes represent downregulated differentially expressed genes.

Construction and analysis of the DEG-miRNA-TF regulatory network

The DEG-miRNA-TF regulatory network was constructed based on TFs, miRNAs and DEGs from the Enrichr database (Fig. 5). There were three TFs (Table II), seven miRNAs (Table III) and 132 target genes in the network. The greatest number of genes was regulated by miR-23a (overlap=18), including BCL2 interacting protein 3 like (BNIP3L), C-X-C motif chemokine ligand 5 (CXCL5) and dual specificity phosphatase 5 (DUSP5). The above results indicated that miR-23a may serve a role in the development of JSA. Furthermore, 22 target genes were regulated by TFs and miRNAs simultaneously, including BNIP3L, DUSP5, E1A binding protein p300 (EP300) and F-box protein 11 (FBXO11).

Figure 5.

DEG-miRNA-TF regulatory network of DEGs. The green V-shaped symbol represents the downregulated TFs, and the pink V-shaped symbol represents the upregulated TFs. Blue color represents miRNAs. Circular symbols represent target genes, including pink nodes (upregulated DEGs) and the green nodes (downregulated DEGs). The arrows represent the regulatory associations between TFs and their targets. TF, transcription factors; miRNA, micro RNA; DEG, differentially expressed gene.

Table II.

TFs in the microRNA-TF regulatory network of differentially expressed genes.

TF	Gene count	P-value	Adjusted P-value
UBTF_ENCODE	64	3.16×10−11	1.61×10−9
FOXM1_ENCODE	13	4.75×10−9	1.21×10−7
CREB1_CHEA	40	7.04×10−4	2.65×10−3

Gene Count represents the number of genes regulated by the corresponding transcription factor. ENCODE, Encyclopedia of DNA Elements; TF, transcription factor; UBTF, nucleolar transcription factor 1; FOXM1, forkhead box protein M1; CREB1, cyclic AMP-responsive element-binding protein 1; CHEA, ChIP-X Enrichment Analysis.

Table III.

miRNAs in the miRNA-transcription factor regulatory network of differentially expressed genes.

miRNA	Overlap	P-value	Adjusted P-value
miR-369-3P	16	3.32×10−7	3.23×10−5
miR-410	11	3.38×10−7	3.23×10−5
miR-374	15	7.38×10−5	4.24×10−3
miR-320	14	8.88×10−5	4.24×10−3
miR-183	11	1.45×10−4	5.56×10−3
miR-23A	18	2.17×10−4	6.93×10−3
miR-23B	18	2.17×10−4	6.93×10−3

miR/miRNA, microRNA.

Target network analysis of chemical small molecules

Following construction of the chemical-target gene interaction network, several small molecules potentially associated with JSA were identified. A total of nine small molecules may be associated with JSA, including methotrexate, flurbiprofen, sulfasalazine and naproxen (Fig. 6).

Figure 6.

Chemical-target interaction network. The blue squares represent the chemicals, the pink nodes represent the upregulated genes and the green nodes represent downregulated genes.

Discussion

Genetic and environmental factors are involved in the occurrence of JSA, although the exact etiology and pathogenesis of JSA remain to be elucidated (5). In the present study, bioinformatics analyses were performed to determine target genes and pathways implicated in the development of JSA. The results of the present study suggested that miR-23a was likely to be involved in the development of JSA. Additionally, arthritis and MAPK signaling pathway-associated terms were demonstrated to serve roles in JSA. Furthermore, the genes BNIP3L, CXCL5, DUSP5 and TP53 were demonstrated to be involved in the pathogenesis of JSA. In addition, a total of nine small molecules were screened, which may either promote or inhibit the development of JSA.

In the present study, the regulatory DEG-miRNA-TF network was constructed, and it included three TFs, seven miRNAs and 132 target genes. miR-23a regulated the greatest number of genes (overlap=18). A previous study provided evidence that the miR-23a cluster was significantly decreased in the PsA synovium (20). PsA and JSA are members of the SpA family of inflammatory diseases (21). Therefore, dysregulation of miR-23a expression may be involved in the development of JSA.

Three target genes were regulated by miR-23a, including two upregulated genes, BNIP3L and CXCL5, and a downregulated gene, DUSP5. The target gene BNIP3L, was regulated by a TF (CREB1) and miR-23a simultaneously. BNIP3L, a member of the BCL2/adenovirus E1B 19 kd-interacting protein family, is associated with apoptosis, autophagy and mitophagy (22). A previous study demonstrated that BNIP3L induces apoptosis via regulation of the mitochondrial membrane permeability (23). Furthermore, BNIP3L has been demonstrated to serve a role in cardiomyocyte apoptosis (24). Lamot et al (2) revealed that blood cells isolated from patients with JSA demonstrated a higher activity of numerous processes associated with apoptosis, which indicated broad cellular dysfunction in blood cells. Thus, it was speculated that BNIP3L may be associated with JSA progression via involvement in the apoptosis process. However, the exact underlying mechanisms remain unclear and require further investigation.

CXCL5, a member of the CXC subfamily of chemokines, is involved in angiogenesis, tumor growth and metastasis (25). Expression of CXCL5 has been identified in animal models and patients with gastrointestinal inflammation (26,27). Mei et al (28) demonstrated that CXCL5 served a role in innate immunity by regulating neutrophil homeostasis. The present study demonstrated that CXCL5 was enriched in the functions of immune response and pathways associated with rheumatoid arthritis (data not shown). Furthermore, in the present study CXCL5 was associated with the CXCR4 gene in the PPI network, and CXCR4 has been previously demonstrated to be overexpressed in patients with JSA (2). Therefore, the results of the present study suggested that inflammatory the mediator CXCL5 may be implicated in the pathogenesis of JSA, by mediating inflammation and the innate immune response.

DUSP5 is a nuclear phosphatase protein and its expression is induced by cytokines, stress and other stimulatory factors (29). In the present study, DUSP5 was enriched in the MAPK signaling pathway, which regulates cell proliferation, apoptosis and the inflammatory response (30,31). A previous report demonstrated that activation of the MAPK pathway may be associated with the inhibition of DUSP5 (32). Phosphoserine and phosphotyrosine residues of MAPK may be specifically dephosphorylated by DUSP5 (32). A previous study indicated that the MAPK signaling pathway is likely to serve a role in the pathogenesis of SpA (2). Therefore, the authors of the present study hypothesize that the MAPK signaling pathway may serve a role in the pathogenesis of JSA, mediated by DUSP5.

It was previously demonstrated that TP53 is downregulated in several autoimmune diseases, including rheumatoid arthritis, systemic lupus erythematosus, multiple sclerosis and type 1 diabetes (33). A previous study demonstrated that activation of TP53 led to numerous responses in cells, including cell cycle arrest (34). TP53 was enriched in four KEGG pathways in the present study, including ‘cell cycle’, ‘MAPK signaling pathway’ and ‘p53 signaling pathway’. Lamot et al (2) demonstrated that patients with JSA exhibited decreased activity of processes associated with the cell cycle. In particular, previous data demonstrated that TP53 is able to functionally interact with the MAPK signaling pathway (34). In the present study, TP53 was associated with the BUB1 gene, which exhibited the highest degree in the PPI network. These results suggested that TP53 may be involved in the pathogenesis of JSA by regulating the cell cycle and MAPK signaling.

The identification of small molecules with potential therapeutic efficacy for the treatment of JSA was the aim of the present study. A total of nine small molecules were associated with JSA, including methotrexate, sulfasalazine, indomethacin and ibuprofen. Previously, a randomized placebo-controlled study indicated that sulfasalazine was effective in patients with JSA (35). Meanwhile, as arthritis in JSA is predominantly peripheral, methotrexate has also been used in the treatment for patients with juvenile idiopathic arthritis (36,37). Multiple small molecules were identified in the present study; however, their roles in the treatment of JSA require further investigation.

In conclusion, the results of the present study indicated that miR-23a may be implicated in the development of JSA. A total of three target genes (BNIP3L, CXCL5 and DUSP5) regulated by miR-23a and associated with MAPK signaling were identified in the present study, and may serve roles in the pathogenesis of JSA. Furthermore, TP53 was identified, which was implicated in four KEGG pathways and may be involved in the pathogenesis of JSA. However, in vitro and in vivo studies are required to confirm the role of the identified genes and pathways in the pathogenesis of JSA.