Novel Autophagy-Related Blood Biomarkers Associated with Immune Cell Infiltration in Ankylosing Spondylitis

Introduction

Ankylosing spondylitis (AS) is a chronic inflammatory bone disease that causes functional impairments and inflammation of the sacroiliac joints and spine. This results in progressive joint ankylosis and can seriously affect a patient’s quality of life.¹ Early diagnosis of AS is crucial to reduce the burden on the patient. Previous reports have indicated that AS is influenced by a combination of genetic, immune, environmental, and metabolic factors.^2–4 Despite this knowledge, the exact causes and development of the disease are still not fully understood. Furthermore, there is a shortage of effective treatments available for AS patients, and a complete cure is currently not possible. Additionally, defining the disease prognosis and progression in AS presents a challenge. Therefore, it is imperative to investigate biomarkers that can facilitate early detection of AS, as this can significantly enhance the prognosis of patients.

Autophagy is a critical component of cellular metabolism that facilitates the breakdown of organelles and cytoplasmic proteins into fatty acids and amino acids for energy recycling.⁵ Extensive research has shown that autophagy is implicated in the development of numerous diseases such as tumor, neurodegenerative diseases, central nervous system diseases, and autoimmune diseases.^6–8 Recent studies have uncovered that mitophagy and autophagy may also have a significant impact on bone diseases.^9–11 In addition, autophagy dysfunction plays a role in the development of AS.¹² Ankylosing spondylitis (AS) patients exhibit defective autophagy activity, which may contribute to spinal injury.¹³ However, the diagnostic usefulness of these autophagy-related genes (ARGs) in AS has not been fully elucidated.

In the field of bioinformatics, there has been notable progress in the identification of new targets and pathways for understanding the pathogenic processes of diseases, such as bone diseases.^14–16 This has been achieved through the use of high-throughput techniques and advanced bioinformatics analyses. As a result, this approach could be useful in investigating the role of autophagy in the development of AS. For this study, we obtained gene expression profiles from the Gene Expression Omnibus (GEO) database and used bioinformatics to identify autophagy-related diagnostic biomarkers. Our findings could aid in the identification of new diagnostic biomarkers for AS patients and enhance our understanding of the pathogenesis involved in AS.

Materials and Methods

Acquisition of Dataset and Data Processing

The human autophagy-dedicated database (HADb) was used to collect the ARGs. The GSE73754 dataset’s raw matrix files were downloaded from the GEO database and then extracted and normalized using the R package “affy” prior to data analysis. The GPL10558 platform annotation file was used to convert the probe expression matrix into the gene expression matrix.

Identification of Differentially Expressed Analysis of ARGs (DEARGs)

The limma package in R software was utilized to identify the differentially expressed genes (DEGs), with a significance threshold of p.adj < 0.05. A volcano plot of the DEGs was generated using the “ggplot2” package. The DEARGs were visualized in a heatmap using the “pheatmap” package.

Functional Enrichment Analysis

To explore the biological functions of the DEARGs, Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) enrichment analyses were conducted in this study. The enrichment analysis was performed using the Cluster Profiler package of R software, with statistical significance being defined as p < 0.05.

Construction of Protein-Protein Interaction Network (PPI) and Identification of Core Biomarkers

In this study, the PPI network of the DEARGs was constructed using the Search Tool for the Retrieval of Interacting Genes (STRING) database. The result of the PPI analysis was visualized with Cytoscape (version 3.8.2). Core genes were identified through the cytoHubba plug-in utilizing 11 algorithms.¹⁷ To evaluate the diagnostic efficacy of the core DEARGs, receiver operating characteristic (ROC) curve analysis was conducted. The pROC package in R software was utilized for this purpose.

Quantitative Real-Time Polymerase Chain Reaction (qRT-PCR)

