The Role and Potential Mechanism of Complement Factor D in Fibromyalgia Development

Introduction

Fibromyalgia (FM) is a syndrome characterized by chronic musculoskeletal pain.^1–3 The diagnostic criteria for FM are still being updated over time, and the current criteria issued by the American College of Rheumatology (ACR) include: chronic (>3 months) widespread or multisite pain, fatigue, sleep disturbances, cognitive problems, other somatic symptoms (headaches, abdominal pain or bloating, dizziness, paresthesia), diffuse and significant soft tissue tenderness on examination.⁴ The cause of FM is still unknown. Researches to date has proven that genetic factors, female sex, and pre-existing pain of the body are the risk factors of FM.^5,6 The principles of treating FM focus on symptomatic treatment and suggest a combination of individualized and multimodal treatment approach, lacking specific therapeutic targets. Study shows some disturbances in pain processing and regulation in patients with FM, with altered levels of pain-related neuropeptides and neurotransmitters.^1,7 Changes in secretory factor levels in FM patients make the exploration of targeted therapeutics meaningful. The current treatment of FM is mostly symptomatic, so the development of new targeted therapies is important to reduce pain and improve the quality of life of patients with FM.

Inflammation of the nervous system is thought to be a potential cause of FM, and pro-inflammatory cytokines such as IL-6, IL-8 are increased in FM patients relative to healthy controls.⁸ Immune cells can act as mediators in the development of inflammation,⁹ for example macrophages in muscle have been shown to contribute to the development of chronic muscle pain.¹⁰ Pro-inflammatory cytokines can activate and sensitize injury receptors, induce pain and trigger nociceptive sensitization.¹¹ Anti-nuclear antibody (ANA) positive FM patients have a more pronounced inflammatory response than ANA negative controls, suggesting that autoimmunity may be involved in the inflammatory pain of FM.¹² Therefore, it is believed that there is an association between FM and neuroinflammation. Besides, some FM patients respond to treatment with N-methyl-D-aspartate (NMDA) glutamate receptor antagonists,¹³ and a low glutamate diet has been shown to reduce FM symptoms.¹⁴ Elevated glutamate levels may be involved in the development of FM. However, FM patients do not always show good tolerance to treatment with glutamate receptor antagonist drugs.¹⁵ Currently, therapeutic drugs for FM include tricyclic antidepressants, serotonin–norepinephrine reuptake inhibitors, gabapentinoids, and analgesics.^1,2 The current clinical treatment for FM is numerous and heterogeneous, with very limited efficacy for some patients, and lack of specific drugs to relieve symptoms. Since neuroinflammation is inextricably linked to FM, it is clinically relevant to explore specific therapeutic factors that inhibit the inflammatory response in FM.

In this study, we screened differentially expressed genes (DEGs) by analyzing and examining RNA-seq data in the fibromyalgia and healthy groups from an open-source data, establishing functional enrichment analysis and protein interaction networks, and further screening for 5 signature genes of FM, which were then validated by cellular experiments. This study explores the specific changes of FM at the molecular level and provides a reliable theoretical basis for the targeted treatment of FM.

Materials and Methods

Data Sources

