Abstract: Background and Objective: Interferon (IFN)-γ plays crucial roles in the development of ankylosing spondylitis (AS) though regulating its sensitive response genes, but the sensitive response genes to IFN-γ dysregulation in AS are partially reported. Thus, the present study explored potential target genes for AS treatment under interferon-γ (IFN-γ) stimulated condition. Materials and Methods: The gene expression profile, GSE11886 was downloaded. The differentially expressed genes (DEGs) in non-INF-γ group and INF-γ group were revealed respectively, followed by function and pathway enrichment analyses, as well as protein-protein interaction network analysis. Furthermore, the small-molecule drugs associated with AS were identified. Finally, the data validation for DEGs in AS with non-INF-γ interference was conducted using another dataset GSE25101. Results: Totally, 317 DEGs and 394 DEGs were revealed in the INF-γ group and non-INF-γ group, respectively. Venn diagram revealed 60 overlapped DEGs. Among the 394 DEGs, nineteen genes were verified in non-INF-γ group. The DEGs in the two groups were mainly enriched in mitotic cell cycle and T cell activation, respectively. CKD1 and BUB1 were outstanding in INF-γ group, while the validated PRKCQ and SMARCA4 might be crucial with AS in non-INF-γ group. In addition, the DL-thiorphan and NS-398 had the highest negative correlation with DEGs in two groups, respectively. Conclusion: Genes like CDK1 and BUB1 might be novel INF-γ sensitive response gene for AS. The PRKCQ and SMARCA might be strongly related to the process of AS. In addition, NS-398 and DL-thiorphan might be ideal small-molecule drugs for AS treatment.
INTRODUCTION
Ankylosing spondylitis (AS) is a chronic inflammatory autoimmune disease of the axial skeleton1. Inflammatory back pain, stiffness and limited mobility are the common outcomes of AS in clinical2. Genetic, chronic infections, autoimmune disorders and endocrine disorders are all proved to be related to the process of AS3. The prevalence is between 0.1 and 0.2% of the general population and an overall increased risk is 50% for vascular mortality4. Although treatments and medications can reduce symptoms and pain to some extent, the AS is incurable till now5,6.
AS is one of the seronegative spondyloarthropathies which has a strong genetic predisposition7. Previous studies indicated that there is a strong relationship between gene expression and AS process8,9. Wang et al.10 indicated that genes down-regulated by Tim-3 are involved in the susceptibility to AS. The association of the ANTXR2 gene polymorphism and AS in China has been revealed11. Actually, various biomarkers have been explored for AS treatment based on gene expressions12. Gene expression profile has been widely used for the pathobiology of AS13,14. Gene expression analysis of macrophages derived from AS patients reveals the dysregulation of interferon-γ (INF-γ)15. IFN-γ is an important activator of macrophages and inducer of Class II major histocompatibility complex molecule expression16. Aberrant IFN-γ expression is associated with a number of autoimmune diseases like AS17. Although IFN-γ binding polypeptides is a useful therapy of IFN-γ related diseases18, the effect of INF-γ on the peripheral blood macrophage cells of AS patients is still unclear.
Based on the gene expression profile (GSE11886) of peripheral blood macrophages tissue samples from AS patients and healthy control subjects, a gene expression analysis has been performed by Smith et al.19. In that study, the author believes that there is a relative defect in IFN-γ gene regulation, with autocrine consequences and implications for disease pathogenesis. However, based on the huge gene expression profile, little information is focused on the genes that are sensitive to IFN-γ interference in AS patients. Furthermore, the effect of IFN-γ on the expression of macrophage in AS patients is still unclear. Thus, in the present study, the significant abnormal expressed genes, pathways as well as protein-protein interaction (PPI) network were investigated, followed by small-molecule drug analysis. With these comprehensive analyses, it was tried to explore several potential genes for AS treatment under normal or IFN-γ stimulated condition.
MATERIALS AND METHODS
Data resource: Gene expression profile dataset GSE1188619 was downloaded from the Gene Expression Omnibus (GEO) database20 and the sequencing data of GSE11886 was produced by using the platform GPL6480 [HG-U133–Plus–2] Affymetrix Human Genome U133 Plus 2.0 Array. The dataset included 33 samples from human peripheral blood macrophages tissue. All these samples were divided into two groups: Non-INF-γ group (8 AS patients and 9 healthy control subjects without INF-γ interference) and INF-γ group (7 AS patients and 9 healthy control subjects with INF-γ interference). In addition, the present study was enrolled at department of Orthopedics, Honghui Hospital and was carried out during December, 2016 to May, 2017.
Data preprocessing and differential expression analysis: After raw data in the format of CEL were obtained, the k-nearest neighbor (KNN) algorithm21 was used for the replacement of missing values. Then the profile data were normalized using the Median clustering method22. The differentially expressed genes (DEGs) between healthy macrophages and AS macrophages in both two groups were identified using the Linear Models for Micro array Data (limma) package in R (v.3.0.0)23. The p-value<0.01 was selected as the threshold for the DEG selection.
Principal component analysis for DEGs: Principal component analysis (PCA) is a one of techniques for clustering gene expression data using multi variate technique that analyzes a data table in which observations are described by several inter-correlated quantitative dependent variables24. This process is defined in such a way: (i) The first principal component has the largest possible variance (that is, accounts for as much of the variability in the data as possible) and (ii) Each succeeding component in turn has the highest variance possible under the constraint25. Principal axes of a set of observed data vectors may be determined via maximum-likelihood estimation of parameters in a latent variable model closely related to factor analysis. To reveal whether there was an ideal distinction of DEGs obtained in either INF-γ group or non-INF-γ group, the PCA analysis was performed on the DEGs obtained above using the prcomp functions in R software26,27.
Venn diagram analysis for DEGs identified in INF-γ group and non-INF-γ group: To reveal the difference of DEGs obtained in INF-γ group and non-INF-γ group, the overlapped or the specific DEGs in the two groups were identified using Venn diagram28.
Function and pathway enrichment analysis of DEGs: The Database for Annotation, Visualization and Integrated Discovery (DAVID)29 is a gene functional classification tool that provides multiple functional annotation tools for investigators to understand biological meaning behind a large list of genes. Gene function includes three functional categories: Molecular function (MF), biological process (BP) and cellular component (CC)30. In the study, BP function enrichment analysis for the DEGs was performed using the DAVID online tool. The BP terms with criterion of p-value<0.05 were considered significant.
Besides, KOBAS2.0 is a web server for annotating an input set of genes with putative pathways and disease31. Based on the hyper geometric distribution, the KOBAS 2.0 was used to perform pathway enrichment analysis. P-value<0.05 was considered as the cut-off value for selection significant enriched pathways.
Protein-protein interaction network construction: The Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) is a biological database and web resource of known and predicted protein interaction data32. According to information in the STRING database, the protein-protein interaction (PPI) networks of DEGs in each group were constructed, with the criterion of combined score >0.8 and then was visualized by the Cytoscape software33. Meanwhile, the nodes possessing a mass of connections were defined as hub nodes.
Small-molecule drug screening: Connectivity map (Cmap) database34 is a novel tool for using gene expression signatures to connect small-molecule drugs, its target genes and disease. In order to screen small-molecule drugs that are closely associated with gene expression and the manifestations of AS, the above DEGs divided into up-regulated and down-regulated groups were further mapped to the Connectivity map, respectively. The absolute value of relevance scores >0.8 was considered as high connection between small-molecule drug and AS.
Validation of the expression level of DEGs in AS samples: Another dataset GSE25101 which was also downloaded from GEO, was utilized for validation the sequencing data of samples in non-INF-γ group and based on the on the platforms of Illumina Human HT-12 V3.0 expression bead chip (GPL6947). This dataset contained 16 whole blood samples from AS patients and 16 normal whole blood samples and all were non-INF-γ treated. After obtaining the data of GSE25101, the data preprocessing and differential expression analysis were used same methods and threshold value used above. Subsequently, the overlapped DEGs in the GSE11886 and GSE25101 datasets were obtained using Venn diagram analysis and those overlapped DEGs were considered as the genes preliminarily validated of the present study.
RESULTS
Selected DEGs in INF-γ and non-INF-γ group: A total of 317 DEGs were obtained from samples between 7 AS patients and 9 healthy control subjects with INF-γ interference, while 394 DEGs were obtained from samples between 8 AS patients and 9 healthy control subjects without INF-γ interference.
Investigation of DEGs based on PCA: The results showed that the AS samples and healthy samples could be well distinguished, which indicated that the identified DEGs in each group were representative for identifying the difference between the samples and there was a good distinction of DEGs to divide samples according to different sample type in either INF-γ group or non-INF-γ group (Fig. 1).
Overlapped DEGs analysis by Venn diagram: Based on the Venn diagram analysis, totally257 and 334 proprietary DEGs were revealed in INF-γ group or non-INF-γ group, respectively. Furthermore, a set of 60 overlapped DEGs were identified in the two groups (Fig. 2).
Pathways and function enrichment investigation of DEGs in the network: To gain insights into DEGs-related pathways, the pathway enrichment analysis was carried out by using the KOBAS software (Table 1). For non-INF-γ group, the DEGs were mainly enriched in three pathways including June 30, 2018Cytokine-cytokine receptor interaction (hsa04060, p = 7.20E-03), Intestinal immune network for IgA production (hsa04672, p = 1.14E-02) and T cell receptor signaling pathway (hsa04660, p = 1.30E-02).For INF-γ group, the DEGs were mainly enriched in Oocyte meiosis pathway (hsa04114, p = 1.05E-02). Furthermore, the overlapped DEGs in two groups were enriched in Hematopoietic cell lineage pathway (hsa04640, p = 4.00E-03). To further explore gene function of the DEGs in PPI network, BP that enriched by DEGs were investigated (Table 2).
| Fig. 1(a-b): | Principal component analysis (PCA) of differentially expressed gene: (a) The result of PCA analysis on INF-γ interference group, (b) The result of PCA analysis on non-INF-γ interference group, Red and orange represent the AS samples, while white and gray represent the normal samples |
With the criterion of p-value<0.05, the DEGs in INF-γ group were mainly enriched in 7 outstanding Bps like mitotic cell cycle (0000278, p = 3.23E-04), fatty acid metabolic process (0006631, p = 4.41E-04) and cell cycle phase (0022403, p = 4.82E-04). On the other hand, the DEGs in INF-γ group were mainly enriched in 8 outstanding BPs such as positive regulation of immune system process (0002684, p = 4.12E-08), positive regulation of leukocyte activation (0002696, p = 2.51E-07) and positive regulation of T cell activation (0050870, p = 1.50E-05).
PPI network analysis of DEGs: According to aforementioned criterion, two PPI networks were constructed. The network of INF-γ group was consisted of 29 connections and 24 DEGs (including 10 up-regulated genes and 14 down-regulated genes) (Fig. 3a), while that of non-INF-γ group was comprised of 69 connections and 59 DEGs (including 25 up-regulated genes and 34 down-regulated genes) (Fig. 3b). Through calculating the degree of each node in each group, BUB1 (degree = 5), ASPM (degree = 5) and CDK1 (degree = 4) were the hub nodes with high degree in the network of INF-γ group, whereas the IFNG (degree = 10), stat1 (degree = 9) and LCK (degree = 9) were the hub nodes with high degree in the network of non-INF-γ group.
| Fig. 2: | Venn diagram for differentially expressed genes (DEGs) in INF-γ interference group and non-INF-γ interference group: Green represents the INF-γ interference group, while purple represents the non-INF-γ interference group |
| Table 1: | Results of pathway enrichment analysis based on differentially expressed gene |
INF(-) represented the non-INF-γ interference group, INF(+) represented the INF-γ interference group, p<0.05 was the cut-off value | |
| Table 2: | Results of functional enrichment analysis based on differentially expressed gene |
INF(-) represented the non-INF-γ interference group, INF(+) represented the INF-γ interference group, Count, the number of gene enriched in target function, p<0.05 was the cut-off value | |
| Table 3: | Result of small-molecule drug investigation |
Predicting small-molecule drugs related to the DEGs: With the criterion of absolute value of relevance score >0.8, a total of 7 and 6 small-molecule drugs were revealed in INF-γ group and non-INF-γ group, respectively (Table 3). Among them, the DL-thiorphan was negatively related to expression of the DEGs in non-INF-γ group (p = 2.17E-03), while NS-398 was negatively associated with those in INF-γ group (p = 2.60E-04). Considering that these gene alterations might be the causative factors for AS development, while the above two molecular drugs could restore these genes expressions, it was speculated that they might have a protective role in AS progression.
Validation of the expression level of DEGs in AS samples: To confirm the DEGs identified in AS samples from non INF-γ treated group, the DEGs were validated in another micro array dataset. Totally 19 overlapped DEGs were found. In addition, among the 59 DEGs in the PPI network (Fig. 3b), two crucial AS related genes (SMARCA4 and PRKCQ) had the same expression trend in the validated dataset.
DISCUSSION
IFN-γ is proved to be associated with the process of AS. In this study, potential genes associated with the development of AS under normal or IFN-γ stimulated condition were investigated through the bioinformatics analysis. A total of 317 DEGs and 394 DEGs were identified in the INF-γ group and non-INF-γ group between the AS patients and healthy controls, respectively. DEGs in non-INF-γ group and INF-γ group were mainly enriched in mitotic cell cycle and T cell activation, respectively. CDK1 and BUB1 were outstanding in INF-γ group. In addition, validated the DEGs identified in the non-INF-γ group were validated by using another dataset GSE25101 and 59 DEGs (e.g., SMARCA4 and PRKCQ) in the PPI network had the same expression trend in the validated dataset. Moreover, the DL-thiorphan and NS-398 had the highest negative correlation with DEGs in non-INF-γ group and INF-γ group, respectively.
CDK1 encoded protein is a highly conserved protein that regulates the cell cycle or general RNA polymerase II functions during gene transcription35. Recent data suggest that CDK1 as an essential Cdk in cell cycle can regulate the proliferative senescence of human diploid fibroblasts36. BUB1 is a serine/threonine protein kinase that plays a key role in the establishment of the mitotic spindle checkpoint and chromosome congression37. It is proved to be associated with abnormal proliferation of cells in various kinds of tumors38.
| Fig. 3(a-b): | Protein-protein interaction (PPI) network of the differentially expressed genes: (a) The PPI network of the INF-γ interference group, (b) The PPI network of the non-INF-γ interference group, Red triangle represents the up-regulated gene, Green triangle represents the down-regulated gene. The node represents a protein encoded by differentially expressed gene and lines represent interactions of the proteins |
In addition, the inhibition of BUB1 may contribute to genomic instability and abnormal proliferation of fibroblast. Actually, the cell proliferation, cycle and apoptosis of fibroblast is remarkably involved with AS39. Besides, it has been suggested that AS is concerned with fibroblasts ossification and cell proliferation and apoptosis of fibroblast might be a crucial role for accelerating cells ossification. In the present study, the CDK1 and BUB1 enriched in the function annotation of mitotic cell cycle were hub genes in the PPI network of INF-γ group. Thus, it was supposed that CDK1 and BUB1 might be closely correlated with the progression of AS via regulating the proliferation of fibroblast. In addition, the level INF-γ in AS patients is significantly lower than normal subject40. Stimulation of INF-γ in macrophage cells for peripheral blood in AS patients might lead to recover the normal expression of CDK1 and BUB1, which further indicated that CDK1 and BUB1 might be novel INF-γ sensitive response gene for AS. However, a further investigation is needed to confirm this speculation.
PRKCQ encodes protein kinase c theta, which is essential for T cell activation via regulating IL-2 production41,42 and it is indicated that IL-2 is significant overproduction in AS patients than healthy person43. Similarly, current result showed that PRKCQ was enriched in the T cell receptor signaling pathway and go term of positive regulation of T cell activation. SMARCA4, also known as BRG1, may function as a pro-inflammatory genes via regulation of SWI/SNF assembly and function44. It is consistent with another study that BRG1 deletion can impair T cell activation which may cause an serious inflammation45. Due to inflammation is a major cause for AS progression, it was supposed that PRKCQ and BRG1 may be closely associated with AS under the non-INF-γ stimulation via regulating T cell activation and inflammation.
Various clinical drugs such as celecoxib, a cyclooxygenase 2-specific (COX-2) inhibitor, are believed to play the important role in inhibiting the process of AS46. Compared with the "large molecule" biologics, one advantage of small-molecule drugs is that many of them can be taken orally, whereas biologics generally require injection or another parenteral administration47. In the present study, the NS-398 and DL-thiorphan had the highest negative correlation with DEGs in INF-γ group and non-INF-γ group, respectively. Reportedly, COX-2 is involved with the pro-inflammatory activities and the level of COX-2 in synovial tissues is proportional with the severity of AS48. Notably, NS-398, as a new COX-2 inhibitor with anti-inflammatory activity that selectively inhibits prostaglandin G/H synthase/COX-2 activity in vitro, can result in the less gastrointestinal toxicity than standard NSAIDs49. DL-thiorphan is a neutral endopeptidase (also known as skin fibroblast elastase or membrane metalloendopeptidase) inhibitor that enhances cell proliferation in bone marrow. A previous study has indicated that the skin fibroblast elastase may increase the number of fibroblast50. Additionally, neutral endopeptidase may play a key role in turning off peptide signaling events at the cell surface which induce the abnormal of the immune and inflammatory systems51. Since that AS is closely associated with the fibroblast proliferation and inflammatory, it was supposed that NS-398 and DL-thiorphan might be a potential drug for AS treatment under INF-γ stimulation. However, there are no related studies exploring the effect of NS-398 and DL-thiorphan for AS treatment till now. Therefore, further validations are required to confirm the efficacy and safety of this drug in combination with other agents.
CONCLUSION
In conclusion, a total of 317 DEGs (eg., CDK1 and BUB1) and 394 DEGs (eg., PRKCQ and BRG1) were identified in the AS patients with and without INF-γ interference groups. CDK1 and BUB1 might be novel INF-γ sensitive response gene for the progression of AS via regulating the proliferation of fibroblast. PRKCQ and BRG1 may be closely associated with AS under the non-INF-γ stimulation via regulating T cell activation and inflammation. In addition, NS-398 as a new COX-2 inhibitor with anti-inflammatory activity and DL-thiorphan associated with the fibroblast proliferation might be a potential drug for AS treatment under INF-γ stimulation. However, further experiment is needed to confirm those speculations.
SIGNIFICANCE STATEMENT
This study discovered several novel INF-γ sensitive response DEGs like CDK1 and BUB1 and non-INF-γ sensitive DEGs PRKCQ and SMARCA4 between the AS patients and healthy controls that can be beneficial for explaining the mechanisms of AS progression. In addition, the predicted NS-398 and DL-thiorphan might be ideal small-molecule drugs for As treatment. This study will help the researcher to uncover the critical areas of AS that many researchers were not able to explore. Thus a new theory on AS treatment targets may be arrived at.
ACKNOWLEDGMENTS
This study was supported by Science and Technology Project of Shaanxi Social Development (Project number: 2016SF-223).