Between 2021 and 2022, we collected 10 health control (HC) blood samples and 10 blood samples from patients with AS at Heilongjiang University of Traditional Chinese Medicine. Approval for this research was obtained from the Ethics Committee of Heilongjiang University of Traditional Chinese Medicine. The informed consent was obtained from the participants. The blood-derived Total RNA was isolated by utilizing the TRIzol reagent (Invitrogen, USA), followed by the synthesis of cDNA using the First-Strand Synthesis System Kit (Invitrogen, USA). The qRT-PCR was subsequently conducted with the utilization of SYBR green qPCR Master Mix (ThermoFisher Scientific, USA) on a Real-Time PCR System (ThermoFisher Scientific, USA). The primers were shown in Table S1.

Immune Infiltration Landscape

The xCell algorithm is a new technique for determining the relative abundance of immune cells based on gene signatures.¹⁸ In this study, xCell was used to compare the proportion of infiltrating immune cells between AS and HC and the results were presented in a box plot. In addition, the association between signature gene expression and immune cell infiltration was examined using the “ggplot2” package and presented in a lollipop plot.

Gene Set Enrichment Analysis (GSEA)

To perform single gene GSEA analysis, we used the GSEA software to classify the samples into two groups of high and low expression based on the median values of gene expression levels. To investigate the molecular mechanisms relevant to the genes based on the phenotypic grouping, we downloaded the subset “c2.cp.kegg.v7.4.symbols.gmt” from the Molecular Signatures Database.

Results

Identification of DEARGs in AS

Initially, a comprehensive analysis of the GSE73754 dataset using R software revealed a total of 2531 DEGs. To visualise these results, volcano plots were generated, where blue dots indicated 1421 genes with down-regulation, red dots represented 1110 genes with up-regulation and grey dots indicated genes that showed no significant difference (Figure 1A). A total of 29 DEARGs were then identified in the intersection of the ARG and DEG datasets using Venn analysis (Figure 1B). The mRNA expression levels of these 29 DEARGs were plotted on a heat map (Figure 1C).

Figure 1 Identification of DEARGs in AS. (A) Volcano plots of DEGs. (B) Venn diagram used to identify 29 DEARGs. (C) Heatmap pattern of the DEARGs.

Functional Enrichment Analysis

As shown in Figure 2A, GO enrichment analysis revealed that the DEARGs were significantly enriched in autophagy, process utilizing autophagic mechanism, macroautophagy, catabolic process, cellular catabolic process, positive regulation of protein metabolic process, etc. The results KEGG pathways revealed that the DEARGs were significantly enriched in the autophagy-animal, apoptosis, Chagas disease, mitophagy-animal, necroptosis, NOD-like receptor signaling pathway, human immunodeficiency virus 1 infection, etc (Figure 2B).

Figure 2 Functional enrichment analysis of 29 DEARGs. The results of GO (A) and KEGG (B) were shown by circle charts.

Identification of Core Biomarkers

The Cytoscape software was used to visualize the PPI of DEARGs (Figure 3A). By overlapping the top 10 DEARGs in 11 algorithms, two hub genes, HSPA5 and SQSTM1, were identified (Figure 3B). The expression of these genes was found to be significantly lower in the AS group compared to the HC group (Figure 4A). The diagnostic area under the curve (AUC) values for HSPA5 and SQSTM1 were 0.739 and 0.846, respectively (Figure 4B). Our findings indicated those genes have the diagnostic value to distinguish the AS from HC. In addition, we collected clinical blood samples to verify the expression of the core genes (Figure 5). The results were consistent with the bioinformatics results, with HSPA5 and SQSTM1 gene expression levels significantly downregulated in AS patients.

Figure 3 Identification of core DEARGs by PPI network analysis. (A) Protein-protein interaction networks of the 29 DEARGs. (B) Intersection diagram of top 10 genes based on 11 topological algorithms. Green outlined section indicates core genes identified through the use of 11 algorithms.

