Introduction

Fibromyalgia (FM) is a pervasive clinical syndrome that features chronic extensive musculoskeletal pain, usually along with fatigue, sleep trouble, and cognitive disorder.1,2 The global average prevalence of FM was 2.7% and varied from 0.4% to 9.3% in different countries, and the morbidity in females is approximately three times that in men.3 Not all physicians have a clear cognition of the clinical features of FM, thus resulting in a missed diagnosis. Previous studies suggested peripheral abnormalities and central sensitization were recognized as pivotal contributors to FM, and other factors such as inflammation, immunity, genetics, endocrine, etc. may also participate in its progression.4 Ankylosing spondylitis (AS) is a chronic inflammatory arthritis, and the main lesion sites were sacroiliac joints, spine, and peripheral joints, some severe cases may occur with spinal deformity and rigidity.5 A nationwide study involving 17,088 participants suggested patients with AS were linked to a higher FM risk, especially those older than 65.6 Notably, the prevalence of FM in all AS patients was 4.11%, while the prevalence of FM among female patients with AS increased to 10.83%.7 Another study involving 419 patients stated that the estimated prevalence of FM in AS was 12.7%.8 The current study supports the genetic underpinnings in the development of FM and AS,1,9 however, the underlying mechanism and pathophysiology link in these two diseases warrants further investigation.

Human leukocyte antigen B27 (HLA-B27) was discovered in 1973 and undoubtedly become the dominating biomarker for diagnosis of AS.5,10 However, HLA-B27 lacks specificity and sensitivity in AS patients complicated by other diseases, such as inflammatory bowel disease.11,12 Improved biomarkers are clinically necessary because we can employ them to measure and evaluate the process of FM and as an indication of early therapeutic intervention.13 Further, concomitant FM is a highlighted problem in AS patients, and its presence often lacks early diagnosis. A previous study reported plasma Pentraxin −3 (PTX-3) could serve as a reliable marker in improving diagnostic certainty of FM in the presence of AS.14 However, the sensitivity and specificity of plasma PTX-3 are insufficient. Therefore, finding more sensitive and specific biomarkers, as well as implementing earlier drug interventions for concomitant FM in AS patients, are of great clinical benefit.

We aimed to screen possible diagnostic biomarkers and unravel shared pathways from the expression profiling datasets in peripheral whole blood in AS and FM patients from the GEO database. To reach this goal, DEG analysis and WGCNA were first conducted, and intersected genes were primarily identified. Meanwhile, KEGG and GO enrichment analyses were implemented to explore the shared pathway in AS and FM. Machine learning algorithms were further implemented to filter out two hub genes with good diagnostic values and validated in clinical samples in two diseases: CETN3 and CACNA1E. Additionally, the ssGSEA method was employed to explore immune cell infiltration in two diseases and correlation analysis between immune cells and hub genes. Ultimately, molecular docking illustrated a potential drug that targets the hub gene of two diseases. Overall, our study offers worthy insight into the underlying mechanisms behind AS and FM and emphasizes the potential of CETN3 and CACNA1E as diagnostic markers for these circumstances.

Method Datasets Information

As is fully delineated in Figure 1, the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/) was obtained to search for expression profiling of AS and FM. Notably, the following search strategies were adopted to ensure the consistency and accuracy of our results: 1) keywords “fibromyalgia” and “ankylosing spondylitis” were used; 2) the samples extracted for gene expression profiles must be peripheral mononuclear cells (PBMCs); 3) the sample size of each group less than 10 were excluded to guarantee the precision of WGCNA.

Figure 1 The working flow chart of this study.

Abbreviations: AS, Ankylosing spondylitis; FM, Fibromyalgia; DEGs, differentially expressed genes; WGCNA, Weighted Genes Correlation Network Analysis; ROC, Receiver Operating Characteristic Curve.

Ultimately, the accession numbers of GSE25101, GSE73754, GSE67311, and GSE221921 were acquired. GSE25101 concluded 32 samples, of which 16 cases were from normal healthy and 16 cases from AS patients. GSE73754 contained 72 samples, comprising 20 healthy controls and 52 samples from AS patients. AS for datasets concerning fibromyalgia, GSE221921 included 189 peripheral blood samples, 93 samples from healthy controls, and 96 samples from FM patients. Additionally, GSE67311 consists of peripheral blood samples from 68 normal controls and 61 FM patients. The platform and grouping information about the datasets brought into this study was exhibited in Supplementary Table 1.

Screening of Differentially Expressed Genes (DEGs)

