Identification of a Novel Gene Expression Signature Associated with Amino Acid Metabolism (AAM) in Ankylosing Spondylitis (AS)

Introduction

Spinal arthropathy (SpA) is a heterogeneous immune mediated chronic inflammatory disease that mainly affects axial bones. It was divided into axial spinal arthritis (axSpA) and peripheral spinal arthritis (pSpA) based on the location of the disease.^1,2 Ankylosing spondylitis (AS) is the most common type of axSpA, which can lead to the development of spinal rigidity and irreversible structural damage.³ At present, our measures to identify AS are still limited to clinical features, including imaging (radiology, computed tomography (CT), and magnetic resonance imaging (MRI)), HLA-B27, and C-reactive protein (CRP), as the main indicators for its clinical diagnosis.^1,4 However, even with modern diagnostic methods, the diagnostic sensitivity and specificity of AS are not higher than about 80%.⁵ Therefore, people eagerly look forward to more highly sensitive and specific methods to find more accurate biomarkers for the diagnosis of AS.

According to the researches, amino acid metabolism plays an important role in AS, especially the level of tryptophan, which is one of the indicators of disease progression in AS patients. Berlinberg et al found that the microbiome metabolizes tryptophan into indole-3-acetate (IAA) and indole-3-acetate (I3Ald), which played an important role in the axSpA epithelial barrier and immune cell function through a joint analysis of multiple omics.⁶ Meanwhile, the pro-inflammatory cytokine interferon-γ could inducing enzyme indoleamine 2,3-dioxygenase (IDO) to convert tryptophan into kynurine.⁷ Eryavuz Onmaz et al also demonstrated that pro-inflammatory cytokines can affect the expression of various kynurine pathway enzymes, leading to metabolic changes, affecting inflammatory response and immunity. They also used the Kynurine/Tryptophan ratio as a diagnostic marker for AS.⁸ In addition, Gao et al also found that the levels of urea, proline, and phenylalanine in the serum of AS patients were significantly increased.⁷ According to reports, the immune effects that occurred in collagen tissue may lead to excessive breakdown of collagen components, thereby increasing the circulating proline in the blood.⁹ It is also worth noticing that HLA-B27 molecules have been proven to be rich in proline,¹⁰ and the HLA-B27 positivity has been widely used in clinical practice to identify SpA.^11,12 What’s more, unpaired cysteine residues in HLA antigens can be directly or indirectly modified by homocysteine and then subsequently disrupted by abnormal autoimmune reactions, which is also considered a mechanism of AS associated with HLA-B27.¹³ These studies have all pointed out the important role of amino acid metabolism in AS. However, it is currently unclear how amino acid metabolism regulates the progression of AS and what role amino acid metabolism-related genes (AAMRGs) play in it.

In this study, we utilized the public database data related to AS and the bio-informatics analysis methods to explore the amino acid metabolism-related genes with diagnostic value for AS, providing a theoretical basis for further revealing the pathogenesis of AS and deepening our understanding of AS.

Materials and Methods

Data Sources

