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 |log 2 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.

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

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

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

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.