The selection of DEGs between the AS and healthy controls, as well as FM patients and their controls, was implemented using the “limma” package in R software. Therefore, DEGs in AS datasets (GSE25101 and GSE73754) and FM datasets (GSE67311 and GSE221921) were filtered upon the screening criteria of p-value < 0.05. Afterward, the expression patterns of DEGs were visualized via the heatmap and volcano plot. Common DEGs of AS and FM were obtained and visualized using the Venn plot employing the Draw Venn Diagram (https://bioinformatics.psb.ugent.be/webtools/Venn/).

Enrichment Analysis of Shared Genes from as with FM

To unravel the potential roles and functions among shared genes from ankylosing spondylitis with fibromyalgia, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were utilized. During the enrichment analysis, “org.Hs.eg.db”, “clusterProfiler”, “enrichplot” and “ggplot2” packages were employed in R software. GO analysis revealed shared gene-related biological processes (BP), molecular function (MF), and cellular components (CC). Meanwhile, KEGG provided the enriched signaling pathways demonstrated by the targeted genes. Afterward, the results were visualized using the bar plot regarded the threshold of p-value < 0.05.

Weighted Gene Co-Expression Network Analysis (WGCNA)

To construct the co-expression network and screen gene module strongly correlated with AS and FM, WGCNA was conducted on GSE25101, GSE73754, GSE221921, and GSE67311 using the “WGCNA” package (version 1.72–5) of R software.15 Initially, hierarchical cluster analysis was carried on by employing the “hclust” function in R software to remove outlier samples. Afterward, based on a scale-free network, the “pickSoftThreshold” function was performed to choose the suitable soft threshold (power, from 1 to 20) for automatic network construction. Consequently, the suitable soft threshold (β=9 for GSE25101 and β=4 for GSE67311). Then, genes with similar patterns were gathered into a module employing the function “dynamic tree cutting” when setting the minimum module size at 30. Next, after obtaining the modules, the different modules eigengenes (ME) were acquired in accordance with the first principal component of the module expression, while the module-trait relationships were assessed complying with the linkage between MEs and clinical phenotypes. Eventually, the modules with the most remarkable positive and negative correlations of module-trait relationships were identified.

Machine Learning

For further refining the filtered candidate genes correlated with FM in AS patients, we utilized three machine learning algorithms. The Least Absolute Shrinkage and Selection Operator (LASSO), is a well-known linear prediction method, that works mainly regarded regression coefficients and provides widespread use in a variety of fields.16 In detail, the “glmnet” package was employed to operate LASSO. cv.glmnet was exploited to optimize lambda. For the parameters, the scale of “λ” was set ranging from 0 to 100 with “binomial” and “class”. Under the minimum λ, glmnet was processed to the LASSO with alpha and a “binomial” method.

Afterward, the support vector machine-recursive feature elimination (SVM-REF) algorithm, widely applied in the field of machine learning, adopts the principle of minimizing structural risk and minimizing empirical error to improve learning performance. SVM was implemented using the “e1071” package (version 1.7.12). tune.svm was adopted to improve the settings parameter with the kernel of “linear”, and the cost ranging from 1 to 20. Following, compliance with the optimal number of support vectors, the SVM model was completed.

Random forest (RF) is an algorithm that integrates multiple trees through the idea of ensemble learning. Its basic unit is the decision tree, which has the advantages of high efficiency of decision tree and accuracy of ensemble learning and has strong model generalization ability. RF was carried out with the “randomForest” package (version 4.7–1.1). Firstly, the tuneRF function was employed to optimize 0–700 trees with one step size. RF was managed according to the minimum error rate to acquire the optimal tree number.

Building of Nomogram and The Assessment of The Prediction Model

The clinically applicable nomogram was developed with the two selected hub genes (CACNA1E and CETN3) by employing the “rms” package. The calibration curves were established to evaluate the predictive power of the nomogram in FM-related AS.

Single-Sample Gene Set Enrichment Analysis (ssGSEA)

To make a thorough exploration of immune cell infiltration states in both AS and FM patients, the “GSVA” function in R software was employed to calculate the immune cell infiltration score of each sample using the gene set composed of immune cell markers.17 Plus, a Spearman correlation method was executed to identify the correlation between the key genes and the immune cells.

Prediction of Therapeutic Drug with Hub Genes via Molecular Docking

Molecular docking experiments were further implemented to explore the interaction between potential therapeutic drugs and the screened hub genes. Firstly, the Drug Signatures Database (DSigDB) on the Enrichr website (https://maayanlab.cloud/Enrichr/) was accessed to acquire potential drugs toward core genes.18,19 Afterward, the protein structures corresponding to the hub genes and small molecular drugs were obtained from the PDB database (https://www.rcsb.org/) and PubChem website (https://pubchem.ncbi.nlm.nih.gov/), respectively. Ultimately, the AutoDuck Vine (version 1.5.6) (https://ccsb.scripps.edu/mgltools/downloads/) and PyMOL software (version 2.5.5) were adopted to simulate the binding of targeted proteins and small molecular drugs.

Collection of Clinical Samples

The clinical trial was implemented strictly complying with the Declaration of Helsinki and obtained approval from the Institute Research Ethics Committee at Shenzhen Nanshan People’s Hospital to make sure there were no ethical issues. From January 1st, 2022 to December 31st, 2023, nine inpatients (including 5 males and 4 females, aged 36.05±3.21 years) clinically diagnosed with AS and FM were enrolled in this study. The AS diagnoses were executed according to the criteria developed by The Assessment of SpondylArthritis International Society (ASAS) in 2009.20 Furthermore, subjects were identified with FM if they satisfied the following criteria: i) Patients suffered chronic pain (> 3 months); ii) pain occurred in ≥4 sites of the body; iii) with fatigue, sleep disorder, and cognitive or somatic symptoms.1 Of note, patients diagnosed with cancer, infection, endocrine system disease, metabolic disorder, etc were ruled out. As a control group, 10 volunteers (including 5 males and 5 females, aged 32.86±4.47 years) undergoing routine health checkups were brought into the study. For each subject, 8 mL of blood was collected from the discarded blood samples after the routine blood test was accomplished. Afterward, the blood samples were labeled for subsequent experiments.

Validation of Hub Gene Expression Between Healthy Control and as-FM Groups via RT-PCR

Total RNA was extracted from whole blood samples strictly following the PAXgene protocol,21 and then reverse transcribed with a reverse transcription kit (Takara, Shija, Japan). Quantitative reverse transcription PCR was executed employing TB Green® Premix Ex Taq™ (Takara, Japan) complying with the instructions. The 2−ΔΔCt method was employed to calculate the relative gene expression values,22 and GAPDH served as the internal reference. The primer sequences are exhibited in Supplementary Table 2.

Statistics

The statistical analysis work of this study was achieved via the R software (Version 4.0.3) and the GraphPad Prism (Version 9.3.2). G power software (version 3.1.9.7, University of Dusseldorf, Germany) was employed to estimate the sample sizes, and statistical power (power,1-β) was set at 0.8. Statistically significant distinctions between the two groups were confirmed through unpaired Student’s t-test or Wilcoxon test. Correlations between hub genes and immune cells were determined by Spearman correlation analysis. Distinctions were regarded as significant at p < 0.05.

Results Data Treatment and Study Procedure