AS-related datasets (GSE25101, GSE73754, and GSE18781) were downloaded from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/). The samples of these datasets are whole blood samples. The GSE25101 dataset consists of 16 normal and 16 AS samples. The GSE73754 dataset includes 20 normal samples and 52 AS samples. The GSE18781 dataset contains 18 AS and 25 normal samples and it was used as external validation set. In addition, 374 amino acid metabolism-related genes (AAMRGs) were downloaded from reaction metabolism of amino acid derivatives of Molecular Signature Database (MSigDB) (http://www.broad.mit.edu/gsea/msigdb/).

Identification of Differential Expression Genes (DEGs) Between Normal and AS Groups

Differential analysis was performed on AS and normal groups in the GSE25101 and GSE73754 datasets using the “limma” R package (version 3.46.0), respectively.¹⁴ The screening conditions for DEGs were P < 0.05 and |log₂FC| > 0.5. The DEGs of the GSE25101 and GSE73754 datasets were noted as DEGs1 and DEGs2, separately. To understand the distribution of DEGs overall, the volcano map and the heat map of up- and down-regulated DEGs were plotted by “ggplot2” R package (version 3.3.6) and “pheatmap” R package, respectively. The intersection of the DRGs1 with the DEGs2 was taken to acquire common DEGs (Co-DEGs) via Veen tool.

Identification of Amino Acid Metabolism-Related DEGs (AAMR DEGs)

Firstly, based on AAMRGs, amino acid metabolism scores were calculated for all samples of the GSE25101 dataset by single sample gene set enrichment analysis (ssGSEA). Then, weighted gene co-expression network analysis (WGCNA) was used to identify AS and amino acid metabolism score-related genes (AS-AAMSRGs).¹⁵ The samples were clustered to remove outliers. And in order to build a scale-free network, the determination of the optimal soft threshold (β) was performed. Thereafter, the adjacency matrix was established and converted to a topological overlap matrix (TOM), and the gene dendrogram and module color were established using the degree of dissimilarity. Correlations between modules and traits (AS, normal and amino acid metabolism score) were then calculated. The modules that met both |cor| > 0.4 and P < 0.05 were selected for subsequent analysis. Next, AAMR DEGs were obtained through the intersection of co-DEGs and AS-AAMSRGs. Finally, we adopted GeneMANIA (http://genemania.org/) to construct the network of AAMR DEGs interactive proteins based on physical interactions, co-expression, predicted, co-localization, pathway, genetic interactions, and shared protein domains.¹⁶

Identification of Biomarkers

Firstly, LASSO was utilized to identified the biomarkers. Then, the diagnostic value of the LASSO model were analyzed via receiver operating characteristic (ROC) curves in GSE25101, GSE73754, and GSE18781 using the “pROC” R package.¹⁷

Creation of a Nomogram

The nomogram was created via “Rms” R package to demonstrate the efficacy of biomarkers for the discrimination of AS. To verify the accuracy of the nomogram, the corresponding calibration curve, clinical impact curve, and decision curve analysis (DCA) were plotted.

Gene Set Enrichment Analysis (GSEA)

The samples of the GSE25101 were sorted into high and low expressed groups according to the median value of the expression of each biomarker, and differential analysis was done. All genes were ranked via log fold change (logFC). The GSEA were performed with |NES| > 1, P < 0.05, and FDR < 0.25.

Immune Analysis

The ssGSEA algorithm was employed to compute the infiltration abundance of 28 immune cells in AS and normal groups in the GSE25101.¹⁸ Differences between AS and normal groups were assessed by Wilcox test. Moreover, the relevance of biomarkers with differential immune cells were analyzed by Spearman algorithm.

Created of a miRNA-mRNA Network

To investigate the mechanisms involved in biomarkers in AS, we created a miRNA-mRNA network. Firstly, we used Mirwalk database (http://mirwalk.umm.uni-heidelberg.de/), TargetScan database (http://www.targetscan.org/vert71/) and miRTarBase (http://mirtarbase.mbc.nctu.edu.tw/index.html) to predict miRNAs that might regulate biomarkers. Then, the miRNAs obtained via the three databases were taken to intersect to obtain the common miRNAs. Finally, Cytoscape software was employed to visualize the miRNA-mRNA network (version 3.7.1).¹⁹

Drug Prediction of Diagnostic Genes

Finally, each biomarker was used as a key word to search for drugs associated with AS in Drug-Gene Interaction Database (DGIdb, https://dgidb.genome.wustl). Afterward, the network diagram of interactions between each biomarker and its corresponding medicines was created via Cytoscape software (version 3.7.1).¹⁹

Expression Profiles of Biomarkers in Public Datasets

To further demonstrate the reliability of the results, we evaluated the expression of biomarkers between AS and normal groups in the GSE25101, GSE73754, and GSE18781.

Sample Collection and Real Time Quantitative-Polymerase Chain Reaction (RT-qPCR)

In our study, the blood samples of 5 AS patients (the disease group) and 5 healthy individuals (the control group) at the Third People’s Hospital of Chengdu, between February 2021 and February 2022, were collected. The study was approved by the Ethics Committee of the Third People’s Hospital of Chengdu and the IRB reference number was [2019]S-103.

The RNA was extracted by TRIzol (Ambion, Austin, USA) in accordance with the instructions provided by the manufacturer. The 2xUniversal Blue SYBR Green qPCR Master Mix (Servicebio, Wuhan, China) was employed for RT-qPCR. The internal reference gene was GAPDH. The 2^−ΔΔCt method was applied to compute the expression of biomarkers²⁰ (Supplementary Tables 1–3).

Results

There Were 3 Co-DEGs Between AS and Normal Groups

There were 115 DEGs (64 up-regulation genes and 51 down-regulation genes) between AS and normal groups in GSE25101 (Figure 1A and B). Similarly, 176 DEGs contained 74 up-regulation genes and 102 down-regulation genes in DEGs2 (Figure 1C and D). 3 Co-DEGs were obtained through the intersection of DEGs1 and DEGs2 (Figure 1E).

Figure 1 Identification of Co-differentially expressed genes (DEGs) between ankylosing spondylitis (AS) and normal groups. (A and B) Volcano plot for DEGs ((A) GSE25101 (B) GSE73754). (C and D) Heatmap of DEGs ((C) GSE25101 (D) GSE73754). (E) Venn diagram of amino acid metabolism related genes (AAMRGs) up-regulated and AAMRGs down-regulated.

Acquisition of AAMR DEGs Between AS and Normal Groups

There were no outliers in GSE25101 (Figure 2A). With regard to the WGCNA parameter setting, β was 9, and 25 modules were identified (Figure 2B and C). The yellow module was markedly relevant to AS and amino acid metabolism score (Figure 2D). And 8401 AS-AAMSRGs were obtained. Subsequently, there were 3 AAMR DEGs that were obtained by the intersection of Co-DEGs and AS-AAMSRGs, namely TP53INP1, RBM47, TUBB (Figure 2E). To explore the genes interactions of three AAMR DEGs, the gene–gene interaction network of AAMR DEGs was created (Figure 2F), and these genes demonstrated a complex interaction network.

Figure 2 Acquisition of AAMR DEGs between AS and normal groups. (A) Sample clustering to detect outliers in GSE25101.(B) Soft threshold estimation for scale-free corepresentation networks. (C) All differences AAMRGs cluster. (D) Correlation between each module in the AS or Normal Groups.(E) Differentially expressed genes and modular genes. (F) miRNA-hub gene Network.

There Were 3 Biomarkers

Furthermore, there were 3 biomarkers via LASSO algorithm. Moreover, the AUC of the LASSO model were all greater than 0.7, suggesting that the 3 biomarkers had diagnostic value for AS (Figure 3A–C). A nomogram including 3 biomarkers was created (Figure 3D). The adjusted C-index was 0.875 (Figure 3E). The DCA curve suggested the nomogram had a higher net benefit (Figure 3F). The clinical impact curve revealed that the “number high risk” curve was nearly identical to the “number high risk with event” curve when the high-risk threshold was between 0.6 and 1.0 (Figure 3G). These results demonstrated that the nomogram had high a predictive accuracy for the AS progression.

Figure 3 Identification of biomarkers. (A–C) Receiver operating characteristic (ROC) for the 3 AAMR DEGs ((A) GSE25101, (B) GSE73754, (C) GSE18781). (D) Nomogram based on characteristic genes. (E) Calibration curve. (F) Decision curve analysis (DCA) curve. (G) Clinical impact curve.

Enrichment Analysis of Biomarkers

To explore potential biological mechanisms associated with biomarkers, we performed GSEA. RBM47 was enriched to 2097 GO terms and 73 KEGG pathways (Figure 4A and B), such as neutrophil activation, neutrophil activation involved in immune response, chemokine signaling pathway. There were 1622 GO terms and 73 KEGG pathways for TP53INP1 (Figure 4C and D), for instance, neutrophil degranulation, neutrophil mediated immunity, growth hormone synthesis, secretion and action. And TUPP was enriched to 914 GO terms and 84 KEGG pathways (Figure 4E and F), such as cotranslational protein targeting to membrane, ribosome.

Figure 4 Gene set enrichment analysis (GSEA) of three biomarkers. (A and B) RBM47 (C and D) TP53INP1 (E and F) TUBB. (A, C and E) Gene Ontology (GO). (B, D and F) Kyoto Encyclopedia of Genes and Genomes (KEGG).

Immuno-Infiltration Analysis Between AS and Normal Groups

In total, the infiltration abundance of 8 immune cells between AS and normal groups was considerably different, namely CD56 bright natural killer cell, activated dendritic cell, effector memory CD4 T cell, etc. (Figure 5A). The correlation analysis revealed that TP53INP1 and RBM47 were significantly positively correlated with gamma delta T cell, effector memory CD4 T cells, and type 17 T helper cell (Figure 5B). TP53INP1 and RBM47 were negatively relevant to natural killer T cell and effector memory CD8 T cell, and RBM47 was also negatively related to MDSC (Figure 5B). The relevance of TUBB to effector memory CD8 T cells, natural killer T cell, and MDSC was positive, and the relevance TUBB to activated dendritic cell, type 17 T helper cell, gamma delta T cell, and effector memory CD4 T cell was opposite (Figure 5B). Interestingly, none of the three biomarkers relevant to CD56 bright natural killer cell (Figure 5B).

Figure 5 Immuno-infiltration analysis between AS and normal groups. (A) Differences of immune cells between AS and normal groups were analyzed by Wilcox. Test and the boxplot. (B) Heatmap of correlations of diagnostic genes with differential immune cells. *P < 0.05, **P < 0.01, ***P < 0.001. Abbreviation: ns, not significant.

The miRNA-mRNA Networks

In total, 27 miRNAs were obtained by prediction. The miRNA-mRNA network was created (Figure 6). In the miRNA-mRNA regulatory networks, the 3 miRNAs (has-miR-429, has-miR-200b-3p, has-miR-200c-3p) could regulate TP53INP1 and TUBB. There was only 1 miRNA (has-miR-122-5p) could regulate RBM47.

Figure 6 The miRNA-mRNA regulatory networks.

Biomarker-Based Potential Drug Prediction

We constructed a drug–gene interaction network from the DGIdb. The results showed that 51 target drugs were related to TUBB (Figure 7), such as colchicine, vinblastine, vincristine.

Figure 7 The drug–gene interaction network from the DGIdb.

The Verification of Biomarkers

In all public datasets, there were considerably different in expression of all biomarkers between AS and normal groups (Figure 8A–C). TP53INP1 and RBM47 were highly expressed and TUBB was lowly expressed in the AS group in GSE25101, GSE73754 and GSE18781. Additionally, to further validate the expression of biomarkers, we performed RT-qPCR. The results revealed that RBM47 and TP53INP1 were considerably highly expressed in the disease group (Figure 8D and E), and the expression trends were in line with other findings based on open-access datasets. However, the expression trend of TUBB was opposite to that of public databases, possibly due to heterogeneity of the samples, as well as variations in the experimental design (Figure 8F).

Figure 8 Expression profiles of diagnostic genes in external validation datasets ((A) GSE25101, (B) GSE73754, (C) GSE18781). The expression of three biomarkers by RT-qPCR ((D) RBM47, (E) TP53INP1, (F) TUBB). *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001. Abbreviation: ns, not significant.

Discussion

Amino acid metabolism has been shown to play an important role in the AS process, which may be related to inflammatory reactions and immunity,⁸ but the specific pathways and mechanisms of action were still unclear. In this study, we obtained three AS diagnostic genes related amino acid metabolism, RBM47, TP53INP1 and TUBB, and they were involved in oxidative phosphorylation, osteoclast differentiation, and a variety of immune or inflammatory related signaling pathways, such as B cell receptor signaling pathway, T cell receptor signaling pathway, chemokine signaling pathway, and neutrophil activation, and etc.

Oxidative phosphorylation could regulate cellular processes, affect intracellular metabolic reactions, and promote the normal resolution of inflammatory reactions, through affecting protein activity. Recent studies by Najafzadeh and Lari have shown significant enrichment of the oxidative phosphorylation pathway in AS patients, which confirmed the importance of oxidative phosphorylation in AS.^21,22 Tumor protein 53-induced nuclear protein 1 (TP53INP1) is one of the main target genes that mediate the p53 signaling pathway of antioxidant activity, and is significantly correlated with oxidative physiology.²³ Perhaps, TP53INP1 exerts antioxidant activity by affecting oxidative physiology, affecting the AS process. In addition, TP53INP1 had been shown to be overexpressed in stress responses including inflammation,²⁴ which was consistent with our research findings.

The pathological changes of AS mainly include osteolytic bone destruction and osteogenic changes. The interaction of mesenchymal stem cells, osteoblasts, osteoclast and other bone cells with T cells, B cells and other immune cells jointly regulated the pathogenesis of AS.²⁵ Piccolo et al pointed out that vertebrate trunk induction requires inhibition of bone morphogenetic protein (BMP) signaling.²⁶ RBM47 is an RNA-binding motif protein that plays multiple roles in bone formation and embryonic development.²⁷ Guan et al also pointed out that RBM47 is involved in bone formation and is associated with BMP signaling and Wnt signaling pathways.²⁶ Our study found a significant positive correlation between RBM47 and osteoclast differentiation function, which may suggest that RBM47 may promote osteoclast differentiation through BMP signaling, thereby affecting AS progression. It is worth noting that RBM47 is also annotated in the B cell receiver signaling path. Moreover, Wei et al also proposed that RBM47 can enhance the production of IL-10 and promote immune suppression of B cells,²⁸ which further deepens our understanding of AS bone immunology. RBM47 can also interact with B cells by participating in the B cell receptor signaling pathway, and jointly regulate the pathogenesis of AS.

TUBB is a type of protein coding gene of β-Tubulin protein. It has been shown to be significantly correlated with the development of the innate immune system²⁹ and differentially expressed in rheumatoid arthritis.³⁰ Besides, Zhang et al also found differential expression of TUBB-peptide in the plasma of AS patients. And in vitro cell experiments showed that TUBB significantly increased the proliferation of fibroblasts,³¹ while fibroblasts may develop into osteoblasts. Numerous studies have shown that fibroblasts are the starting point for AS ligament ossification.^32,33 Moreover, Joos et al confirmed IL-1β is significantly negatively correlated with TUBB, and participates in regulating various cytoskeleton components in human chondrocytes, and is related to the pathogenesis of osteoarthritis.³⁴ We have also mentioned that pro-inflammatory cytokines (such as TNF-α) can affect the expression of various kynurenine pathway enzymes, leading to metabolic changes and affecting inflammatory response and immunity.⁸ Excitingly, when we further analyzed the upstream regulatory mechanisms of diagnostic genes, we found that has-miR-429, has miR-200c-3p, and has miR-200b-3p can simultaneously regulate TP53INP1 and TUBB. Among them, has miR-200b-3p has been reported to be involved in epithelial mesenchymal transition (EMT).³⁵ Research has shown that cysteine rich domains 1 could promote EMT induced by IL-1β.³⁶ These results might suggest that has-miR-200b-3p is involved in EMT induced by IL-1β to regulates TUBB, thereby affecting the expression of various kynurenine pathway enzymes, and participates in the regulation of cytoskeleton components, which may also be the key mechanism of amino acid metabolism involved in AS. In future drug predictions, we predicted that TUBB has 51 targeted drugs, such as colchicine, vinblastine, vincristine, etc. Among them, colchicine is a common anti IL-1 and TNF-α inhibitors, and have been applied in the clinical treatment of AS.³⁷ Our study further confirms the clinical value of TUBB and provides a theoretical basis for the clinical treatment of colchicine in AS.

We further explored the immune microenvironment of AS, and the correlation results showed that the three diagnostic genes were highly significantly related to gamma delta T cell and type 17 T helper cell. According to reports, gamma delta T cells that produce IL-17 promote bone regeneration by acting on mesenchymal stem cells during fracture healing,³⁸ and the response seems to depend on nonantigenic driving mechanisms, including cytokines (IL-1, IL-23, and IL-28) and chemokines (Th1 and Th2).³⁹ Hull pointed out that the increase in type 17 T helper cells with IL-17 phenotype play a key role in AS treatment, and this mechanism is associated with chemokines Th1 and Th2.⁴⁰ Our study not only pointed out that the functions of all three diagnostic genes were annotated in the chemokine signaling pathway, but also found that all three diagnostic genes were associated with NF-κB signaling pathway. What is interesting is that the balance of Th1/Th2 is maintained by NF-κB signal.⁴¹ These evidence all suggested that RBM47, TP53INP1, and TUBB regulate the balance of Th1/Th2 through the chemokine signaling pathway or NF-κB signaling pathway, and participate in the regulation of bone regeneration in AS by gamma delta T cells and type 17 T helper cells.

In summary, this study screened three diagnostic genes for AS related to amino acid metabolism, and further understood their molecular mechanisms of action on AS through functional enrichment and molecular regulation prediction. However, our conjecture about the regulatory mechanism of diagnostic genes has not been confirmed yet. We will further verify the mechanism of has-miR-200b-3p regulating TUBB through cell experiments, and explore the roles of pro-inflammatory factors such as IL-1, IL-17, IL-1β, and TNF-α in it.