Figure 4 Analysis core DEARGs in GSE73754 dataset. (A) Gene expression level of HSPA5 and SQSTM1. (B) ROC curve of HSPA5 and SQSTM1. **p < 0.01 and ***p < 0.001.

Figure 5 Validation of core DEARGs by clinical blood samples. *p < 0.05 and ***p < 0.001.

Immune Cell Infiltration Analysis

Differences in immune cell infiltrates between the HC and AS groups were assessed using xCell. As depicted in Figure 6A, a notable reduction in the levels of basophils, CD8+ T cells, CD8+ Tcm, CD8+ Tem, macrophages M2, NK cells, Tgd cells, Th1 cells, and Th2 cells was observed in the AS group. Conversely, the level of neutrophils was significantly increased in the AS group. In addition, correlation analysis indicated that HSPA5 gene was positively correlated with Th1 cells, CD8+ T cells, CD4+ memory T cells, macrophages, CD8+ Tem, and plasma cells (Figure 6B). SQSTM1 gene was positively correlated with basophils (Figure 6C).

Figure 6 Immune infiltration analysis. (A) The landscape of the immune cells between HC and AS groups. (B and C) Correlation analysis between the hub DEARGs (HSPA5 and SQSTM1) and immune cell infiltration levels. *p < 0.05, **p < 0.01, and ***p < 0.001; ns indicates no statistical difference.

Single-Gene GSEA

We conducted an analysis of potential signaling pathways linked to the core genes by utilizing GSEA. As presented in Figure 7, cysteine and methionine metabolism, citrate cycle TCA cycle, lysine degradation, tyrosine metabolism, oxidative phosphorylation, and T cell receptor signaling pathway were mainly enriched in the HSPA5 high-expressed subgroup. The starch and sucrose metabolism and aldosterone regulation sodium reabsorption were enriched in SQSTM1 low-expressed subgroup (Figure 8).

Figure 7 Single-gene GSEA analysis based on the median value of HSPA5 gene expression level.

Figure 8 Single-gene GSEA analysis based on the median value of SQSTM1 gene expression level.

Discussion

Ankylosing spondylitis (AS) is a inflammatory condition that affects the sacroiliac and spinal joints, leading to dysfunction of both the immune and skeletal systems.¹⁹ Although AS was described 200 years ago, its pathogenesis is still not fully understood. Previous reports indicated that the pathogenesis of AS involves genetic susceptibility, immune dysfunction, intestinal mucosal inflammation, and environmental factors.^20,21 Autophagy has been shown to have an impact on the immune system and has therefore been implicated in certain rheumatic disease processes. There is increasing evidence to support the hypothesis that autophagy may also play a role in the pathophysiology of AS.²² As early diagnostic indicators are currently inadequate, patients with AS often miss the optimal window for treatment, leading to an unfavourable prognosis. Therefore, identifying unique diagnostic markers and characterising the immune cell infiltration pattern in AS has significant implications for improving the prognosis of AS patients.

