Manual Therapy Improves Fibromyalgia Symptoms by Downregulating SIK1

Fibromyalgia (FM), classified by ICD-11 with code MG30.0, is a chronic debilitating disease characterized by widespread pain, fatigue, cognitive impairment, sleep, and intestinal alterations, among others. FM affects a large proportion of the worldwide population, with increased prevalence among women. The lack of understanding of its etiology and pathophysiology hampers the development of effective treatments. Our group had developed a manual therapy (MT) pressure-controlled custom manual protocol on FM showing hyperalgesia/allodynia, fatigue, and patient’s quality of life benefits in a cohort of 38 FM cases (NCT04174300). With the aim of understanding the therapeutic molecular mechanisms triggered by MT, this study interrogated Peripheral Blood Mononuclear Cell (PBMC) transcriptomes from FM participants in this clinical trial using whole RNA sequencing (RNAseq) and reverse transcription followed by quantitative Polymerase Chain Reaction (RT-qPCR) technologies. The results show that the salt-induced kinase SIK1 gene was consistently downregulated by MT in FM, correlating with improvement of patient symptoms. In addition, this study compared the findings in a non-FM control cohort subjected to the same MT protocol, evidencing that those changes in SIK1 expression with MT only occurred in individuals with FM. This positions SIK1 as a potential biomarker to monitor response to MT and as a therapeutic target of FM, which will be further explored by continuation studies.