The brief research procedure is illustrated in Figure 1. Two raw datasets of AS and healthy controls comprising the discovery dataset GSE25101 and validation dataset GSE73754, as well as two raw datasets of FM and control samples including the discovery dataset GSE221921 and validation group GSE67311, were obtained from GEO databases (Supplementary Table 1). Each dataset was normalized respectively and then subjected to subsequent analysis.

Filtration of DEGs from the Discovery and Validation Cohorts of as and FM

For acquiring as many common genes as possible in peripheral blood samples of AS and FM patients, the DEGs were recognized with cut-off criteria of p-value < 0.05. For GSE25101, 2546 DEGs were identified with 1449 down-regulated and 1097 up-regulated genes (Figure 2A and C). Besides, 6220 DEGs of which 3401 down-regulated and 2819 up-regulated were screened in GSE73754 (Supplementary Figure 1A and C). For FM datasets, GSE221921 identified 9722 DEGs including 4608 down-regulated and 5114 up-regulated genes (Figure 2B and D). Besides, there were 1348 DEGs confirmed in GSE67311, of which 644 were down-regulated and 704 were up-regulated (Supplementary Figure 1B and D). Afterward, the up-regulated and down-regulated DEGs of these four datasets were intersected, respectively. Eventually, 179 shared down-regulated and 419 shared up-regulated DEGs in AS-FM patients were identified (Figure 2E and F).

Figure 2 Shared DEGs in AS and FM and functional analysis. (A and B) Illustration of DEGs in GSE25101 (AS) and GSE221921 (FM) using the heatmap. (C and D) Exhibition of DEGs in GSE25101 (AS) and GSE221921 (FM) using volcano plots. (E and F) The shared down-regulated and up-regulated DEGs between AS and FM. AS1, GSE25101; AS2, GSE73754; FM1, GSE221921; FM2, GSE67311.

Functional Enrichment Analysis

We then conducted GO and KEGG enrichment analysis on the DEGs from the above results. The finding unraveled that 419 shared up-regulated DEGs were principally enriched in terms of positive regulation of the cell cycle process and mitophagy (Figure 3A and B). Likewise, 179 shared down-regulated DEGs primarily enriched in BP terms like negative regulation of phosphate metabolic process, phosphorus metabolic process, and protein modification process. For CC, they were mainly clustered in the mitochondrial inner membrane, coated vesicle, and coated vesicle membrane. For MF, protein-folding chaperone binding and phosphoprotein binding were chiefly identified. Meanwhile, we observed pathways of neurodegeneration-multiple diseases, Parkinson’s disease, Prion disease, and Huntington’s disease were recognized via the KEGG analysis (Figure 3C and D).

Figure 3 Enrichment analysis. (A and B) GO and KEGG enrichment analysis of 419 intersecting up-regulated DEGs in AS and FM. (C and D) GO and KEGG enrichment analysis of 181 intersecting down-regulated DEGs in AS and FM.

Establishment of Co-Expression Network and Deciphering Hub Gene Modules

From the peripheral blood sample expression profile data, the WGCNA network was then built to investigate the connected co-expression network in AS and FM patients. In the discovery cohort of AS, the most suitable soft threshold power (β) was set at 9 based on scale independence and mean connectivity (Figure 4A and C). Similar gene modules were then merged and thirteen modules were recognized, with the MEgreen (r=0.42, p=0.02) and MEbrown (r=0.41, p=0.02) modules exhibiting a significant positive link to AS (Figure 4E). Analogously, in AS validation dataset, was analyzed using the same procedure. The soft threshold power was determined at 4 and the MEblue module (r=0.49, p=1e-05) was identified as the positive-related module (Supplementary Figure 2A, C and E).

Figure 4 Weighted genes correlation network analysis (WGCNA) of AS and FM. (A) Determination of soft threshold power for GSE25101 (AS). (B) Calculation of soft threshold power for GSE221921 (FM). (C) Cluster dendrogram of AS highly connected genes in key modules. (D) Cluster dendrogram of FM highly connected genes in key modules. (E) Relationships between modules and traits in AS. (F) Module-trait relationships in FM.

In the discovery cohort of FM, the suited soft threshold power was confirmed at 6. Afterward, ten gene modules were recognized, of which the MEmidnightblue module (r= 0.3, p=3e-05) and MEmagenta module (r=0.15, p=0.04) were confirmed as FM-related gene modules (Figure 4B, D and F). As for the validation group, the soft threshold power was set at 6, and seven gene modules were determined, of which the MEyellow module was recognized as an FM-related gene module (Supplementary Figure 2B, D and F).

Identification of Shared Hub Genes via Machine Learning Algorithm

The intersection of genes linked to AS and FM was recognized by comparing the genes that appeared in positively related modules for both diseases. A total of 143 common genes were discovered across six relevant modules in AS and FM by employing the Venn diagram (Figure 5A). Then we used WGCNA-derived significant module genes and DEGs yielded 6 potential hub genes for AS and FM (Figure 5B). Following, machine learning algorithms were adopted to screen underlying diagnostic biomarkers from 6 shared key genes previously determined in both AS and FM. In the discovery cohort of AS, the LASSO algorithm recognized two biomarkers as depicted in Figure 5C.

Figure 5 Machine learning in screening candidate diagnostic markers for AS with FM. (A) Venn plot of the shared genes between the MEblue, MEbrown, and MEgreen modules of AS and MEmidnightblue, MEmagenta, and MEyellow modules of FM by overlapping them. (B) Venn diagram of 6 shared genes identified from the intersection of up-regulated DEGs and significant positive modules in AS and FM. (C) Biomarkers of AS screening in the LASSO model. (D) Four crosstalk genes were selected by using the SVM-RFE algorithm for AS. (E) The error rate confidence intervals and the relative importance of genes for random forest algorithm in AS. (F) The Venn diagram of the intersection of LASSO, SVM-REF, and random forest signature genes. (G) Biomarkers of FM screening in the LASSO model. (H) Six crosstalk genes were selected by using the SVM-RFE algorithm for AS. (I) The error rate confidence intervals and the relative importance of genes for random forest algorithm in FM. (J) The Venn diagram of the intersection of LASSO, SVM-REF, and random forest signature genes.

