Identification of key candidate genes and pathways in axial Spondyloarthritis through integrated bioinformatics analysis

Background Radiographic axial Spondyloarthritis (r-axSpA) is the prototypic form of seronegative spondyloarthritis (SpA). In the present study, we evaluated the key genes related with r-axSpA, and then elucidated the possible molecular mechanisms of r-axSpA. Material/Methods The gene expression GSE13782 was downloaded from the GEO database contained five proteoglycan-induced spondylitis mice and three naïve controls. The differentially expressed genes (DEGs) were identified with the Bioconductor affy package in R. Gene Ontology (GO) enrichment and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis were built with the DAVID program followed by construction of a protein-protein interaction (PPI) network performed with Cytoscape. WebGestalt was performed to construct transcriptional regulatory network and microRNAs-target regulatory networks. RT-PCR and immunohistochemical staining were performed to testify the expression of hub genes, transcription factors (TFs) and microRNAs. Results A total of 230 DEGs were identified. PPI networks were constructed by mapping DEGs into STRING, in which 20 hub proteins were identified. KEGG pathway analyses revealed that the chemokine, NOD-like receptor, IL-17, and TNF signalling pathways were altered. GO analyses revealed that DEGs were extensively involved in the regulation of cytokine production, the immune response, the external side of the plasma membrane, and G-protein coupled chemoattractant receptor activity. The results of RT-PCR and immunohistochemical staining demonstrated that the expression of DEGs, TFs and microRNAs in our experiment were basically consistent with the predictions. Conclusions The results of this study offer insight into the pathomechanisms of r-axSpA and provide potential research directions.

Ankylosing spondylitis (axSpA) is a chronic autoimmune disease that involves the axial 3 skeleton (1). AxSpA with spine(2) attachment point inflammation ultimately develops into 4 pathological ossification and rigidity, which is the main cause of disability and the declining 5 quality of life in patients with axSpA. Patients with ankylosing spondylitis have 6 inflammatory pain and stiffness in the joints and spine, severely affecting range of activity. 7 As a complicated chronic autoimmune inflammatory disease, in addition to invading the 8 sacroiliac joint, spine, and spine soft tissue, axSpA can also cause extra-joint manifestations 9 such as anterior uveitis and gut inflammation (3). Its pathological mechanism may be related neural pain signal, which are defined as neuro-immune interactions (5). 16 The influencing factors of axSpA such as gut dysbiosis and the HLA-B27 gene, all 17 primarily participate in its pathological processes (6). The HLA-B27 gene is relevant to the confirmed the high therapeutic effect of TNF-α blockers (golimumab, adalimumab, and 2 etanercept) in the prevention and arrest of onset in patients with axSpA (4), although their 3 prevention ability remains controversial. About 60% of axSpA patients respond to antibiotics 4 (9). Numerous studies have verified the curative effects of non-steroid anti-inflammatory 5 drugs (NSAIDs) on axSpA, although some may increase the risk of 6 developing stomach ulcers (10). NSAIDs are most often used for the treatment of 7 inflammatory pain symptoms because they block Cox enzymes and prostanoids (11). 8 Recently, mesenchymal stem cells were utilised to relieve inflammatory responses and 9 promote tissue regeneration through cell-cell interactions and release of cytokines, which is a 10 novel therapeutic direction (12).

11
In this study, we obtained the candidate DEGs with R language. And then Kyoto 12 Encyclopedia of Genes and Genomes (KEGG) pathway enrichment and PPI network were 13 analysed and visualized. Subsequently, transcription factor (TF)-mRNA regulatory network 14 and miRNA-target genes regulatory network were constructed to identify the key regulatory 15 genes. In conclusion, a total of 230 DEGs, 20 hub genes, 10 miRNAs and 10 TFs were 16 identified.  The raw expression data (Affymetrix CEL files) were first pre-processed into expression 3 values at the gene symbol level. Quantile normalisation was performed with the Robust 4 Multiarray Average normalisation approach of the Bioconductor affy package and the 5 preprocessCore package in R (13, 14) . Finally, the limma package with multiple testing 6 corrections was used to identify significant DEGs among the five PGISp mice and three naïve 7 controls. The classic paired t-test was used to screen the statistically significant DEGs, and p 8 < 0.05 and absolute value of log-fold change |logFC| > 1.5 were set as cut-off thresholds.  19 The STRING search tool was performed to build the PPI network (http://string-db.org/) 20 (16). The DEG list was uploaded to STRING with a interaction score > ≥ 0.4 to build the PPI 21 network. We analysed the numbers of DEGs between PGISp mice and blank control mice 22 and further visualised the PPI network using the Cytoscape software. The nodes, edges of the 23 network and the hub protein, were performed with the cytoHubba plugin (MCODE) with the 6 1 criteria as follows: degree cutoff = 2, node score cutoff = 0.2, max Depth = 100 and k score=2. to build the DEG-associated transcriptional regulatory network. In this TF-targeted gene 9 network, an adjusted P < 0.05 was considered the statistical standard for sorting out the most 10 significant TF-targeted genes.    Table 2). 14 2.9. Statistics 15 All experiments were repeated at least twice. Data are presented as mean±standard deviation. 16 Statistical analyses were performed using Student's t-test. p<0.05 was considered to be 17 significant.  and 372 edges) using the Cytoscape software (Fig. 4A). In the PPI network, the hub genes 16 were sorted by the nodes number. In this study, twenty genes with highest degree scores were 17 proposed as hub genes (Table 1). MS4A6D (degree=28), CTSS (degree=26), FCGR1 18 (degree=25), CD14 (degree=23), and GPR65 had higher degrees in the PGISp mice, and 19 were mainly associated with the NOD-like receptor, IL-17, and TNF signalling pathways.

20
These core genes may play key regulatory roles in the pathogenesis of axSpA disease. In 21 addition, sub-network clustering analyses were performed using the CytoHubba plugin, and 22 two DEG-centred sub-networks were identified (Fig. 4B,C). The the DEGs-related TFs regulatory network was identified with Cystoscope and is 2 presented in Figure 5. A total of 10 TFs including AP1, LSCBP, IRF1, RUNX1, and STAT3 3 and 48 DEGs were involved in the TF regulatory network. Several targeted regulatory genes 4 overlapped between the TF series. The TF-DEG regulatory network consisted of 58 nodes 5 and 77 lines. These TFs exhibited the ability to regulate DEGs which play vital roles in 6 inflammatory response, immune response, and the TNF signalling pathway.

QPCR and histology validation 17
To testify the results of bioinformatics analysis, 9 hub gene, 3 microRNAs and 3 TFs were 18 selected for RT-PCR and IHC respectively. Our findings indicated that the expression levels 19 of the 7 hub genes, 2 microRNAs and 3 TFs were in accordance with the microarray data.  that miR142, miR23, miR452, miR509, and miR365 are closely associated with the immune 6 response, inflammatory response, and bone metabolism. Therefore, further investigations will 7 be beneficial in examining these microRNAs with respect to axSpA.

Conflicts of interest 2
The authors declare that there was no conflict of interest.