Introduction
The latest version of the International Classification of Diseases (ICD-11) adopted by the WHO (World Health Organization) in May 2019 [1] classifies fibromyalgia (FM) as a multifactorial chronic primary widespread pain syndrome (code MG30.0) presenting diffuse pain in at least four of the five body regions, anxiety, depression, and overall functional disability [1].
Diagnosis is based on clinical criteria defined by the ACR (American College of Rheumatology) 1990 case definition with revisions [2,3].Appropriate diagnosis should ensure that pain is not directly attributable to a nociceptive process but consistent with nociplastic pain [1], which is caused by poorly understood mechanisms, rather than local nerve damage (neuropathic pain) [4].Additional symptoms include non-restorative sleep, fatigue, cognitive impairment, and intestinal problems, overlapping with symptoms present in Myalgic Encephalomyelitis/Chronic Fatigue Syndrome (ME/CFS) [5,6].Presentation of FM peaks between 20 and 55 years with marked an increased prevalence in women [7,8].Epidemiology reports vary across countries and regions, with worldwide impact showing is pre-pandemic (NCT04174300 completed before 03/2020).The average age for the FM cohort was 55.6 ± 7.2 years (range 43-71), and the time from primary FM diagnosis over 3 years was 10.3 ± 7.5 years (range [3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21].A subcohort of six participants composed entirely of women, with an average age of 54 years ± 8.44 and a range of between 43 and 69 years, was selected for PBMC RNAseq analysis to evaluate the effects of the MT program on the immune system of FM (see Section 2.2 for details).

Participant Phenotyping
As mentioned in the study design section above, all participants were finely phenotyped with the use of the FIQ [19,20], MFI [21], and SF-36 (Likert scale) [22] validated questionnaires.As shown in Table 1, participating patients presented severe FM (total FIQ > 59) and moderate fatigue, with most domains showing scores well above the non-FM control cohort, except for general fatigue, where scores were very similar.Quality of life (SF-36) was much superior (>50 in all domains) in the non-FM cohort than the FM cohort, except for the role emotional domain.No major baseline (pretreatment) differences were found between the subcohort of FM (n = 6) subjected to RNAseq analysis and the previously described complete FM cohort (n = 38) [18], while the non-FM cohort statistically differed from the complete, as well as the sequenced subcohorts of FM patients.Individual participant scores are available in Supplementary Table S1.
In addition, a comparison of baseline PPTs showed that low cervical points were the most sensitive to pressure-induced pain, while gluteus, trochanters, and knees were the least sensitive to pain in both FM cohorts, with a tendency to increased sensitivity in the trapezius right point for the RNAseq FM subcohort (n = 6) (p-value = 0.052, Table 2).As expected, non-FM participants showed marked differences in resistance to pressuretriggered pain (Table 2).Individual participant PPT values are available in Supplementary Table S2.

Differential Gene Expression in PBMCs of FM with Therapy
Differential gene expression in PBMCs of FM with MT was assessed by RNAseq analysis of total RNA prepared from an RNAseq subcohort of FM patients (discovery phase) before and after the MT treatment (n = 12 paired samples from six patients).The results showed an overexpression of 72 transcripts and an underexpression of 256 transcripts (p < 0.05, FDR < 0.1) (Figure 1, and Supplementary Table S3).
At the individual level, significant changes with treatment varied across participants, with the upregulation of at least 22 genes and the downregulation of at least 42 in all participants (see Supplementary Figure S1 for individual volcano plots and Supplementary Table S4 for RNAseq DE analysis at the individual level).

Gene Enrichment and Pathway Analysis with MT in the Immune System of FM
To understand the biological significance of the genes differentially expressed (DE) with MT in PBMCs of FM, we performed enrichment analysis using the gene ontology (GO) knowledgebase (http://geneontology.org, accessed on 31 July 2024) [24,25] (Figure 2, left panel) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (https: //www.genome.jp/kegg/pathway.html,accessed on 31 July 2024) (Figure 2, right panel), as detailed in Methods.[19,20], MFI [21], and SF-36 (Likert scale) [22] questionnaire scores by cohort, as indicated.

Differential Gene Expression in PBMCs of FM with Therapy
Differential gene expression in PBMCs of FM with MT was assessed by RNAseq analysis of total RNA prepared from an RNAseq subcohort of FM patients (discovery phase) before and after the MT treatment (n = 12 paired samples from six patients).The results showed an overexpression of 72 transcripts and an underexpression of 256 transcripts (p < 0.05, FDR < 0.1) (Figure 1, and Supplementary Table S3).At the individual level, significant changes with treatment varied across participants, with the upregulation of at least 22 genes and the downregulation of at least 42 in all participants (see Supplementary Figure S1 for individual volcano plots and Supplementary Table S4 for RNAseq DE analysis at the individual level).

Gene Enrichment and Pathway Analysis with MT in the Immune System of FM
To understand the biological significance of the genes differentially expressed (DE) with MT in PBMCs of FM, we performed enrichment analysis using the gene ontology (GO) knowledgebase (http://geneontology.org, accessed on 31/07/2024) [24,25] (Figure 2, left panel) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (https://www.genome.jp/kegg/pathway.html,accessed on 31/07/2024) (Figure 2, right panel), as detailed in Methods.
The findings included responses to stress and infectious processes, with the involvement of cytokine, chemokine, MAPK, and NFkB signaling processes (Figure 2 and Supplementary Table S5).

RT-qPCR Validation of Protein-Coding Genes Differentially Expressed in Response to MT in FM
To validate top relevant functions among the 328 differentially expressed (DE) transcripts with treatment (p < 0.05, FDR < 0.1) that are now in the complete FM cohort of (n = 38), we selected only RNAs with protein-coding potential (Supplementary Table S3, coding potential tab) by removing 149 TCONS (Transcript Cluster Noncoding RNAs), 11 pseudogenes, one miRNA cluster, two lncRNAs (RP4-737E23.5 and RP11-213H15.1),three lincRNAs (LINC00861, LINC00877, and AP001046.5),and seven novel transcripts.Therefore, RNAs with coding potential DE by MT included 161 DE RNAs (22 upregulated and 139 downregulated).Then, we evaluated DE at the individual level and imposed the condition of having DE in at least 50% of the samples (n ≥ 3), (Table 3), which reduced the list to only six genes being consistently overexpressed and eighteen underexpressed by MT.They were then considered top candidate effectors of the therapy.Individual DE data are provided in Supplementary Table S4 and summarized in Table 3 (Indiv.pval< 0.05, FDR < 0.1).The findings included responses to stress and infectious processes, with the involvement of cytokine, chemokine, MAPK, and NFkB signaling processes (Figure 2 and Supplementary Table S5).

RT-qPCR Validation of Protein-Coding Genes Differentially Expressed in Response to MT in FM
To validate top relevant functions among the 328 differentially expressed (DE) transcripts with treatment (p < 0.05, FDR < 0.1) that are now in the complete FM cohort of (n = 38), we selected only RNAs with protein-coding potential (Supplementary Table S3, coding potential tab) by removing 149 TCONS (Transcript Cluster Noncoding RNAs), 11 pseudogenes, one miRNA cluster, two lncRNAs (RP4-737E23.5 and RP11-213H15.1),three lincRNAs (LINC00861, LINC00877, and AP001046.5),and seven novel transcripts.Therefore, RNAs with coding potential DE by MT included 161 DE RNAs (22 upregulated and 139 downregulated).Then, we evaluated DE at the individual level and imposed the condition of having DE in at least 50% of the samples (n ≥ 3), (Table 3), which reduced the list to only six genes being consistently overexpressed and eighteen underexpressed by MT.They were then considered top candidate effectors of the therapy.Individual DE data are provided in Supplementary Table S4 and summarized in Table 3 (Indiv.pval< 0.05, FDR < 0.1).GO and KEGG analyses were then reassessed for these top 24 DE RNAseq genes of individual responses to MT in at least 50% of the FM cases, finding that responses to bacteria and chemokine signaling were among top cell functions affected by MT (Figure 3).
To orthogonally validate DE RNAseq genes by the alternative RT-qPCR method, primer sets were designed with even more restrictive selection criteria, as the amount of material for this step was limited.To select the top six DE genes out of the twenty-four listed in Table 3, we chose the least individually altered with MT in at least four out of the six FM cases from the RNAseq subcohort.With this, we designed primers for CD3E and CX3CR1 from the upregulated group.However, for the downregulated group, the number of genes fulfilling the applied criteria (being DE in at least 4 out of the 6 patients) exceeded four genes (Table 3, downregulated), among those the top 4 candidates were SIK1, HBEGF, and EGR2, all individually downregulated in five of the six FM cases.In addition, we also included the top downregulated gene (EREG) in this RT-qPCR validation step (Table 3).Individual details for DE of these genes are provided in Supplementary Table S6.Primer sets to measure the expression levels for these six DE-selected genes (Supplementary Table S7), together with a set to detect the housekeeping GAPDH gene, were then used for technical validation of the RNAseq results in our FM subcohort (n = 6), as well as in the complete FM cohort (n = 38) (population validation of our RNAseq data).
GO and KEGG analyses were then reassessed for these top 24 DE RNAseq genes of individual responses to MT in at least 50% of the FM cases, finding that responses to bacteria and chemokine signaling were among top cell functions affected by MT (Figure 3).To orthogonally validate DE RNAseq genes by the alternative RT-qPCR method, primer sets were designed with even more restrictive selection criteria, as the amount of material for this step was limited.To select the top six DE genes out of the twenty-four listed in Table 3, we chose the least individually altered with MT in at least four out of the six FM cases from the RNAseq subcohort.With this, we designed primers for CD3E and CX3CR1 from the upregulated group.However, for the downregulated group, the number of genes fulfilling the applied criteria (being DE in at least 4 out of the 6 patients) exceeded four genes (Table 3, downregulated), among those the top 4 candidates were SIK1, HBEGF, and EGR2, all individually downregulated in five of the six FM cases.In addition, we also included the top downregulated gene (EREG) in this RT-qPCR validation step (Table 3).Individual details for DE of these genes are provided in Supplementary Table S6.Primer sets to measure the expression levels for these six DE-selected genes (Supplementary Table S7), together with a set to detect the housekeeping GAPDH gene, were then used for technical validation of the RNAseq results in our FM subcohort (n = 6), as well as in the complete FM cohort (n = 38) (population validation of our RNAseq data).S4) and the stringent selection criteria described above, while the lower panel in Figure 3 shows the results obtained by the alternative RT-qPCR approach, performed to validate the RNAseq results in the FM subcohort (n = 6).As observed, only SIK1 and HBEGF downregulation could be validated.None of the four remaining selected DE genes appeared significantly changed with MT by RT-qPCR (Figure 3, lower panel).
However, RT-qPCR analysis of the complete FM cohort (n = 38) confirmed the upregulation of CX3CR1 and the downregulation of all four selected genes among those downregulated by RNAseq.No significant change with MT was found for CD3E by RT-qPCR (Figure 5, upper panel).
To investigate whether DE with MT were exclusive to FM or their change with MT corresponded to a generalized response, DE was also measured in PBMCs from non-FM matched volunteers subjected to the same eight-session MT treatment program previously described [18] (n = 12).As can be observed in the lower panel of Figure 5, while the upregulation of CX3CR1 seemed also upregulated and HBEGF appeared downregulated by MT (interpreted as potential general responders to MT), all the three other genes downregulated in FM by MT (EREG, EGR2, and SIK1) did not show significant changes in non-FM individuals, interestingly, suggesting that these genes may constitute sensors of the response to MT therapy in FM.It should be noted that basal expression levels vary greatly across the two groups being compared, with an approximately 10-fold difference for CX3CR1 and HBEGF, which increased in the FM group and decreased in the indicated group for EREG.The potential significance of these large basal differences is unknown.However, SIK1 basal levels appear similar between the FM and the non-FM cohorts but are only significantly inhibited by MT in the FM group.
Furthermore, since our previous work, in the context of clinical trial NCT04174300, identified differences in response to MT among FM patients co-diagnosed with Myalgic Encephalomyelitis/Chronic Fatigue Syndrome (ME/CFS) [18], we reassessed our RT-qPCR results, taking into account whether or not patients with FM had also received the ME/CFS diagnostic (Supplementary Table S1).The results indeed point out that the FM group with ME/CFS co-diagnosis (n = 19) does not seem to respond to MT by increasing their CX3CR1 levels, while DE of HBEGF and EGR2 appears more related to this group (Figure 6).EREG and SIK1 DE seem to specifically associate with both patient groups, without changes in the control non-FM group (Figure 5).
The upper panel of Figure 4 illustrates individual relative DE values of the proteincoding genes selected for RT-qPCR validation, according to RNAseq individual data (Supplementary Table S4) and the stringent selection criteria described above, while the lower panel in Figure 3 shows the results obtained by the alternative RT-qPCR approach, performed to validate the RNAseq results in the FM subcohort (n = 6).As observed, only SIK1 and HBEGF downregulation could be validated.None of the four remaining selected DE genes appeared significantly changed with MT by RT-qPCR (Figure 3, lower panel).However, RT-qPCR analysis of the complete FM cohort (n = 38) confirmed the upregulation of CX3CR1 and the downregulation of all four selected genes among those downregulated by RNAseq.No significant change with MT was found for CD3E by RT-qPCR (Figure 5, upper panel).
To investigate whether DE with MT were exclusive to FM or their change with MT corresponded to a generalized response, DE was also measured in PBMCs from non-FM matched volunteers subjected to the same eight-session MT treatment program previously described [18] (n = 12).As can be observed in the lower panel of Figure 5, while the upregulation of CX3CR1 seemed also upregulated and HBEGF appeared downregulated by MT (interpreted as potential general responders to MT), all the three other genes downregulated in FM by MT (EREG, EGR2, and SIK1) did not show significant changes in non-FM individuals, interestingly, suggesting that these genes may constitute sensors of the response to MT therapy in FM.It should be noted that basal expression levels vary greatly across the two groups being compared, with an approximately 10-fold difference for CX3CR1 and HBEGF, which increased in the FM group and decreased in the indicated group for EREG.The potential significance of these large basal differences is unknown.However, SIK1 basal levels appear similar between the FM and the non-FM cohorts but are only significantly inhibited by MT in the FM group.Furthermore, since our previous work, in the context of clinical trial NCT04174300, identified differences in response to MT among FM patients co-diagnosed with Myalgic Encephalomyelitis/Chronic Fatigue Syndrome (ME/CFS) [18], we reassessed our RT-qPCR results, taking into account whether or not patients with FM had also received the ME/CFS diagnostic (Supplementary Table S1).The results indeed point out that the FM group with ME/CFS co-diagnosis (n = 19) does not seem to respond to MT by increasing

Correlation of Genes Differentially Expressed in Response to MT with Patient Symptoms and Sensitivity to Pain (PPTs)
MT seems to provide improvement of patient symptoms to a certain extent, as shown by the changes detected for overall and total FIQ scores, as well as for the SF-36 "Bodily pain" domain for the complete FM group (n = 38), as previously reported [18].By contrast, changes were not appreciated in the FM RNAseq RNAseq cohort (n = 6), perhaps due to the small sample size, and were absent in the non-FM cohort (n = 12), as expected (Table 4).In the latter, scores seem to rather be associated with improved fatigue and mental health, perhaps related to the relief of some ailments unrelated to FM.
On another side, the monitoring of PPT changes with MT also had shown significant changes with treatment for the most sensitive "Low cervical" tender points (n = 38) [18].Again, no differences were detected for the FM subcohort (n = 6), and unexpectedly, thresholds changed with treatments for the non-FM cohort in different anatomic locations (Table 5).

Correlation of Genes Differentially Expressed in Response to MT with Patient Symptoms and Sensitivity to Pain (PPTs)
MT seems to provide improvement of patient symptoms to a certain extent, as shown by the changes detected for overall and total FIQ scores, as well as for the SF-36 "Bodily pain" domain for the complete FM group (n = 38), as previously reported [18].By contrast, changes were not appreciated in the FM RNAseq RNAseq cohort (n = 6), perhaps due to the small sample size, and were absent in the non-FM cohort (n = 12), as expected (Table 4).In the latter, scores seem to rather be associated with improved fatigue and mental health, perhaps related to the relief of some ailments unrelated to FM.
On another side, the monitoring of PPT changes with MT also had shown significant changes with treatment for the most sensitive "Low cervical" tender points (n = 38) [18].Again, no differences were detected for the FM subcohort (n = 6), and unexpectedly, thresholds changed with treatments for the non-FM cohort in different anatomic locations (Table 5).
To determine if the DE genes with MT play a role in improvement of symptoms, we evaluated the potential correlations between gene expression and symptom differences across our validated results.Figure 7, interestingly, shows that low levels of SIK1 negatively correlate with higher scores for the SF-36 "Bodily pain" domain, which indicates better health (reduced pain), supporting a potential role of the MT treatment in reducing pain in FM by decreasing SIK1 levels.Positive correlations of SIK1 levels with FIQ "Overall" and "Symptoms" domains also support their potential participation in treatment-associated patient improvement.By contrast, high MFI "Physical fatigue" domain scores associated with lower SIK1 levels would indicate worsening of fatigue in FM with treatment (Figure 7, upper left panel).Importantly, no relevant correlations of SIK1 levels with questionnaire score changes after treatment were detected in the control group (Figure 7, upper right panel), supporting a specific role of SIK1 in FM symptom relief in FM with MT.Although the correlations with total FIQ and SF-36 vitality with the upregulation of CX3CR1 and those of HBGEF with the same domains correlating with SIK1, except for the SF-36, seem to mainly associate with symptoms improvement with a change in expression, the observed changes in these two genes are not specific to FM (Figure 5).
On another end, lower levels of SIK1 seem to correlate with increased threshold values in the lower cervical left tender point of FM patients (n = 38) (Figure 7, lower left panel), indicating that the changes in SIK1 may be associated with improvement of patient allodynia.No correlations of the changes in SIK1 with PPT value improvement were detected in the control non-FM group (n = 12) (Figure 7, lower right panel).Statistically significant correlations of PPT changes in other DE genes (CX3CR1 with low cervical right and EREG with gluteal left) were considered spurious since DE of CX3CR1 and EREG were found to be non-significant by RT-qPCR in this control group (Figure 5), and, therefore, their detailed analysis was not further pursued.
As differences in response to MT in patients having received ME/CFS co-diagnosis had been previously reported [18], we set to evaluate potential correlations between DE genes and symptoms or between DE genes and PPT changes with ME/CFS co-diagnosis.As shown in Figure 8 (upper left), patients not fulfilling diagnosis criteria for ME/CFS show a clear benefit of the MT program applied, mostly reproducing the results obtained with the full cohort.By contrast, the subgroup of FM participants that had received a diagnosis of ME/CFS as well seemed to present a more reduced benefit of symptoms associated with decreased levels of SIK1 (Figure 8, upper right), with significant improvement of the symptom domain of the FIQ questionnaire for low SIK1 levels only.Again, significant associations between symptoms and other DE genes, different from SIK1, were no further pursued, as DE of the gene had not passed the test of significance (e.g., CX3CR1 in the ME/CFS co-diagnosed group or EGR2 in the FM subgroup not having received co-diagnosis of FM) (see Figure 6) or had shown lower correlation values (e.g., CX3CR1 with mental health of the SF-36 questionnaire) (Figure 8, upper left).
With respect to changes in PPT values with reduced SIK1 levels, once more, the FM subgroup not fulfilling ME/CFS diagnosis criteria seems to behave as the complete FM group, indicating that 50% of patients in the group (FM with ME/CFS diagnosis) were mostly irresponsive to MT, at least symptom wise with SIK1 changes.Statistically significant associations unrelated to SIK1 were not further examined, as the affected PPTs (e.g., great trochanters or lateral epicondyle humerus) were not among the affected by MT in FM (Table 5).

Discussion
This study expands our previous knowledge on the improvement of FM symptoms by a self-designed controlled-pressure MT protocol (NCT04174300) [18] by evidencing molecular changes in the immune system of FM participants with MT.This study comprised two phases: a discovery phase of genome-wide transcriptomic profiling to detect changes in expression levels with MT in the immune system of an RNAseq FM subcohort (n = 6) and a validation phase, extending the main molecular findings to the complete cohort (n = 38).This later validation phase also examined changes in expression levels with MT in the immune system of a non-FM control cohort (n = 12) treated with the same self-designed controlled-pressure MT protocol as the FM cases [18].The objective was to find out if the observed findings in the immune system with MT were specific to FM or, by contrast, corresponded with a general mechanism triggered by MT in all individuals.Although limitations associated with the selection process of an RNAseq subcohort of FM and the selection of DE genes with MT leave room for further findings, the results strikingly show that the downregulation of SIK1 correlates with patient symptom improvement, particularly with some FIQ domains ("Symptoms" and "Overall"), as well as with the SF-36 "Bodily pain" subdomain, with the latter two domains having shown most improvement with MT [18].They also show that this correlation is specific for FM, as SIK1 levels do not seem to change with MT in the immune system of control non-FM participants.Whether and how SIK1 changes affect FM immune transcriptome [16,17] will require further work.However, our finding that basal levels between the FM and non-FM cohorts appear similar in our RT-qPCR data and the absence of a direct connection of SIK1 with FM phenotype, to the best of our knowledge, may indicate an indirect effect of SIK1 inhibition to mediate MT therapeutic action.SIK1, initially identified for its role in sodium sensing, belongs to the salt-inducible kinases (SIKs) family, which includes three homologous serine-threonine kinases (SIK1, SIK2, and SIK3) that regulate multiple aspects of the human physiology in response to extracellular signals, including feeding/fasting metabolic responses, inflammation and immune responses, and sleep (circadian rhythms), among others [26].Of the three kinases, only SIK1 expression is upregulated at the transcriptional level through a consensus CREB (cAMP response element) present in its promoter, as shown in myocytes and the suprachiasmatic nucleus of the brain (SCN) [26][27][28].SIK activity regulates innate immunity responses by suppressing the production of the IL-10 anti-inflammatory cytokine in macrophages.In fact, pharmacological inhibition of SIK activity increases the levels of IL-10 while suppressing the levels of the proinflammatory IL-6, IL-12, and TNF-α after TLR (toll-like receptor) stimulation by LPS [27,29].However, conflictive data regarding the production of proinflammatory cytokines and activation of the transcription factor NFκB with increased SIK activity exist [30,31].
The current intense research in the development of member-specific inhibitors of SIK activity [32][33][34] should eventually help to ascertain the precise attributes and contributions of each member of this family of proteins, in particular cell and environmental scenarios, leading to the development of novel pharmacological treatments.For example, to increase the production of IL-10 in the gut, Sundberg et al. screened a library of kinase inhibitors after challenging murine bone-marrow-derived dendritic cells (DCs) with the yeast cell wall preparation zymosan, finding that the protective effects involved SIK activity inhibition in a subpopulation of CD11c (+) CX3CR1(hi) cells isolated from murine gut tissue [32].Thus, SIK activity seems relevant in still other immune system compartments, including mast cell IL-33 cytokine release [35], and it even modulates adaptive immunity through the regulation of T-cell lineage commitment, differentiation, and survival [36,37], therefore offering a multitude of potential indirect mechanisms to exert MT therapeutic effects in FM.Although drastic SIK1 downregulation may not be desirable because of its role in blood pressure or its tunning in certain cell types or diseases could constitute valued therapeutic options [30,38,39].
Whether SIK1 transcription downregulation by MT is mediated through the conserved CREB element in its promoter or through alternative mechanisms seems an important question for future work in the field of physiotherapy.Other possibilities worth exploring after this initial finding are the potential impact of MT on the muscle, blood pressure, and cell metabolism or on the circadian system through changes in SIK activity.
By contrast, CX3CR1 levels appeared significantly increased in both study groups (FM and non-FM individuals), indicating that MT triggers this change in all individuals, with the exception of those FM patients co-diagnosed with ME/CFS.CX3CR1 is a Gprotein-coupled receptor and the only binder of fractalkine present on a subset of immune cells, including monocytes and macrophages, as well as DCs, T helper (Th) 1, CD8+T effector/memory and γδ T lymphocytes, and NK cells [40].Its main role in immune cells is to detect and migrate toward inflamed tissue, "crawling and patrolling" from blood vessel endothelium to different destinies according to fractalkine's gradient, the objective being to initiate innate immune responses followed by adaptive responses [40,41].In the brain, it is mainly expressed in astrocytes and microglia regulating cellular communication between neurons, in addition to providing protection from the neurotoxicity induced by the HIV-1 envelope protein gp120 [42].In the gut, CX3CR1-positive macrophages produce the IL-10 immunoregulatory cytokine [43], and a lack of CX3CR1 expression is associated with an altered microbiome and impaired intestinal barrier [44].Regulatory mechanisms of CX3CR1 expression and the implications of its overexpression are complex and require further research to understand their impact on health and disease.Why patients codiagnosed with ME/CFS do not respond to MT with increased CX3CR1 levels is unknown at present, but it seems to support differential response to MT with ME/CFS co-diagnosis, as previously shown [18].Molecular basal differences between patients fulfilling only FM or ME/CFS (codiagnosis) have been found by our group [45] and by others [46], seemingly demanding a review of the case definition for patients fulfilling both clinical criteria [45].Our previous report, NCT04174300 [18] showed differences in response to MT between patients that had or had not received a co-diagnosis of ME/CFS.The results of this study further confirm differences across these two FM subgroups, not only for a lack of upregulation of CX3CR1 levels in response to MT but also for the downregulation of EGR2, occurring only in the co-diagnosed group.EGR2 or early growth response 2 is a transcription factor with an essential epigenetic regulatory role (DNA methylation turnover) for the differentiation of human monocytes [47] and a novel regulator of the senescence of fibroblast and epithelial cells [48].Together with EGR3, it is needed for T-and B-cell development and activation [49].Whether MT preferential upregulation of EGR2 in patients with an ME/CFS status relates to increased EGR2 basal levels in these patients (as shown in Figure 6), coinciding with Dr. Kerr's previous findings [50,51], and whether this relates to EBV infection history of the patient, seems like a possibility to be further explored.
Finally, the downregulation of EREG, also known as epiregulin, seems to discriminate responses to MT between both FM subgroups and the control.Being a soluble peptide hormone involved in inflammation and wound healing, upregulated by LPS induction and by stress of the endoplasmic reticulum [52], its downregulation by MT may relate to patient improvement.However, correlations with questionnaire scores did not detect such a link.
The fact that RT-qPCR did not validate the increased levels of CD3E with MT detected by RNAseq does not serve to refute its findings, as the methodological differences may indeed constitute the reason for the discrepancy.
On the question of what could the mechanisms that exert changes in the molecular profiles of immune cells by MT be, we are far from being able to give a detailed response.However, elucidation of the immunomodulatory effects of massage either by direct pressure/mechanotransduction or by indirect pathways effected by MT, such as cytokine, chemokine, microRNA release [12], sleep improvement [53], or others, are on their way.

Study Design and Intervention
This is an observational study consisting of the analysis of the molecular changes taking place in the circulating immune cells of FM patients (n = 38) and non-FM patients (n = 12) by a physiotherapy program of manual therapy.The program consisted of eight sessions (twice weekly for four weeks) with a 25 min custom protocol, including pressure maneuvers of about 4.5 N each of 10 out of the 18 FM tender points and surrounding areas (NCT04174300), as described in our previous publication [18].This study also included a pilot interventional non-randomized single-arm trial, replicating the treatment program previously applied to FM and now using a matched non-FM cohort as a control.Participants could not be enrolled in pharmacological CTs or receive additional physiotherapy treatment while participating in this study and agreed to withdraw medication 12 h before blood draws.Comparison across groups was performed to find out whether the MT program triggers similar or distinct changes in FM vs. non-FM individuals.Before vs. after FIQ [19,20], MFI [21], and SF-36 [22] questionnaire scores, as well as PPT values of the 18-FM tender points, were registered for all participants.The studies were approved by the Universidad Católica de Valencia San Vicente Mártir Ethics Committee with study codes UCV/2018-2019/076 and UCV/2020-2021/167, respectively.All participants signed an informed consent before they were included in this study.For the analysis of molecular changes in the immune system of FM, the whole transcriptome of an RNAseq subcohort of FM patients (n = 6) was obtained before and after the program.Top differentially expressed genes were validated in the complete FM cohort (n = 38) by RT-qPCR and studied in the non-FM cohort (n = 12) for their comparison.

Total RNA Preparation and Quality Assessment
Total RNA was prepared from a previous collection of PBMC pellets (−150 • C) (≥10 6 cells) registered at the Institute of Health Carlos III National Biobank, Madrid, Spain (Ref.C.0006924) [12] with the RNeasy Mini Kit (Qiagen, Venlo, The Netherlands cat.74104) following the manufacturer's protocol.Cell lysis included the addition of 1/100 β-mercaptoethanol (Sigma, Burlington, MA, USA cat.63689) before a 5 min vortex to favor lysis.Removal of contaminant DNA was performed in a column with the RNasefree DNase set (Qiagen, cat.79254).Elution RNase-free ddH 2 O was supplemented with RNasin (Promega, Madison, WI, USA cat.N2118) to a final concentration of 0.4 U/µL before use.RNA yields and quality were obtained with an Agilent TapeStation 2100 (Agilent Technologies, Santa Clara, CA, USA).RNA integrity was double checked by agarose gel electrophoresis and individual electropherograms (Supplementary Figure S2).Only samples with RIN ≥ 7 were subjected to downstream analysis.

RNAseq
Upon ribosomal RNA depletion, 1 µg of RNA was used for whole transcriptome sequencing (Illumina NovaSeq 6000 PE150, 50M reads) (Novogene, Cambridge, UK).After the removal of filtered to adapter sequences and low-quality reads, sequences were aligned to the human GRCh38/hg38 genome using HISAT2 software v.2.2.0 [54].Sequence assembly analysis was performed using Cufflinks, converted to BAM format, and sorted and indexed with samtools [55].Cuffdiff was used to calculate differential expression, expressed as FPKM (fragments per kilobase of exon per million mapped fragments) numbers for each sample (https://cole-trapnell-lab.github.io/cufflinks/cuffdiff/index.html,accessed on 30 March 2020).Differences were given as the absolute value log2 fold change in the ratio between the pretreated and the posttreatment samples > 1, p-adjust < 0.05, and FDR < 0.1 = TRUE.Heatmaps and volcano plots were generated with the CummeRbund R package (R version 4.2.1, cummeRbund 2.38.0) [56] and GraphPad Prism 8.0.2 software.

Enrichment Analysis
Gene ontology (GO) enrichment analysis of differentially expressed genes was performed with the gene ontology (GO) knowledgebase (http://geneontology.org, accessed on 31 July 2024) [24,25] or with the Goseq, version 1.56.0 and R package, version 4.4 [57].After gene length bias correction, p-values less than 0.05 were considered significantly enriched by differentially expressed genes.
For pathway analysis of differentially expressed genes, the Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis of differentially expressed genes was performed with the https://www.genome.jp/kegg,accessed on 31 July 2024, online tool.

RT-qPCR Validation
Reverse transcription was performed using the High Capacity Complementary DNA (cDNA) Reverse Transcription Kit (Applied Biosystems, Waltham, MA, USA, cat.4308228) with 1-2 µg of total RNA, according to the manufacturer's guidelines.Relative gene expression was assessed by qPCR of triplicates per sample, with the primer sets in Supplementary Table S7.qPCR was performed with the PowerUP Sybr Green Master Mix (Applied Biosystems, Foster City, CA, USA cat.100029283) on a Lightcycler LC480 instrument (Roche, Penzberg, Germany) with the following amplification conditions: a single preactivation cycle of the hotstart polymerase at 94 • C for 15 min, followed by 40 amplification cycles, each of which consisted of three steps: 95 • C for 15 s, 54 • C for 30 s, and extension at 70 • C for 30 s. Gene expression levels were normalized to GAPDH and quantified by the 2 −∆∆Ct method [58].

Statistics
Continuous data are expressed as means ± SD (standard deviation) and range values.Statistical differences were assessed by paired t-tests for normal value distributions and either nonparametric Mann-Whitney or Wilcoxon analysis if values did not follow a normal distribution.Normality was determined with the Shapiro-Wilk test.Differences between groups were considered significant when p ≤ 0.05.Analysis was conducted with Excel, the SPSS package 13.0 (SPSS Inc., Chicago, IL, USA), and R v4.2.1 [59].Pearson correlations were evaluated with the WGCNA R package v1.72.5 [60].Plots were drawn using the GraphPad Prism 5.0 program (San Diego, CA, USA) and the package ggplot2 [61].

Conclusions
In conclusion, the molecular data obtained by comparing the immune system transcriptome of FM before and after a pressure-controlled MT protocol (NCT04174300) identifies the downregulation of SIK1 as a prominent and specific target of MT in FM, which is associated with patient symptom improvement.In addition to SIK1's potential biomarker value for monitoring the response to MT in FM, pending external validation in extended cohorts, this pioneer finding opens the interesting possibility of therapeutic applications of SIK1 inhibitors on FM.Future research efforts in the direction of improving FM healthcare seem granted.

Institutional Review Board Statement:
The study was conducted in accordance with the Declaration of Helsinki and was approved by the Institutional Review Board (or Ethics Committee) of UNIVER-SIDAD CATÓLICA DE VALENCIA SAN VICENTE MÁRTIR ETHICS COMMITTEE (protocol code UCV/2018-2019/076, 21 February 2019; and protocol code UCV2020-2021/167, 22 December 2021).All personal data were anonymized in fulfillment of Spanish data protection laws.All tasks were performed by collegiate professionals or qualified trained personnel.
Informed Consent Statement: Informed consent was obtained from all subjects involved in the study.

Figure 1 .
Figure 1.Volcano plot representation of differential gene expression in PBMCs of FM with therapy.|Log2FoldChange| > 1 values (X axis) are displayed with respect to -log10 of their p-values (Y axis), and the significance is set at p < 0.05, FDR < 0.1.

Figure 1 .
Figure 1.Volcano plot representation of differential gene expression in PBMCs of FM with therapy.|Log2FoldChange| > 1 values (X axis) are displayed with respect to log10 of their p-values (Y axis), and the significance is set at p < 0.05, FDR < 0.1.Int.J. Mol.Sci.2024, 25, x FOR PEER REVIEW 12 of 33

Figure 2 .
Figure 2. GO (left) and KEGG (right) pathways targeted by MT in the immune system of FM.Function significance (color palette, padj < 0.05), DE gene count in each pathway are indicated by dot thickness for each panel.

Figure 2 .
Figure 2. GO (left) and KEGG (right) pathways targeted by MT in the immune system of FM.Function significance (color palette, padj < 0.05), DE gene count in each pathway are indicated by dot thickness for each panel.

Figure 3 .
Figure 3. Top pathways targeted by MT in the immune system of FM, as determined by top RNAseq DE genes.Function significance (color palette, padj < 0.05) and approximate DE gene count in each pathway are indicated by dot thickness for each panel.

Figure 3 .
Figure 3. Top pathways targeted by MT in the immune system of FM, as determined by top RNAseq DE genes.Function significance (color palette, padj < 0.05) and approximate DE gene count in each pathway are indicated by dot thickness for each panel.The upper panel of Figure 4 illustrates individual relative DE values of the proteincoding genes selected for RT-qPCR validation, according to RNAseq individual data (Supplementary TableS4) and the stringent selection criteria described above, while the lower panel in Figure3shows the results obtained by the alternative RT-qPCR approach, performed to validate the RNAseq results in the FM subcohort (n = 6).As observed, only SIK1 and HBEGF downregulation could be validated.None of the four remaining selected DE genes appeared significantly changed with MT by RT-qPCR (Figure3, lower panel).However, RT-qPCR analysis of the complete FM cohort (n = 38) confirmed the upregulation of CX3CR1 and the downregulation of all four selected genes among those downregulated by RNAseq.No significant change with MT was found for CD3E by RT-qPCR (Figure5, upper panel).To investigate whether DE with MT were exclusive to FM or their change with MT corresponded to a generalized response, DE was also measured in PBMCs from non-FM matched volunteers subjected to the same eight-session MT treatment program previously described[18] (n = 12).As can be observed in the lower panel of Figure5, while the upregulation of CX3CR1 seemed also upregulated and HBEGF appeared downregulated by MT (interpreted as potential general responders to MT), all the three other genes downregulated in FM by MT (EREG, EGR2, and SIK1) did not show significant changes in non-FM individuals, interestingly, suggesting that these genes may constitute sensors of the response to MT therapy in FM.It should be noted that basal expression levels vary greatly across the two groups being compared, with an approximately 10-fold difference for CX3CR1 and HBEGF, which increased in the FM group and decreased in the indicated group for EREG.The potential significance of these large basal differences is unknown.However, SIK1 basal levels appear similar between the FM and the non-FM cohorts but are only significantly inhibited by MT in the FM group.Furthermore, since our previous work, in the context of clinical trial NCT04174300, identified differences in response to MT among FM patients co-diagnosed with Myalgic Encephalomyelitis/Chronic Fatigue Syndrome (ME/CFS)[18], we reassessed our RT-qPCR results, taking into account whether or not patients with FM had also received the ME/CFS diagnostic (Supplementary TableS1).The results indeed point out that the FM group with ME/CFS co-diagnosis (n = 19) does not seem to respond to MT by increasing their CX3CR1 levels, while DE of HBEGF and EGR2 appears more related to this group (Figure6).EREG

Figure 5 .
Figure 5. DE-expressed genes with MT on PBMCs from FM patients, as determined by RT-qPCR.Relative expression by ΔΔCt values upon GAPDH normalization for each sample (triplicates) in each study group (n = 38 for the FM cohort in blue, upper panel, and n = 12 for the non-FM control cohort in black, lower panel) are shown.A statistical paired two-Wilcoxon test with Benjamin-Hochberg p-value correction.(* p < 0.05, ** p < 0.01, *** p < 0.001, **** p < 0.0001) was applied to assess the significance of DE.

Figure 5 .
Figure 5. DE-expressed genes with MT on PBMCs from FM patients, as determined by RT-qPCR.Relative expression by ∆∆Ct values upon GAPDH normalization for each sample (triplicates) in each study group (n = 38 for the FM cohort in blue, upper panel, and n = 12 for the non-FM control cohort in black, lower panel) are shown.A statistical paired two-Wilcoxon test with Benjamin-Hochberg p-value correction.(ns, non-significant, * p < 0.05, ** p < 0.01, *** p < 0.001, **** p < 0.0001) was applied to assess the significance of DE.

Figure 6 .
Figure 6.DE-expressed genes with MT on PBMCs from FM patients with or without ME/CFS codiagnosis, as determined by RT-qPCR.Relative expression by ΔΔCt values upon GAPDH normalization for each sample (triplicates) in each study group (n = 19 for the FM-only group, dark blue, upper panel; and n = 19 for the FM group with ME/CFS co-diagnosis, light blue, lower panel) are shown.A statistical paired two-Wilcoxon test with Benjamin-Hochberg p-value correction.(* p < 0.05, ** p < 0.01, *** p < 0.001, **** p < 0.0001) was applied to assess the significance of DE.

Figure 6 .
Figure 6.DE-expressed genes with MT on PBMCs from FM patients with or without ME/CFS co-diagnosis, as determined by RT-qPCR.Relative expression by ∆∆Ct values upon GAPDH normalization for each sample (triplicates) in each study group (n = 19 for the FM-only group, dark blue, upper panel; and n = 19 for the FM group with ME/CFS co-diagnosis, light blue, lower panel) are shown.A statistical paired two-Wilcoxon test with Benjamin-Hochberg p-value correction.(ns, nonsignificant, * p < 0.05, ** p < 0.01, *** p < 0.001, **** p < 0.0001) was applied to assess the significance of DE.

Figure 7 .
Figure 7. Symptom improvement with validated DE genes with MT in FM (n = 38) (upper left) and non-FM controls (n = 12) (upper right), and the correlation of PPT ratios (post and pre) with DE genes in FM (n = 38) (lower left) and non-FM controls (n = 12) (lower right).Pearson correlation values and associated p-values (*, p < 0.05; **, p < 0.01; ***, p < 0.001) between gene expression levels and symptom scores or PPT ratios are shown.

Figure 8 .
Figure 8. Symptom improvement with validated DE genes with MT in FM (n = 19) (upper left) and FM with co-diagnosis of ME/CFS (n = 19) (upper right), and the correlation of PPT ratios (post and pre) with DE genes in FM (n = 19) (lower left) and FM with the co-diagnosis of ME/CFS (n = 19) (lower right).Pearson correlation values and associated p-values (*, p < 0.05; **, p < 0.01; ***, p < 0.001) between gene expression levels and symptom scores or PPT ratios are shown.

Figure 8 .
Figure 8. Symptom improvement with validated DE genes with MT in FM (n = 19) (upper left) and FM with co-diagnosis of ME/CFS (n = 19) (upper right), and the correlation of PPT ratios (post and pre) with DE genes in FM (n = 19) (lower left) and FM with the co-diagnosis of ME/CFS (n = 19) (lower right).Pearson correlation values and associated p-values (*, p < 0.05; **, p < 0.01) between gene expression levels and symptom scores or PPT ratios are shown.

Table 3 .
DE -expressed genes with MT on PBMCs from FM patients, as determined by RNAseq.Transcript Post.Value Pre.Value FC (Post/Pre) Log2FC p-Value Indiv.pval< 0.05 Gene_Name UPREGULATED

Table 3 .
DE-expressed genes with MT on PBMCs from FM patients, as determined by RNAseq.
Numbers in the Indv.pvalcolumn indicate the FM participants who showed DE of the indicated genes with MT by RNAseq analysis (p < 0.05; FDR < 0.1), up-or downregulated, as indicated.FC: fold change.Genes selected for RT-qPCR validation appear bolded.

Table 5 .
Patient response to MT as evidenced by PPT differences by the studied cohort.