By contrast, the SVM-REF algorithm illustrated that the model involved four genes reached the best accuracy (Figure 5D). Furthermore, the RF algorithm identified five underlying biomarkers and ranked them by their significance, as shown in Figure 5E. By intersecting the results of three algorithms, CETN3 and CACNA1E were affirmed as biomarkers for AS (Figure 5F).

The machine-learning algorithms were then conducted to discern possible biomarkers in FM. In the discovery group, the LASSO algorithm figured out six hub genes (Figure 5G). In the meantime, the SVM-REF method exhibited six genes that achieved the highest accuracy (Figure 5H). Eventually, the RF algorithm decided top five biomarkers according to their importance (Figure 5I). The results of three algorithms were intersected to obtain the candidate biomarkers: CETN3, CACNA1E, OGT, QRFPR, and SCOC (Figure 5J).

Selection and Validation of Two Core Diagnostic Biomarkers

We eventually obtained two core diagnostic biomarkers in AS-FM patients, CETN3 and CACNA1E, by intersecting the machine learning results of two conditions (Figure 6A). Thus, these two core biomarkers were further investigated to expand their clinical significance in the onset of FM among patients suffering from AS.

Figure 6 Selection and validation of the two shared diagnostic genes. (A) The Venn plot exhibited the two shared diagnostic genes. (B) Differential expression of CETN3 in the discovery cohort for AS and FM. (C) ROC curve of CETN3 in the discovery cohort for AS and FM. (D) Differential expression of CACNA1E in the discovery cohort for AS and FM. (E) ROC curve of CACNA1E in the discovery group for AS and FM. (F) Differential expression of CETN3 in the validation cohort for AS and FM (G) ROC curve of CETN3 in the validation cohort for AS and FM. (H) Differential expression of CACNA1E in the validation cohort for AS and FM. (I) ROC curve of CACNA1E in the validation group for AS and FM. * p <0.05; ** p <0.01; *** p < 0.001.

Abbreviations: AS, Ankylosing Spondylitis; FM, Fibromyalgia; ROC, Receiver Operating Characteristic.

First of all, we examined CETN3 and CACNA1E expressions in the discovery datasets of AS and FM. Figure 6B revealed that CETN3 was significantly elevated in the AS (p < 0.01) and FM groups (p < 0.001). CACNA1E expressed higher both in the AS (p < 0.05) and the FM groups (p < 0.001) (Figure 6D). ROC analysis was then applied to inspect the specificity and sensitivity of the two hub genes for the diagnosis of AS and FM. Regarding AS biomarkers, CETN3 and CACNA1E obtained the AUCs of 0.785 and 0.727. The FM group was then also subjected to ROC analysis. Diagnostic value was high for both biomarkers: CETN3 (AUC= 0.770), and CACNA1E (AUC= 0.720) (Figure 6C and E).

The dependability of CETN3 and CACNA1E as core diagnostic genes for AS and FM was further verified in the validation cohorts. In the two validation cohorts, CETN3 expression was remarkedly increased in the AS (P < 0.001) and FM (P < 0.05) groups (Figure 6F). Meanwhile, ROC curves indicated CETN3 had good diagnostic performance in the validation groups of AS (AUC = 0.720) and FM (AUC = 0.750) (Figure 6G). Moreover, CACNA1E was significantly elevated in the AS (P < 0.05) and FM (P < 0.001) groups (Figure 6H). Analogously, CACNA1E also properly diagnosed AS (AUC = 0.751) and FM (AUC = 0.766) (Figure 6I). In summary, the results verified the potential of CETN3 and CACNA1E to serve as key discriminatory genes for AS and FM.

Validation of CETN3 and CACNA1E by RT-PCR in Human Whole Blood Samples

To validate the mRNA abundance of the two diagnostic markers in clinical samples, RT-PCR was then conducted on the blood samples from the healthy controls and AS-FM patients. Compared with the healthy controls, quantification of the mRNA abundance of CETN3 (p < 0.05) and CACNA1E (p < 0.01) were both markedly up-regulated in the AS-FM patients (Figure 7A and B).

Figure 7 Transcriptomic validation of CETN3 (A) and CACNA1E (B) in AS-FM patients (n=9) as well as in healthy controls (n=10). A nonparametric Student’s t-test with Welch’s correction was adopted when comparing the two groups. * p <0.05; ** p <0.01.

Abbreviations: AS-FM, ankylosing spondylitis- fibromyalgia; HC, healthy controls.

Construction of Nomogram and Prediction of Clinical Risk

The nomogram was built in the light of the two hub genes, and a calibration plot was set up to evaluate the predictive value of the model (Figure 8A and 8B). The Calibration plot illustrated the minimal deviation between the actual event risk and the predicted event risk, manifesting the acceptable reliability of the models (Figure 8C and D).

Figure 8 Nomogram construction and diagnostic performance validation. (A) The nomogram was visualized using 2 candidate biomarkers to predict the risk of AS. (B) The calibration plot showing the accuracy of the nomogram in AS. (C) Visualization of FM risk model using nomogram. (D) The calibration plot shows the accuracy of the nomogram in FM.

Abbreviations: AS, Ankylosing Spondylitis; FM, Fibromyalgia.

Immune Infiltration Analysis and Correlation Between Hub Genes and Immune Cells