In this study, we obtained the AS expression profile dataset from the GEO database and identified 29 DEARGs. GO enrichment analysis revealed that these DEARGs were primarily associated with autophagy, process utilizing autophagic mechanism, macroautophagy, catabolic process, etc. KEGG results showed that these DEARGs were significantly enriched in autophagy-animal, apoptosis, Chagas disease, necroptosis, NOD-like receptor signaling pathway, human immunodeficiency virus 1 infection, etc. These findings indicated that the immune response and autophagy may be a crucial factor in the development and progression of AS. A previous study highlighted the critical role played by TNFAIP3-DEPTOR complex-mediated autophagy in suppressing the immune response in AS disease.²³ The elevated systemic immune-inflammation index observed in AS patients appears to have potential value as a novel biomarker to assess disease activity in AS.²⁴ The findings of this research are in line with our analysis, indicating the accuracy of the study’s results. In addition, PPI network construction revealed HSPA5 and SQSTM1 as key hub genes. Analysis of clinical samples showed that the expression levels of these genes could effectively differentiate AS from normal samples. Furthermore, the ROC curves generated for these genes showed remarkable AUCs, highlighting their diagnostic significance for AS. This study therefore recognised the importance of further investigation of these genes. Heat-shock-protein family A member 5 (HSPA5), a heat shock protein, serves as a crucial regulator of the Unfolded Protein Response (UPR) and is known to interact with CFC1, an oncoprotein that is developmentally regulated and anchored by glycosylphosphatidylinositol.²⁵ HSPA5 can protect cells from ER stress and damage caused by reactive oxygen species.²⁶ Specifically, through N-terminal arginylation in the presence of cellular stressors such as reactive oxygen species (ROS), HSPA5 can regulate the activity of the SQSTM1/p62 macroautophagic receptor.²⁷ HSPA5 is crucial for cell survival as it helps clear misfolded proteins in lysosomes by interacting with SQSTM1.²⁷ In head and neck cancer, HSPA5 was found to negatively regulate lysosomal activity by ubiquitinating MUL1.²⁸ In addition, HSPA5 has been identified as a promising target for treating osteoporosis. The use of HA15 to target HSPA5 has been found to impede the progression of osteoporosis, making it a potential therapeutic option for postmenopausal osteoporosis.²⁹ The p62 protein, also known as sequestosome 1 (SQSTM1), functions as a scaffolding protein that binds to ubiquitin. It is often found together with aggregates of ubiquitinated proteins in a number of neurodegenerative diseases and proteinopathies affecting the liver.³⁰ The discovery of p62/SQSTM1 as an autophagy receptor in mammals has advanced our comprehension of macroautophagy as a discerning mechanism.^31,32 Impairment of the ability of SQSTM1 to transport its target cargo for degradation can potentially disrupt the balance of cellular signalling pathways and ultimately contribute to a range of cardiovascular and cardiometabolic disorders.³³ In addition, a previous study has suggested that animal-related factors play an important role in the development of Paget disease of bone and may interact with SQSTM1 mutations to influence disease severity.^34,35 On the basis of our own findings as well as previous research, we hypothesised that HSPA5 and SQSTM1 are likely to play an important role in the progression of AS.

AS is characterised as an inflammatory disease and a recent study found evidence of immune cell involvement in its progression.^36,37 Analysis of immune cell infiltration in our study showed that NK cells, neutrophils, CD8+ T cells, Tgd cells and Th1 cells are the main culprits in the development of AS. In addition, our findings suggested that HSPA5 is critical for most of the immune cells involved. A key component in the progression of AS may be the increase in neutrophil count.³⁸ The peripheral circulation of AS patients showed alterations in NK cell subsets as demonstrated by single-cell RNA sequencing.³⁹ According to a recent report, the frequency of CD8+ T cells was significantly reduced in blood samples from AS patients. This suggests that CD8+ T cells may indeed play a key role in the development of AS.⁴⁰ In AS patients, peripheral blood mononuclear cells showed reversible expression of T-regs under anti-TNF-α therapy.⁴¹ The changes in Th1/2 cytokine levels observed in people with AS may be indicative of their disease activity and prognosis.⁴² Additionally, HSPA5 gene was positively correlated with Th1 cells, CD8+ T cells, CD4+ memory T cells, macrophages, and CD8+ Tem. It is suggested that HSPA5 plays a role in upregulating immune cells, ultimately contributing to the onset and progression of AS. Further investigation is needed to fully understand the complex relationship between genes and immune cells.

Conclusion

In conclusion, this study analysed ARGs in AS and identified two ARGs with promising diagnostic potential. Furthermore, immune infiltrates were found to be less frequent in the AS group compared to the HC group, while HSPA5 showed a positive correlation with immune cell infiltrates. These findings suggest that the identified ARGs could serve as valuable indicators and provide insight into potential immunotherapy against AS.