The GSE67311 dataset was obtained from the GEO database (https://www.ncbi.nlm.nih.gov/geo/) and contained 142 samples, including 67 FM group and 75 control group whole blood gene expression data. All datasets were subjected to background elimination to remove missing, low-expressed, and non-corresponding gene probes,¹⁶ resulting in a total of 23,306 genes annotated.

Cell Culture

Mouse microglial cell line BV2 and human renal epithelial cell line 293T were purchased from the China Typical Culture Collection (Wuhan, China). The cells were inoculated in DMEM medium containing 10% fetal bovine serum and 1% penicillin – streptomycin and cultured in a constant temperature incubator at 37 °C with 5% CO₂. When the cell fusion reached 70–80%, the cells were passaged, and cells in logarithmic growth phase were taken for cell experiments.

Main Instruments and Reagents

DMEM culture medium, PBS solution, 0.25% Trypsin-EDTA trypsin were purchased from HyClone Laboratories (Logan, USA). Adenosine triphosphate (ATP) was obtained from Sigma-Aldrich (St. Louis, USA). Human recombinant CFD protein was obtained from Solarbio Science & Technology (Beijing, China). miRNA-mimic, primer sequences for IL4 and IL6 were brought from Sangon Biotech (Shanghai, China). Extraction reagents, reverse transcription kits and fluorescent quantitative PCR reagents were purchased from Takara Bio (Kyoto, Japan).

Differentially Expressed Genes

In this study, the RNA-seq data of GSE67311 (fibromyalgia group vs healthy group) were subjected to t-test and Bayesian test using the limma package in R language. Differentially expressed genes (DEGs) were screened for |Log2(FC)| >0.25 and P < 0.05, with negative numbers representing down-regulation and positive numbers representing up-regulation. Pearson correlation coefficients were calculated to determine the gene-to-gene correlations, with |R|>0.2 and P<0.05 as filtering correlation conditions.

Functional Enrichment Analysis of Differentially Expressed Genes and Protein Interaction Network (PPI)

The software packages “org.Hs.eg.db”, “clusterProfiler”¹⁷ and “enrichplot”¹⁸ of R were applied Gene ontology functional analysis (GO) and Kyoto Encyclopedia of Genes and Genomes pathway enrichment analysis (KEGG) analysis were performed and the results were visualized. Prediction of differentially expressed protein interactions was performed using STRING online database, setting the effective binding fraction >0.7.

Molecular Typing and Weighted Co-Expression Network Analysis of Fibromyalgia Samples

The NMF molecular typing model was constructed using the “Consensus ClusterPlus” analysis package in R language. The NMF hierarchical clustering was performed using the adjusted and unified dataset with the clusters k-value ranging from 2 to 8, and the k-value with better clustering stability were selected for the molecular typing according to the clustering effect.¹⁹ Similarly, the “WGCNA” toolkit in R language was used to construct a scale-free weighted gene co-expression network to identify co-expressed genes and modules associated with clinical traits.²⁰

Characteristic Gene Screening

Characterization pre-selection of FM differentially expressed genes by LASSO regression using R package glmnet.¹⁷ LASSO combines feature selection and model construction by adding penalty constraints to the algorithm. In this study, all FM differential gene expression data were applied as feature selection inputs to the LASSO logistic regression model. We used gradient feature selection to evaluate the expression values of each gene. Features were selected when a t-test result of P < 0.05. The regularization parameter λ is selected by 10 times cross-validation so that the error of the selected model is within one standard deviation from the minimum error.²¹ The obtained differentially expressed characteristic genes were crossed with the secreted protein library (https://www.proteinatlas.org/humanproteome/proteinclasses) and the resulting ones were the characteristic genes of secreted protein.

Determination of BV2 Cell Viability by CCK-8

Logarithmically grown BV2 cells were inoculated in 6-well plates and divided into blank group, ATP model control group and treatment group (ATP+human recombinant protein CFD), with three parallel sub-wells of 1.5×10³ cells per well in each group, and the cells were placed in a 37°C, 5% CO₂ incubator overnight. On the second day, 1 µL/mL DMSO solution was added to the control group, 300 µmol/L ATP to the ATP model control group and 300 µmol/L ATP and 1µL/mL rCFD to the treatment group, respectively. After incubation for 48 h, the supernatant was aspirated and 100 µL of culture medium was added to each well, followed by 10 µL of CCK-8 solution. The cells were then incubated in a cell incubator protected from light. After incubation for 2 h in the incubator, the absorbance was read at 450 nm.

Detection of IL4 and IL6 Gene Expression by Quantitative Real-Time Fluorescence Polymerase Chain Reaction (qRT-PCR)

BV2 cells were inoculated in 6-well plates with 1×10⁵ cells per well, and the cells were incubated at 37 °C in a 5% CO₂ incubator. After cell fusion reached 60–70%, the medium was changed and DMSO solution (1 µL/mL), ATP solution (300 µmol/L) and ATP solution (300 µmol/L) + rCFD (1 µL/mL) were added to each group, respectively. Incubation was continued for 24–48 hours. PCR operations were then performed.

Total cellular RNA was extracted by TRIzol method, reverse transcribed into cDNA and then subjected to fluorescent quantitative PCR. Reaction system: SYBR 5 µL, PCR upstream and downstream primers 0.8 µL each, DNA template 1 µL, ddH₂O 3.2 µL, total volume 10 µL. Reaction conditions: preheated denaturation at 95 °C for 30s; cyclic system 95 °C for 15s, 58 °C for 1 min, 40 cycles in total. The relative expression levels of mRNA were normalized by GAPDH and the relative expression was calculated by the 2^−ΔΔCt method.²² The primer sequences are shown in Table 1.

Table 1 List of Primers for qRT-PCR

Glucose, ATP and Lactate Content Testing

Measure key products of glycolysis, including glucose, ATP and lactate, to understand changes in glycolysis. Cell precipitates (1×10⁷ cells or more) were collected for subsequent sonication and used for glucose, ATP and lactate micro-titration assays according to the kit procedure (Solarbio Science & Technology, Beijing, China).

Dual Luciferase Reporter Gene Assay

Firstly, the 3′UTR sequence of CFD gene was obtained from the gene library to construct the wild-type fragment. The wild-type fragment included the binding site of miRNA-3193. The purified PCR product was cloned into the psiCHECK-2 vector named psiCHECK-2-CFD-3′UTR. Subsequently, the firefly luciferin signal and sea kidney luciferin signal were measured in the cells using a fluorescence luminescence detector according to the instructions of the dual luciferase reporter gene assay kit (Hanheng Bio, Shanghai, China).

Statistical Analysis

GraphPad Prism 8 (GraphPad Software, San Diego, USA) software was used for data analysis. The mean data conforming to normal distribution were expressed as mean ± standard deviation (x ± s). The comparison of the means between 2 groups was performed by independent sample t-test, and the comparison between multiple groups was performed by ANOVA analysis. The difference was considered statistically significant at P<0.05.

Results

Screening for FM-Related Differentially Expressed Genes

By preprocessing the GEO-GSE67311 transcriptional profile and normalizing the data, a total of 23,306 mRNAs were identified. 70 differentially expressed genes were screened in the FM group compared with the healthy population (|Log2(FC)| >0.25 and P < 0.05), of which 38 genes were significantly upregulated and 32 genes were downregulated (Figure 1A). Next, GO function annotation and KEGG pathway enrichment analysis of 70 differentially expressed genes revealed that the GO functions of differentially expressed genes were mainly enriched in antioxidant response, tyrosine phosphorylation, myofibrillogenesis and contraction, actin assembly, endopeptidase activity, growth factor activity regulation and receptor ligand activation (Figure 1B–D). KEGG was mainly enriched in IL-17 signaling pathway, tumor necrosis factor signaling pathway, cytokine receptor interaction, PI3K-Akt signaling pathway and complement system (Figure 1E). Then, to further clarify the role of these genes in fibromyalgia, this study used the STRING website to construct a protein-protein interaction network for these 62 differentially expressed genes. The results showed that IFIT1, RSAD2, EPB42, IFI44, IFI44L and SLC4A1 were the key node genes (Figure 2A and B).

Figure 1 Identification and enrichment analysis of FM differentially expressed genes. (A) GEO-GSE67311Heat map of metabolically differentially expressed genes in the dataset. (B) Biological process (BP). (C) Cellular component (CC). (D) Molecular function (MF). (E) KEGG.

Figure 2 Molecular typing and weighted co-expression network analysis of FM samples. (A) PPI network. (B) Number of core genetic link nodes. (C) The consistency matrix for K=2,3,4. (D) Cumulative distribution function when K=2-9. (E) Plot of correlation coefficient as a function of k value. (F) Correlation diagram between modules and traits. (G and H) Soft threshold determination.

Molecular Typing and Weighted Co-Expression Network Analysis of FM-Related Differentially Expressed Genes

Traditionally, the diagnosis of disease relies on pathology,²³ but this hardly reflects the biological nature of the disease. New approaches to disease typing based on molecular typing can provide global characterization of disease at the genetic level, greatly deepen the understanding of molecular pathological information, thus this typing approach has an important clinical role. To verify whether the above 70 FM-associated differentially expressed genes have clinical translational potential, we used the GEO-Merge dataset for molecular typing studies of FM.

Based on the enrichment scores of 70 genes, we performed a clustering analysis of 68 FM specimens using K-Means consistency clustering. However, stable clustering results could not be obtained when k=2-9 (Figure 2C–E), ie, the 70 differentially expressed genes could not be used for typing of dilated cardiomyopathy samples, and their mechanisms of action need further investigation.

Next, we performed WGCNA construction of the GEO-GSE67311 dataset. The weighted co-expression network was constructed using the WGCNA package in the R toolkit, however, our further study revealed that the gene co-expression modules (Figure 2F–H) were not obtained by WGCNA analysis, which also suggests that there may be more complex mechanisms of action for these genes that need further investigation.

FM-Related Trait Gene Screening

The 70 FM-related differentially expressed genes were characterized by the LASSO algorithm,¹⁷ and 20 factors were screened as candidate feature genes (Figure 3A). Subsequently, we intersected these genes with the secreted proteins in the secreted protein library and finally obtained five FM-related differentially secreted genes, namely CPA3, IL3RA, CFD, OCLN and IGKV6-21 (Figure 3B). We found that secreted factors CPA3, IL3RA, CFD, OCLN and IGKV6-21 were all down expressed in FM (Figure 3D). Further ROC curves also suggested AUC values of 0.745, 0.670, 0.640, 0.578, and 0.648 for CPA3, IL3RA, CFD, OCLN, and IGKV6-21, respectively (Figure 3C). Based on this, we screened 5 genes of secreted proteins by bioinformatics analysis, which suggested that these genes may have a role in inhibiting the process of FM.

Figure 3 Screening for FM-related trait genes. (A) The LASSO algorithm screened 20 glycolysis differentially expressed signature genes. (B) Characterized genes are intersected with secreted proteins in the secreted protein library. (C) Analysis of diagnostic efficiency (AUC values) of CPA3, IL3RA, CFD, OCLN and IGKV6-21 genes. (D) Differential expression analysis of CPA3, IL3RA, CFD, OCLN and IGKV6-21 in FM. ns means no significant difference, **P<0.01, ***P<0.001, compared with the control group.

Construction of ceRNA Network Upstream of 5 Secretory Factors

In this study, the miRNAs and lncRNAs upstream of CPA3, IL3RA, CFD, OCLN and IGKV6-21 were analyzed by TargetScan, miRanda and miRDB databases. A total of 36 CPA3-associated miRNAs, 2 IL3RA-associated miRNAs, 1 CFD-associated miRNA, 31 OCLN-associated miRNAs and 0 IGKV6-21-associated miRNAs were intersected by the three databases (Figure 4). Next, lncRNA prediction of miRNAs upstream of CPA3, IL3RA, CFD and OCLN was performed using the spongeScan database, and a total of 36 lncRNAs were obtained. We performed ceRNAs network construction of these predicted miRNAs and lncRNAs using Cytoscape software. The results showed that hsa-miR-590-3p and hsa-miR-377-3p were the core node genes of OCLN (Figure 4).

Figure 4 Construction of ceRNA network upstream of 5 secretory factors. Construction of ceRNA networks for CPA3, IL3RA, CFD and OCLN genes.

Correlation Analysis of 5 Characteristic Secretory Genes with IL-4, IL-6, IL-10 and TNF Genes

The correlations of CPA3, IL3RA, CFD, OCLN and IGKV6-21 with IL4, IL6, IL10 and TNF genes were analyzed using the GEO-GSE67311 dataset (Figure 5A–D). We found that complement factor D (CFD) was negatively correlated with the pro-inflammatory factor IL-6 (Figure 5A) and positively correlated with the anti-inflammatory factor IL-4 (Figure 5C). Therefore, we suggested that CFD may play an important role in the pain-causing of FM.

Figure 5 Correlation analysis of 5 secreted protein genes with IL-4, IL-6, IL-10 and TNF genes. (A) Correlation analysis of CPA3, IL3RA, CFD, OCLN and IGKV6-21 with IL-6 genes. (B) Correlation analysis of CPA3, IL3RA, CFD, OCLN and IGKV6-21 with TNF genes. (C) Correlation analysis of CPA3, IL3RA, CFD, OCLN and IGKV6-21 with IL-4 genes. (D) Correlation analysis of CPA3, IL3RA, CFD, OCLN and IGKV6-21 with IL-10 genes. Red-coloured text means the difference is significant (P<0.05).

Effect of Recombinant Protein CFD (rCFD) on the Viability of ATP-Induced BV2 Cells

Based on the findings, we considered that the secreted factor CDF may have an important role in neuropathic pain of FM. First, the expression of CFD in healthy individuals and FM groups was compared using the GEO-GSE67311 dataset. The results showed that the CFD was significantly decreased in the FM group (Figure 3C), suggesting that CFD may be a pain suppressor in FM. Then, we applied 1ug/mL and 2ug/mL rCFD to BV2 cells and calculated the cell viability by CCK8 assay. The results showed that the concentration of rCFD less than 2ug/mL was not toxic to BV2 cells (Figure 6A). Third, ATP-induced BV2 cells were applied to establish a neuropathic pain model,²⁴ and the cell survival rate was examined after the addition of rCFD. We found that 1ug/mL rCFD had an anti-ATP-induced change in neuronal cell viability compared to the control (Figure 6B).

Figure 6 Effect of recombinant protein CFD (rCFD) on ATP-induced BV2 cell viability. (A) Toxic effects of rCFD on BV2 cells. (B) Effect of rCFD on ATP-induced changes in BV2 cell viability. (C) 1ug/mL rCFD inhibited the expression of the pro-inflammatory factor IL-6. (D) 1ug/mL rCFD promoted the expression of the anti-inflammatory factor IL-4. (E and F) CFD with miRNA-3193 luciferase assay. *** P<0.001, compared with the negative control group. Gold highlighting means the base pairing region.

Next, we studied the effect of secreted factor rCFD on the expression of IL-4 and IL-6 genes in BV2 cells induced by ATP. The results showed that ATP could significantly increase the expression of IL-6 in BV2 cells; the addition of 1ug/mL rCFD significantly inhibited the expression of IL-6 compared with the ATP group (Figure 6C). Similarly, ATP significantly decreased the expression of IL-4 in BV2 cells; while the addition of 1ug/mL rCFD significantly promoted the expression of IL-4 compared with the ATP group (Figure 6D). These findings suggested that secreted factor CFD can affect neuropathic pain and inflammatory responses of FM.

Subsequently, to investigate the upstream regulatory mechanism of CFD, combining the results of TargetScan, miRanda and miRDB database analysis, we verified the targeting relationship between CFD and miRNA-3193. The results showed that CFD was a direct target gene of miRNA-3193 (Figure 6E and F).

Effect of rCFD on Glucose Metabolism in BV2 Cells

There are few studies on glucose metabolism and FM, and it is unknown whether fibromyalgia can regulate glucose metabolism in glial cells through secreted factors and thus affect inflammatory response and pain. To clarify the pathogenic mechanism of FM, we extracted genes related to glucose metabolism from the MSigDB website (Systematic name: M1870) combined with the GEO-GSE67311 dataset for differential analysis. The results showed that 26 genes showed differential expression in FM compared to the controls (Figure 7A). Next, we analyzed the correlation of CFD with these 26 genes. Our study revealed a negative correlation of CFD with MIOX, CHPF and ESRRB (Figure 7B), and a positive correlation with GOT1, RAE1, ZNF292, ZNF292, NSDHL, ECD, EIF6, VCAN, PDHA1, FAM162A, CHST12, PMM2 and CXCR4 (Figure 7B). In conclusion, our findings suggested that CFD may have a function in regulating glucose metabolism in FM.

Figure 7 Effect of recombinant protein CFD (rCFD) on glucose metabolism in BV2 cells. (A) Identification and analysis of 325 differentially expressed genes for gluconeogenesis in the GEO-GSE67311 dataset. (B) Correlation of CFD with 26 differentially expressed genes. (C–E) Effect of rCFD on changes in glucose, ATP and lactate expression in glucose metabolism of BV2 cells. * P<0.05, ** P<0.01, *** P<0.001, compared with the control group. Red-coloured text means the difference is significant (P<0.05).

To verify this result, the levels of glucose, ATP and lactate in BV2 cells were explored in this study. The levels of glucose (Figure 7C), ATP (Figure 7D) and lactate (Figure 7E) were significantly decreased in BV2 cells after the addition of rCFD compared to the control. These results suggested that CFD may affect neuropathic pain and inflammatory responses through glucose metabolism.

Discussion

Fibromyalgia can have a significant impact on patients both physically and psychologically.²⁵ At present, FM lacks specific treatment and is often treated symptomatically, requiring consideration of both somatic and psychological factors, which may be related to the fact that the etiology and pathogenesis of fibromyalgia are still unclear.²⁶ The search for molecular markers specific to FM has facilitated the development of its therapeutic modalities.²⁷

Bioinformatics combined with big biological data analysis has allowed researchers to gain a deeper understanding of disease at the genetic level, especially in cancer, where a range of biomarkers have been identified to predict cancer recurrence, metastasis, chemosensitivity, and prognosis,²⁸ and these technological advances have provided unprecedented opportunities to explore the pathogenesis of FM. We performed differential expression gene analysis using the GSE67311 dataset derived from the GEO database (67 FM cases versus 75 control cases) and screened 70 differentially expressed genes, which provides a basis for further focusing and screening of new molecular targets.

Secreted proteins are mainly a group of proteins with important physiological functions such as cytokines, growth factors, complement, degrading enzymes, antibodies, peptide hormones and immunoglobulins. More attention is currently being paid to secreted proteins in the field of oncology, as they have become a major source of potential tumor markers that can be used as targets to guide tumor diagnosis, treatment and prognosis.^29,30 Studies on the correlation between FM and secreted proteins have been mainly limited to cytokines,³¹ IL-4, IL-6, IL-10 and TNF factors are involved in the polarization of microglia, where IL-6 and TNF factors have the ability to promote the inflammatory response, while IL-4 and IL-10 have the ability to inhibit the inflammatory response.² However, there are relatively few studies to further investigate whether secreted proteins can be used as potential markers or therapeutic targets for FM. In this study, we screened FM differentially characterized genes by LASSO and obtained a total of 5 FM differentially expressed secreted proteins. We analyzed the correlation between 5 secreted factors and IL-4, IL-6, IL-10 and TNF factors. Since CFD was negatively correlated with the pro-inflammatory factor IL-6 and positively correlated with the anti-inflammatory factor IL-4, this suggested that CFD may be an anti-inflammatory factor and may have a suppressive effect on the development of FM.

CFD, also known as adipsin, is mainly synthesized by adipose tissue and sciatic nerve tissue. Current studies report that CFD plays an important regulatory role in the physiological and pathological processes of systemic sclerosis, diabetes (diabetic nephropathy), obesity, dyslipidemia, and hypertension.^32,33 Recent studies have also suggested the possible involvement of CFD in respiratory disorders,³⁴ neurological disorders³⁵ and even preterm abortion.³⁶ However, the role of CFD in FM has not been studied. As mentioned above, we observed that CFD is involved in the regulation of the levels of inflammation-related factors IL-4 and IL-6. Besides, pro-inflammatory cytokines, such as IL-1β, IL-6, and TNF-α, can activate and sensitize injury receptors, and another potential source of them is adipose tissue. There is also evidence that diffuse or multifocal pain is more common in obese individuals.³⁷ Weight loss can improve symptoms of FM,³⁸ and weight loss of at least 10% is associated with greater improvement in FM symptoms and an increase in IL-10.³⁹ Besides, researches have indicated that FM patients showed reduced responses to adrenocorticotropic hormone, epinephrine, cortisol and growth hormone compared to healthy controls, suggesting insulin resistance,⁴⁰ and that FM can be prevented by improving the control of blood glucose levels.⁴¹ Based on the above findings, we suggest that CFD is likely to be involved in regulating the occurrence of FM through the regulation of energy metabolism. Through cellular assays, we verified that CFD reversed the decrease in IL-4 expression in BV2 cells caused by ATP; similarly, CFD reversed the increase in IL-6 expression in BV2 cells caused by ATP. Therefore, it is reasonable to believe that CFD can be used as a therapeutic target for FM or as an inhibitor of inflammatory factors in combination with conventional FM drugs for the treatment of FM. At the last, though we have specified that CFD was negatively correlated with the pro-inflammatory factor IL-6 and positively correlated with the anti-inflammatory factor IL-4, further experiment needs to be performed to elucidate the detailed mechanism and signal pathway underlying the effect of CFD on anti-inflammation.

Conclusion

In summary, based on bioinformatic methods and big data mining we obtained a new target CFD for FM, and further experiments verified that CFD has significant inhibition of ATP-induced neuropathic pain. These findings provide a new theoretical basis for the clinical translation of CFD.