Immune and inflammatory responses have been demonstrated to play vital roles in the progression of AS and FM, and the ssGSEA algorithm was employed to calculate the constitution of immune cells. The heatmap exhibited the relative content of distinct immune cells in each sample of AS and FM (Figure 9A and B). Further, we compared the proportion of immune cells between AS and controls, as well as FM and their corresponding controls. In comparison with healthy controls, AS patients had a higher content of activated dendritic cells, CD56+ natural killer cells, γΔ.T cells, macrophage, regulatory T cells, T follicular helper cells, and Type 17 T helper cells (Figure 9C). However, a comparison of FM and controls indicated many significantly decreased cell types such as CD56± natural killer cell, γΔ.T cell, MDSC, macrophage, monocyte, etc (Figure 9D). Interestingly, similar differences were discovered in T cell subtypes such as Type 17 T helper cells and Natural killer T cells.

Figure 9 Immune cell infiltration analysis with the ssGSEA algorithm. (A) The heatmap exhibited the difference between healthy control and AS of 28 types of immune cells. (B) The heatmap exhibited the difference between healthy control and FM of 28 types of immune cells. (C) Comparison of 28 types of immune cells between healthy control and AS using the violin plot. (D) Comparison of 28 types of immune cells between healthy control and FM using the violin plot. (E) Correlation analysis between hub genes and immune cells in AS patients. (F) Correlation analysis between hub genes and immune cells in FM patients. * p <0.05; ** p <0.01; *** p < 0.001.

Abbreviations: AS, Ankylosing Spondylitis; FM, Fibromyalgia.

The linkage between the immune cells and hub genes was further investigated. In both AS and FM groups, CETN3 expression was negatively correlated with regulatory T cell and CACNA1E was negatively related with five T cell subsets such as effector memory CD8 T cell, effector memory CD4 T cell, central memory CD8 T cell, central memory CD4 T cell and activated CD8 T cell (Figure 9E and F).

Visualization of Molecular Docking Simulation Toward Hub Genes

To explore the potential therapeutic drugs directed at AS patients with FM, we tried to screen small molecule drugs toward CETN3 and CACNA1E with molecular docking. Unfortunately, the protein structure coded by CETN3 was not found in the RCSB database therefore we only downloaded the crystal structure of CACNA1E with high resolution. Three potential small molecular drugs, gabapentin, gabapentin enacarbil, and pregabalin were filtered out with binding energy were −5 kcal/mol, −6.1 kcal/mol, and −4.6 kcal/mol respectively (Figure 10A–C). The binding poses and sites are illustrated in Figure 10, where the blue dotted lines and number mean the hydrogen-binding and hydrogen bond length, and the green and red represent the small molecular drugs.

Figure 10 Molecular docking patterns of CACNA1 complexed with (A) gabapentin, (B) gabapentin enacarbil, and (C) pregabalin.

Discussion

Due to the high prevalence of fibromyalgia among patients with AS,7,23 it has become an urgent need to probe the shared pathogenesis of the two disorders. In the study, we employed the DEGs analysis, WGCNA, and three machine learning algorithms to excavate the shared genes and two diagnostic genes of CETN3 and CACNA1E which were both markedly associated with AS and FM. Furthermore, we discovered that common pathways via enrichment analysis, and potential drugs were predicted for treating FM in AS patients. To sum up, the newly discovered diagnostic genes and potential molecular mechanisms in the present research offered fresh clinical insights and advice for the diagnosis and therapy of FM and AS patients.

Pain is thought as one of the vital clinical diagnostic criteria for AS and FM, however, the manifestations and mechanisms behind it are distinct. For AS-associated pain, it is intermittent early on and gradually aggravates into persistent.24 It is reported that inflammatory pain and neuropathic pain exist in AS patients.25,26 In an inflammatory state, immune cells gather around nociceptors in nerve endings, particularly in conditions of chronic pain.27 Meanwhile, inflammation biomarkers have been demonstrated to be linked with other diseases such as obstructive sleep apnea syndrome (OSAS), and provided avenues for future study.28 Of note, previous studies suggested HLA-B27 and PTX-3 as effective diagnostic markers between AS and FM patients.10,14 In this study, we identified CACNA1E and CETN3 as reliable blood biomarkers in predicting AS and FM. The pain in patients with FM is characterized by long duration (> 3 months) and generalized.1 Growing evidence suggests that neuroinflammation happening in the peripheral nociceptors, and central nervous system (CNS) is closely associated with the pathophysiology of FM,29–31 and abnormal activity of immune cells such as neutrophils, has been demonstrated to mediate the pain in FM.6,32 Centrin-3, which is encoded by CETN3, is a vital component of Ca2+ binding proteins in the conserved EF-hand family.33 Interestingly, CACNA1E, the other hub gene identified in this study, encodes the α1-subunit of the voltage-gated CaV2.3 channel.34 Of note, calcium channel dysfunction has been reported to result in multisystem disorders, eg arrhythmia and autism,35 and targeting calcium channels has been demonstrated as an ideal strategy for controlling neuropathic pain.36 Thus, we presumed that aberrantly high expression of CETN3 and CACNA1E in AS-FM patients might relate to the activation of calcium channels in CNS and peripheral nociceptors, thereby causing chronic pain in patients.

Up to now, the effective pharmacotherapy for the of AS-FM is still insufficient, under this case, screening the potential drugs is of great significance.

In recent years, impressive achievements have been made in the area of identifying therapeutic small molecular compounds in various disorders. Robust evidences suggest that small molecular compounds on the strengths of high tissue penetration, an adjustable half-life, and ideal oral bioavailability, which enable them efficient in treating.37 Pregabalin, depicted as a structurally semblable of γ- aminobutyric acid (GABA), has been exploited as a promising first-line agent for treating neuropathic pain.38,39 Previous studies have noted that pregabalin could alleviate chronic pain via diminishing glutamatergic activity in the central nervous system and of FM patients.40 However, there was little literature that reported screening of potential small molecular compounds as an option for therapeutic strategies for AS-related FM based on the selection of gene signatures via gene expression profiling in the patient’s whole blood. Herein, via molecular docking, this study offered a fresh insight linking AS-FM pathogenic genes to discover the potential drugs, and we eventually got 3 small molecular compounds (gabapentin, gabapentin enacarbil (GpEn), and pregabalin). Strikingly, GpEn, an extended-release prodrug of gabapentin, showed the lowest binding energy of −6.1 Kcal/mol, implying that it may reduce the pathogenic gene expression in AS-FM patients to the maximum extent. Although no direct evidence is revealed about the effect of GpEn therapy in treating AS or FM, evidence has demonstrated that GpEn can ameliorate restless legs syndrome and postherpetic neuralgia (PHN).41–43 Interestingly, the effect of GpEn has been validated at spine, supraspinal, and peripheral sites and it works via binding calcium channels and also blocks synaptogenesis.43,44 According to the previous findings above, the pharmacological mechanisms of GpEn make it possible to be an ideal drug for treating neuroinflammatory pain, and early intervention with GpEn may improve prognosis and quality of life in AS-FM patients.

FM is widely regarded as a disease caused by a central sensitization phenomenon featured by the disorder of neuro-circuits, which covers the complex process of afferent nociceptive stimuli, with prevailing pain symptoms in the locomotor system.4 In recent years, an increasing amount of proof has started to realize that increased levels of prooxidative factors like mitophagy can lead to pain sensitization in FM.45 Interestingly, oxidative stress-mediated mitochondrial dysfunction was demonstrated to promote mesenchymal stem cell aging in AS,46 which suggests mitochondria-related dysfunction may play a key role in the pathogenesis of AS and FM.47,48 For this research, the results of GO and KEGG enrichment analysis illustrated the shared DEGs for AS and FM were principally engaged in mitophagy and cell cycle-related pathways, manifesting that mitophagy and cell cycle-related pathways might play critical roles in the pathophysiology of AS-related FM.

When we investigated the immune cell infiltration states of AS and FM, distinct modes were uncovered. Multiple immune cells were activated in AS, resulting in various immune responses, consisting of nascent vascularization and osteogenesis.49 However, the research on immune cells in FM is very limited when compared to AS. Interestingly, Yao et al identified a subtype of FM featured by micro-inflammation status and this FM subgroup correlated with depression and activated innate immune response. They further uncovered that the proportion of immune infiltrating cells such as activated dendritic cells (DCs) and natural killer (NK) cells were remarkedly elevated in the micro-inflammation FM group.50 We observed the large difference in immune infiltrating cells between these two diseases, moreover, CACNA1E was positively correlated with macrophages and MDSC in AS but they are negatively linked in FM. However, it is worth noting that CETN3 demonstrated opposite correlations with effector memory CD8 T cells and central memory CD4 T cells in AS and FM, indicating these linkages do not necessarily reflect a causal association between hub genes and immune response between these two diseases. In the context of FM, neuroinflammation is regarded as a possible pathological mechanism with the activation of the immune system.51,52 As a result, further research is urgently required to explore the underlying function of immune infiltrating cells and correlations of two hub genes in the pathogenesis of AS and FM.

Despite this, this study is subjected to some constraints. Firstly, although we identified CETN3 and CACNA1E as hub genes in AS-FM patients using bioinformatics various methods, the mechanistic experiments including in vivo and in vitro are still required to be further investigated. Secondly, this study analyzed the gene expression data focusing only on whole blood, a more comprehensive study should include transcriptome sequencing datasets targeted tissues and potential heterogeneity in patient populations should be taken into account. Thirdly, although gabapentin enacarbil is identified as a potential therapeutic agent, clinical trials should be implemented to substantiate the curative effect.

To sum up, the study aimed at whole blood gene profiling data of patients with AS and FM, combined DEG analysis, WGCNA, and machine learning method, eventually determined CETN3 and CACNA1E as shared diagnostic genes of AS and FM. Plus, we validated the expression levels of hub genes in AS-FM patients. Moreover, ssGSEA unraveled the correlation between two shared genes and immune cells. Eventually, molecular docking analysis revealed that GpEn has the potential to become an effective therapeutic drug for AS-FM patients. This study offered potential diagnostic biomarkers and a fresh view for elucidating the pathophysiology behind AS-associated FM.

Ethics Approval

Experiments involving human samples in this study were authorized and approved by the Shenzhen Nanshan People’s Hospital. All participants were informed of their involvement in this study and signed informed consent.

Consent for Publication

All authors read the final version of this manuscript and endorsed it for publication.

Data Sharing Statement

The datasets involved in this study are available from the GEO (Gene Expression Omnibus) database (https://www.ncbi.nlm.nih.gov/geo/) under accession number GSE25101, GSE73754, GSE67311, and GSE221921.

Funding

This study was financed by the Medical Research Foundation of Guangdong Province (No. B2024010) and Shenzhen Nanshan District Science and Technology Plan Project (NS2024038).

Disclosure

The authors declare this research was implemented without any conflicts.

References

1. Bair MJ, Krebs EE. Fibromyalgia. Ann Internal Med. 2020;172(5):Itc33–itc48. doi:10.7326/AITC202003030

2. Sarzi-Puttini P, Giorgi V, Marotto D, Atzeni F. Fibromyalgia: an update on clinical characteristics, aetiopathogenesis and treatment. Nat Rev Rheumatol. 2020;16(11):645–660. doi:10.1038/s41584-020-00506-w

3. Queiroz LP. Worldwide epidemiology of fibromyalgia. Curr Pain Headache Rep. 2013;17(8):356. doi:10.1007/s11916-013-0356-5

4. Siracusa R, Paola RD, Cuzzocrea S. Impellizzeri D: fibromyalgia: pathogenesis, Mechanisms, Diagnosis and Treatment Options Update. Int J Mol Sci. 2021;22(8):1.

5. Ranganathan V, Gracey E, Brown MA, Inman RD, Haroon N. Pathogenesis of ankylosing spondylitis - recent advances and future directions. Nat Rev Rheumatol. 2017;13(6):359–367. doi:10.1038/nrrheum.2017.56

6. Gau SY, Lee YH, Tsou HK, et al. Patients With Ankylosing Spondylitis Are Associated With High Risk of Fibromyalgia: a Nationwide Population-Based Cohort Study. Front Med. 2021;8:618594. doi:10.3389/fmed.2021.618594

7. Almodóvar R, Carmona L, Zarco P, et al. Fibromyalgia in patients with ankylosing spondylitis: prevalence and utility of the measures of activity, function and radiological damage. Clin Exp Rheumatol. 2010;28(Suppl 63):S33–39.

8. Salaffi F, De Angelis R, Carotti M, Gutierrez M, Sarzi-Puttini P, Atzeni F. Fibromyalgia in patients with axial spondyloarthritis: epidemiological profile and effect on measures of disease activity. Rheumatol Int. 2014;34(8):1103–1110. doi:10.1007/s00296-014-2955-9

9. Sieper J, Poddubnyy D. Axial spondyloarthritis. Lancet. 2017;390(10089):73–84. doi:10.1016/S0140-6736(16)31591-4

10. Brewerton DA, Hart FD, Nicholls A, Caffrey M, James DC, Sturrock RD. Ankylosing spondylitis and HL-A 27. Lancet. 1973;1(7809):904–907. doi:10.1016/S0140-6736(73)91360-3

11. Zhang M, Zhou J, Wang H, et al. Exploration of the shared pathways and common biomarker PAN3 in ankylosing spondylitis and ulcerative colitis using integrated bioinformatics analysis. Front Immunol. 2023;14:1089622. doi:10.3389/fimmu.2023.1089622

12. Jacques P, Elewaut D, Mielants H. Interactions between gut inflammation and arthritis/spondylitis. Curr Opin Rheumatol. 2010;22(4):368–374. doi:10.1097/BOR.0b013e3283393807

13. Hackshaw KV. The Search for Biomarkers in Fibromyalgia. Diagnostics. 2021;11(2):156.

14. Baraka E, Balata M, Ahmed S, El-Blbehisy M, Elattar E. Concomitant Diagnosis of Fibromyalgia and Ankylosing Spondylitis: relation to Clinical Features and Plasma Pentraxin −3 Level. Curr Rheumatol Rev. 2021;17(3):331–341. doi:10.2174/1573397117666210114110823

15. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinf. 2008;9(1):559. doi:10.1186/1471-2105-9-559

16. Motamedi F, Pérez-Sánchez H, Mehridehnavi A, Fassihi A, Ghasemi F. Accelerating Big Data Analysis through LASSO-Random Forest Algorithm in QSAR Studies. Bioinformatics. 2022;38(2):469–475. doi:10.1093/bioinformatics/btab659

17. Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinf. 2013;14(1):7. doi:10.1186/1471-2105-14-7

18. Yoo M, Shin J, Kim J, et al. DSigDB: drug signatures database for gene set analysis. Bioinformatics. 2015;31(18):3069–3071. doi:10.1093/bioinformatics/btv313

19. Kuleshov MV, Jones MR, Rouillard AD, et al. Enrichr: a comprehensive gene set enrichment analysis web server 2016 update. Nuc Acids Res. 2016;44(W1):W90–97. doi:10.1093/nar/gkw377

20. Sieper J, Rudwaleit M, Baraliakos X, et al. The Assessment of SpondyloArthritis international Society (ASAS) handbook: a guide to assess spondyloarthritis. Ann Rheumatic Dis. 2009:68. Suppl 2:ii1–44.

21. Pimentel-Santos FM, Ligeiro D, Matos M, et al. Whole blood transcriptional profiling in ankylosing spondylitis identifies novel candidate genes that might contribute to the inflammatory and tissue-destructive disease aspects. Arthritis Res Ther. 2011;13(2):R57. doi:10.1186/ar3309

22. Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2−ΔΔCT method. Methods. 2001;25(4):402–408. doi:10.1006/meth.2001.1262

23. Ablin JN, Eshed I, Berman M, et al. Prevalence of Axial Spondyloarthritis Among Patients With Fibromyalgia: a Magnetic Resonance Imaging Study With Application of the Assessment of SpondyloArthritis International Society Classification Criteria. Arthr Care Res. 2017;69(5):724–729. doi:10.1002/acr.22967

24. Bidad K, Gracey E, Hemington KS, Mapplebeck JCS, Davis KD, Inman RD. Pain in ankylosing spondylitis: a neuro-immune collaboration. Nat Rev Rheumatol. 2017;13(7):410–420. doi:10.1038/nrrheum.2017.92

25. Wu Q, Inman RD, Davis KD. Neuropathic pain in ankylosing spondylitis: a psychophysics and brain imaging study. Arth Rheumat. 2013;65(6):1494–1503. doi:10.1002/art.37920

26. Gok K, Cengiz G, Erol K, Ozgocmen S. Neuropathic Pain Component in Axial Spondyloarthritis and the Influence on Disease Burden. J Clin Rheumatol. 2018;24(6):324–327. doi:10.1097/RHU.0000000000000711

27. Kim CF, Moalem-Taylor G. Interleukin-17 contributes to neuroinflammation and neuropathic pain following peripheral nerve injury in mice. J Pain. 2011;12(3):370–383. doi:10.1016/j.jpain.2010.08.003

28. Lavalle S, Masiello E, Iannella G, et al. Unraveling the Complexities of Oxidative Stress and Inflammation Biomarkers in Obstructive Sleep Apnea Syndrome: a Comprehensive Review. Life. 2024;14(4):425. doi:10.3390/life14040425

29. Slade GD, Conrad MS, Diatchenko L, et al. Cytokine biomarkers and chronic pain: association of genes, transcription, and circulating proteins with temporomandibular disorders and widespread palpation tenderness. Pain. 2011;152(12):2802–2812. doi:10.1016/j.pain.2011.09.005

30. Sturgill J, McGee E, Menzies V. Unique cytokine signature in the plasma of patients with fibromyalgia. J Immunol Res. 2014;2014:938576. doi:10.1155/2014/938576

31. Mendieta D, De la Cruz-Aguilera DL, Barrera-Villalpando MI, et al. IL-8 and IL-6 primarily mediate the inflammatory response in fibromyalgia patients. J Neuroimmunol. 2016;290:22–25. doi:10.1016/j.jneuroim.2015.11.011

32. Caxaria S, Bharde S, Fuller AM, et al. Neutrophils infiltrate sensory ganglia and mediate chronic widespread pain in fibromyalgia. Proc Natl Acad Sci USA. 2023;120(17):e2211631120. doi:10.1073/pnas.2211631120

33. Ying G, Frederick JM, Baehr W. Deletion of both centrin 2 (CETN2) and CETN3 destabilizes the distal connecting cilium of mouse photoreceptors. J Biol Chem. 2019;294(11):3957–3973. doi:10.1074/jbc.RA118.006371

34. Carvill GL. Calcium Channel Dysfunction in Epilepsy: gain of CACNA1E. Epilepsy Curr. 2019;19(3):199–201. doi:10.1177/1535759719845324

35. Splawski I, Timothy KW, Sharpe LM, et al. Ca(V)1.2 calcium channel dysfunction causes a multisystem disorder including arrhythmia and autism. Cell. 2004;119(1):19–31. doi:10.1016/j.cell.2004.09.011

36. Patel R, Montagut-Bordas C, Dickenson AH. Calcium channel modulation as a target in chronic pain control. Br J Pharmacol. 2018;175(12):2173–2184. doi:10.1111/bph.13789

37. Zhang B, Dömling A. Small molecule modulators of IL-17A/IL-17RA: a patent review (2013-2021). Expert Opin Ther Patent. 2022;32(11):1161–1173. doi:10.1080/13543776.2022.2143264

38. Taylor CP, Angelotti T, Fauman E. Pharmacology and mechanism of action of pregabalin: the calcium channel alpha2-delta (alpha2-delta) subunit as a target for antiepileptic drug discovery. Epilepsy Res. 2007;73(2):137–150. doi:10.1016/j.eplepsyres.2006.09.008

39. Zhang Z, Jiang J, He Y, et al. Pregabalin mitigates microglial activation and neuronal injury by inhibiting HMGB1 signaling pathway in radiation-induced brain injury. J Neuroinflam. 2022;19(1):231. doi:10.1186/s12974-022-02596-7

40. Harris RE, Napadow V, Huggins JP, et al. Pregabalin rectifies aberrant brain chemistry, connectivity, and functional response in chronic pain patients. Anesthesiology. 2013;119(6):1453–1464. doi:10.1097/ALN.0000000000000017

41. Imamura S, Kushida C. Gabapentin enacarbil (XP13512/GSK1838262) as an alternative treatment to dopaminergic agents for restless legs syndrome. Exp Opin Pharmacother. 2010;11(11):1925–1932. doi:10.1517/14656566.2010.494598

42. Winkelman JW, Bogan RK, Schmidt MH, Hudson JD, DeRossett SE, Hill-Zabala CE. Randomized polysomnography study of gabapentin enacarbil in subjects with restless legs syndrome. Movement Disorder. 2011;26(11):2065–2072. doi:10.1002/mds.23771

43. Thomas BM, Farquhar-Smith P. Gabapentin enacarbil extended release for the treatment of postherpetic neuralgia in adults. Ther Clin Risk Manage. 2013;9:469–475. doi:10.2147/TCRM.S50212

44. Bennett MI, Simpson KH. Gabapentin in the treatment of neuropathic pain. Palliative Med. 2004;18(1):5–11. doi:10.1191/0269216304pm845ra

45. Cordero MD, De Miguel M, Moreno Fernández AM, et al. Mitochondrial dysfunction and mitophagy activation in blood mononuclear cells of fibromyalgia patients: implications in the pathogenesis of the disease. Arthritis Res Ther. 2010;12(1):R17. doi:10.1186/ar2918

46. Ye G, Xie Z, Zeng H, et al. Oxidative stress-mediated mitochondrial dysfunction facilitates mesenchymal stem cell senescence in ankylosing spondylitis. Cell Death Dis. 2020;11(9):775. doi:10.1038/s41419-020-02993-x

47. Ramírez-Tejero JA, Durán-González E, Martínez-Lara A, et al. Microbiota and Mitochondrial Sex-Dependent Imbalance in Fibromyalgia: a Pilot Descriptive Study. Neurol Int. 2023;15(3):868–880. doi:10.3390/neurolint15030055

48. Liu Z, Gao L, Wang P, et al. TNF-α Induced the Enhanced Apoptosis of Mesenchymal Stem Cells in Ankylosing Spondylitis by Overexpressing TRAIL-R2. Stem Cells Int. 2017;2017:4521324. doi:10.1155/2017/4521324

49. Cai C, Huang Y, Li L, et al. Angiogenesis-related immune response may be the prelude to the syndesmophyte formation in Ankylosing spondylitis. Int Immunopharmacol. 2024;133:112040. doi:10.1016/j.intimp.2024.112040

50. Yao M, Wang S, Han Y, et al. Micro-inflammation related gene signatures are associated with clinical features and immune status of fibromyalgia. J Transl Med. 2023;21(1):594. doi:10.1186/s12967-023-04477-w

51. Giorgi V, Bazzichi L, Batticciotto A, et al. Fibromyalgia: one year in review 2023. Clin Exp Rheumatol. 2023;41(6):1205–1213. doi:10.55563/clinexprheumatol/257e99

52. Littlejohn G. Neurogenic neuroinflammation in fibromyalgia and complex regional pain syndrome. Nat Rev Rheumatol. 2015;11(11):639–648. doi:10.1038/nrrheum.2015.100