Gut microbiome and serum metabolome analyses identify molecular biomarkers and altered glutamate metabolism in fibromyalgia.

BACKGROUND: Fibromyalgia is a complex, relatively unknown disease characterised by chronic, widespread musculoskeletal pain. The gut-brain axis connects the gut microbiome with the brain through the enteric nervous system (ENS); its disruption has been associated with psychiatric and gastrointestinal disorders. To gain an insight into the pathogenesis of fibromyalgia and identify diagnostic biomarkers, we combined different omics techniques to analyse microbiome and serum composition.
METHODS: We collected faeces and blood samples to study the microbiome, the serum metabolome and circulating cytokines and miRNAs from a cohort of 105 fibromyalgia patients and 54 age- and environment-matched healthy individuals. We sequenced the V3 and V4 regions of the 16S rDNA gene from faeces samples. UPLC-MS metabolomics and custom multiplex cytokine and miRNA analysis (FirePlex™ technology) were used to examine sera samples. Finally, we combined the different data types to search for potential biomarkers.
RESULTS: We found that the diversity of bacteria is reduced in fibromyalgia patients. The abundance of the Bifidobacterium and Eubacterium genera (bacteria participating in the metabolism of neurotransmitters in the host) in these patients was significantly reduced. The serum metabolome analysis revealed altered levels of glutamate and serine, suggesting changes in neurotransmitter metabolism. The combined serum metabolomics and gut microbiome datasets showed a certain degree of correlation, reflecting the effect of the microbiome on metabolic activity. We also examined the microbiome and serum metabolites, cytokines and miRNAs as potential sources of molecular biomarkers of fibromyalgia.
CONCLUSIONS: Our results show that the microbiome analysis provides more significant biomarkers than the other techniques employed in the work. Gut microbiome analysis combined with serum metabolomics can shed new light onto the pathogenesis of fibromyalgia. We provide a list of bacteria whose abundance changes in this disease and propose several molecules as potential biomarkers that can be used to evaluate the current diagnostic criteria.

Evidence before this study

Fibromyalgia is a complex disease with chronic pain as its primary symptom. To date, no molecular biomarkers exist for it, leaving its diagnosis up to subjective questionnaires. Several alterations in fibromyalgia patients have pointed towards the central nervous system as the origin of this pathology. The gut microbiome can influence the CNS through the gut-brain axis.

Added value of this study

Employing microbiome and metabolomics analysis along with cytokine and miRNA profiling we identified several alterations between healthy controls and fibromyalgia patients that could be used as potential biomarkers. We also studied how the microbiome and metabolomics datasets correlated with each other to elucidate the role of microbiome alterations in host metabolism.

Implications of all available evidence

Taken together, this study provides candidate molecular biomarkers for fibromyalgia, and supports an alteration of neurotransmitter levels in fibromyalgia patients.

Alt-text: Unlabelled Box

Background

Fibromyalgia is a complex disease of unknown pathophysiology, for which no specific molecular biomarkers or biochemical alterations have been identified. In 1990, the American College of Rheumatology (ACR) recognised this syndrome as a disease and proposed the Widespread Pain Index (WPI), determined by measuring tenderness on pressure at 18 defined points, as a major diagnostic indicator. In 2010, the ACR introduced the Severity Score (SS), which also takes into account the associated symptoms and their severity [102]. Thus, the diagnosis of fibromyalgia is currently based on subjective pain evaluation and a set of associated signs and symptoms, which are used to assess the severity of the disease.

Even though the fibromyalgia is a complex disease with a multitude of signs and symptoms associated with many organs, the participation of the Central Nervous System (CNS) in its pathogenesis is broadly acknowledged [33]. Some studies have tried to identify molecular signatures that could explain some of the features of fibromyalgia and have provided some potential biomarkers. Several polymorphisms linked to the metabolism and breakdown of neurotransmitters involved in pain modulation have been identified as specific markers of increased risk of fibromyalgia [2]. Such polymorphisms have been found for the serotonin transporter gene 5-HTT [14,68] and the catechol-O-methyl-transferase (COMT) gene [30,106]. Some environmental factors, such as viral and bacterial infections [10], e.g. HCV infection [9,78] and psychological stressors [32], known to produce alterations in the hypothalamic-pituitary-adrenal (HPA) axis, have been associated with this disease. Fibromyalgia is prevalent in individuals with chronic pain attributable to peripheral pain generators, such as rheumatoid arthritis [1]. At the molecular level, glutamate is elevated in the cerebrospinal fluid of fibromyalgia patients [26,71,85]. A decrease in insular levels of γ-aminobutyric acid (GABA) has also been described [21]. An inflammatory component in the pathogenesis of this disease has also been proposed: certain cells might trigger and perpetuate chronic pain by releasing chemokines and cytokines, such as IL-6 and IL-8, whose levels are elevated in the sera of fibromyalgia patients [62,95].

The microbiome has a significant role in maintaining health [37,47]. Alterations in the gut microbiome have been linked to a large number of diseases, including intestinal bowel disease (IBD) [45] and metabolic [43] and neurological [84,89] disorders [40]. The microbiome has been recurrently associated with the CNS [89], indicating the existence of a gut-brain axis [16,22]. Disturbances in the microbiome might lead, in some cases, to neural disorders such as depression or autism. Some changes linked to microbial gut dysbiosis, understanding dysbiosis as those differences between healthy individuals and disease-specific patients [35], are also associated with symptoms used to determine the SS2 score in the diagnosis of fibromyalgia. The gut-brain axis has been proposed as a bidirectional communication system between the gastrointestinal tract and the brain, involving both neural and humoral mechanisms (reviewed in [15]). The intestinal GABA produced by the bacteria from glutamate might affect the behaviour of the host, and it might be involved in anxiety and depression [8,34,57,88]. Alterations in the microbiome composition can escalate the interactions between bacteria and the gut immune system due to the breakage of the intestinal barrier, promoting the release of pro-inflammatory molecules. Such events have been reported in IBD, where a release of IL-2, IL-17, interferon and/or TNFβ has been observed [41]. Interestingly, several pro-inflammatory cytokines can increase the permeability of the blood-brain barrier [16]. The microbiome also has metabolic, immunological and gut-protecting functions in the host. The fermentation of dietary carbohydrates by gut bacteria, for example, results in the production of short-chain fatty acids (SCFAs). These molecules are essential for the maintenance of the integrity of intestinal barrier [40] and other health-related functions [77], including the correct development and maintenance of the blood-brain barrier [7].

These interactions between the microbiome and other functional systems of the organism has been widely studied. Microbiome data have been scrutinised in conjunction with host's genome, epigenome, transcriptome and metabolome [99]. The integration of different omics data relies mostly on dimension reduction approaches and is not specific to any omics technology, except for the metabolomics data. Correlation, regression and network-based approaches have also been implemented to integrate microbiome data with other omics analyses. As a result, the role of the of the host genome in regulating microbiome composition has been revealed [28]. Combination of Genome Wide Association Studies (GWAS) and microbiome-GWAS has been applied also to assess the impact of diet on microbiome composition. For example, associations between lactase [5] and variations of vitamin D receptor [98] genes with specific bacteria have been reported. Metabolomics-microbiome integration studies using correlation approaches have shown the effect of microbiome on host's insulin sensitivity [70] and on the development and progression of colorectal cancer [66,100]. Metabolomics – microbiome integration studies employing a mix of correlation and network methods have obtained a comprehensive profile of the existent interactions between intestinal mucosa and gut microbiome [58]. The authors of these studies have used standard statistical methods but suggested that new, specific methods are needed for omics integration, to take into account the particular omics data characteristics [99].

The aim of this work was to identify potential molecular biomarkers for fibromyalgia diagnosis and characterisation, employing different omics technologies: the analysis of microbiome from faeces samples and metabolomics, cytokine and miRNA profiling using serum samples.

Methods

Cohort recruitment

Individuals included in the study were recruited in two different hospitals in the Basque Country. Both fibromyalgia patients and healthy individuals were given a form with questions concerning several lifestyle variables (diet, smoking, alcohol consumption, physical exercise, other diseases and mood),. Blood samples were obtained from fibromyalgia patients and control individuals. Stool samples were collected from all participants, stored the samples at 4 °C until they could be delivered to the biobank. Blood and stool samples collected in each hospital were then sent to the Basque Biobank. Samples were aliquoted samples and frozen at −80 °C. The hospitals clinicians (neurologists and rheumatologists) were responsible for the fibromyalgia diagnosis. The following criteria were used:

Fibromyalgia group: WPI ≥ 7 and SST (Severity Score) ≥ 5 or WPI between 3 and 6 and SST ≥ 9. Patients with other diseases with similar symptoms were discarded.

Control group: healthy individuals without any clinical manifestation of fibromyalgia and/or any other similar disease. To reduce the potential confounding factors associated with lifestyle, they also were age-paired with the patient group and came from the same environment.

All donors signed the informed consent form, and the study was approved by the appropriate ethical committee (CEIC-PI2016037). DNA from faeces was extracted using PSP Spin Stool DNA Plus kit (STRATEC Molecular®), following the manufacturer's protocol. Lysis buffer was added to the frozen samples, to ensure the preservation of nucleic acids. DNA extractions were then aliquoted into samples of 2.5 μg of DNA at the concentration of 100 ng/μL and then frozen until sequencing. All sample processing and distribution were managed by the Basque Biobank. The summary of the collection workflow can be found in Fig. 1.

Fig. 1

Experimental design workflow, from patient recruitment and sample collection to the arrival of processed samples into the research centre and their examination using distinct omics techniques.

V3–V4 16S rDNA sequencing

DNA amplicon libraries were generated and sequencing performed following the recommendations of Illumina Inc. Sequencing was conducted at the FISABIO Sequencing Core Facility, as were the quality assessment using prinseq-lite [87] and the sequence joining, employing FLASH software [53] with default parameters. The complete protocol can be found in the Supplementary Methods file.

Microbiome sequences bioinformatics analysis

QIIME2 package (v. 2017.10) [12] was employed to perform the Operational Taxonomic Units (OTU) clustering and identification, using de novo methodology at 97% similarity threshold. Diversity analysis was performed, and the OTUs were annotated using GreenGenes 13_8 database. The OTU table was exported to SIMCA-9+ 12.0.1 (Umetrics AB, Umeå, Sweden) to perform multivariate analysis and to R programme (R Development Core Team [108]; http://cran.r-project.org) to conduct the statistical analysis using phyloseq [60], microbiome [48] and DESeq2 [52] R packages. CORBATA [49] approach was used to identify and plot the bacteria in the core microbiome. SIAMCAT tool [107] was used to assess the potential effects of confounding factors such as sex, different hospitals and distinct drug types. The adjusted p-value < .05 was considered statistically significant unless stated otherwise. The complete protocol can be found in the Supplementary Methods file.

qPCR validation

From the glutamate cytoplasmatic incorporation and degradation pathways we selected four genes (gadC, glnA, glsA and glsB) to validate our findings related to glutamate and microbiome interaction. We designed specific primer pairs with Primer-BLAST from NCBI webtool (https://www.ncbi.nlm.nih.gov/tools/primer-blast/) indicating specificity for five bacterial families: Bacteroides, Bifidobacterium, Eubacterium, Lachnospiraceae and Ruminococcaceae. Complete procotol can be found in the Supplementary Methods file.

Metabolomics methodology

To 40 μL aliquots of human serum, 40 μL of water/0.15% formic acid (FA) was added. Then, the proteins were precipitated by addition of 120 μL of acetonitrile. To achieve the optimum extraction, after the addition of acetonitrile, the samples were sonicated for 15 min and agitated at 1400 rpm for 30 min (at 4 °C). Next, they were centrifuged at 14,000 rpm for 30 min at 4 °C. The supernatants were transferred to vials. Samples were examined using a UPLC system (Acquity, Waters Inc., Manchester, UK) coupled to a time-of-flight mass spectrometer (ToF MS, SYNAPT G2, Waters Inc.). A 2.1 × 100 mm, 1.7-μm BEH amide column (Waters Inc.), kept at 40 °C, was used to separate the analytes before the MS. The MS was operated in positive electrospray ionisation full scan mode. Spectral peaks were automatically corrected for deviations in the lock mass. The complete specifications can be found in the Supplementary Methods.

Scaled and normalised data were uploaded to R. Principal Component Analysis (PCA) was performed to check whether the differences between sample metabolomes were due to sample origin and to account for the autoclaving process used by one of the hospitals. We excluded the metabolites whose expression differed between the hospitals, to avoid the bias introduced by the sample origin. Metabolomic features with >30% of missing values in either hospital were removed from the analysis. Fold changes and p-values (adjusted using the Bonferroni method) were computed. Afterwards, differential peaks were selected for further annotation and metabolite identification using the METLIN database [29]. The identification was confirmed using commercial standard injection.

MetScape [44] and Ingenuity Pathway Analysis® were used to map the identified metabolites to corresponding functionalities in humans.

Cytokine and miRNA profiling

The cytokine and miRNA profiling was performed by Abcam FirePlex Service (Boston, USA). The cytokine analysis was conducted using the FirePlex Human Discovery Cytokine Panel (Abcam, ab227936), allowing simultaneous profiling of 70 targets in a single well. Each sample was analysed in duplicate, following the manufacturer's instructions. The flow cytometer output was analysed using the FirePlex™ Analysis Workbench software (http://www.abcam.com/FireflyAnalysisSoftware). Cytokine concentration in a sample was interpolated from the standard curve obtained in duplicate for each plate. The data was log-normalised, and then the fold changes and Bonferroni-adjusted p-values were computed to assess the differences between the cytokine profiles.

The miRNAs were profiled using the FirePlex miRNA Assay Core Reagent Kit (Abcam, ab218342) employing a custom multiplex panel with 68 miRNAs selected on the basis of literature review. Each sample was run in singlicate, as previously described [93]. Data analysis was performed using the FirePlex™ Analysis Workbench software. Three miRNAs used for normalisation, hsa-miR-17-5p, hsa-miR-320b and hsa-let-7i-5p, were selected employing the geNorm algorithm [96]. The data was log-normalised, and then the fold changes and Bonferroni-adjusted p-values were computed to evaluate the differences between the miRNA profiles.

Omics integration

Microbiome and metabolomics

Spearman's correlation coefficients were computed for relationships between relative abundances of microbiome bacteria with the identified genus and normalised individual metabolomic features. A scaled heatmap was constructed for the correlation matrix, including cladogram classification of the variables, using the default clustering method.

Integration of all datasets

We employed the Data Integration Analysis for Biomarker Discovery (DIABLO) using Latent cOmponents implementation in the mixOmics R package [79,90]. Thirty-six fibromyalgia and 35 control samples were used. Microbiome data was normalised using DESeq2 counts function. The mixOmics block.splsda function, with full weighted design and 10 components, was primarily used to identify the optimal number of components, which was defined in 3 methods using the centroid distance technique. To decide which variables to keep in each component, models with 10, 5, 5 and 5 randomly selected variables were tested for the microbiome, metabolomics, cytokines and miRNAs, respectively. Finally, different model features were obtained and the results were plotted using mixOmics predefined and ad-hoc functions. This procedure was followed for both the identified-metabolite dataset and the full dataset of unidentified metabolomics features.

Results

Clinical samples

One hundred and five confirmed fibromyalgia patients (ACR 2010 modified criteria) [102] and 54 age- and environmentally-paired healthy individuals were recruited. The latter group consisted of individuals who did not present any disease or symptoms related to fibromyalgia and came from the same environment as the fibromyalgia patients. The characteristics of the study cohort are shown in Table 1.

Table 1

Cohort characteristics. The number of individuals included in each group is given in parentheses. For Age, WPI and SST, mean values ± standard deviation are shown.

	Controls (n = 54)	Fibromyalgia-diagnosed patients (n = 105)
Sex	48.15% ♀, 51.85% ♂	69.52% ♀, 30.48% ♂
Age (years)	53.5 ± 12.4	52.52 ± 10.3
Age at diagnosis (years)	NA	48.2 ± 11.1
Time since diagnosis (years)	NA	3.4 ± 6
WPI	NA	13.28 ± 3.91
SST	NA	8.62 ± 2.15
SS1	NA	6.6 ± 1.8
SS2	NA	2.1 ± 0.4

During WPI evaluation, >90% of the patients reported pain in the back, shoulder girdle and the abdomen. Neck pain was described by 85% of patients, while upper and lower arm, hip and upper and lower leg pain were reported by 70% of fibromyalgia patients. At least 50% of the patients were affected by jaw and chest pain. The SST index is the combination of two sub-indexes, SS1 (the severity of 3 main symptoms in fibromyalgia: fatigue, sleep quality and cognitive problems) and SS2 (the list of associated fibromyalgia symptoms). Approximately 90% of patients reported moderate to severe scores for the 3 main symptoms for the SS1 sub-index in the week preceding the collection of the samples. In the evaluation of associated fibromyalgia symptoms (SS2), 70.7% of fibromyalgia patients presented at least 4 symptoms from the neurological sphere (muscle pain, fatigue, thinking or memory problems, headache, numbness/tingling, etc.). Among them, 70% used painkillers, while approximately 55% were taking antidepressants and benzodiazepines and approximately 30%, antiepileptic drugs (SUPPLEMENTARY TABLE S1). Half of the patients reported some physical exercise and some alcohol consumption, while 23% identified themselves as smokers.

V3 + V4 16S rDNA sequencing

We obtained 6,110,564 reads, of which 99.56% passed the quality check. Of the cleaned reads, the 81.91% (4,982,956) were joined. To decide on the number of reads to which the samples should be rarefied; we computed the rarefaction curves for both observed OTUs and Shannon indices (Supplementary Fig. S1A). After rarefying at 12,000 reads/sample, the median coverage was 96.35 ± 2.33%. Rarefaction step did not reduce diversity (Supplementary Fig. S1B). Sequencing data was uploaded to ENA under Project Accession code PRJEB27227.

Microbiome analysis

The multivariate unsupervised PCA (Fig. 2A) did not show any differences between the control and the fibromyalgia samples. The supervised Partial Least Squares Discriminant Analysis (PLS-DA), however, provided two sample groups (Fig. 2B) (p-value, 0.0019). In the specific diversity analysis for 4 alpha-diversity indices (Faith's Phylogenetic Distance, ace, chao1 and observed OTUs) we observed a discrete decrease in bacterial diversity in fibromyalgia patients although only the Faith's PD index showed a statistically significant difference (Fig. 2C). This reduction in bacterial diversity was also observed in the analysis of the core microbiome at the taxonomic family level. We used CORBATA default parameters (80% ubiquity, 1% abundance) to identify which bacteria families present in both fibromyalgia and control core microbiomes. The two core microbiomes contained the same 4 bacteria families (C. Ruminococcaceae, C. Lachnospiraceae, B. Rikenellaceae and B. Bacteroidaceae). We observed that the control group presented a more diverse bacterial community. The comparison of the two sample groups revealed that Clostridiales Ruminococcaceae was more abundant in the healthy control group than in fibromyalgia patients, although the differences were not statistically significant (Fig. 3A). After reducing the cut-off to 50% ubiquity, we observed differences between the core microbiomes of the two groups. Specifically, two bacteria families that were absent in the fibromyalgia core microbiome, the Bifidobacteriales Bifidobacteriaceae and the Bacteroidales Prevotella, which were represented in the control core microbiome (Supplementary Fig. S2A).

Fig. 2

Microbiome multivariate analysis. (A) Principal Component Analysis (PCoA) of the complete cohort. (B) Supervised Partial Least Squares Discriminant Analysis (PLS-DA) analysis, showing the discrimination between the sample groups. (C) Alpha-diversity indexes for each sample group, showing the adjusted p-value computed using Student's t-test.

Fig. 3

Core microbiome and genus-discriminant analyses. (A) The composition of core microbiome for each sample group and the comparison of bacterial ubiquity in the two groups. (B) Genera significantly different (adj p > .05) between the control and fibromyalgia samples, obtained using the protocols described in the Methods. Positive log2 fold changes (x-axis) indicate genera with positive fold difference between fibromyalgia and control. Negative log2 fold changes are shown as negative x values. Each point represents a single OTU, coloured by phylum. On the y-axis, the taxonomic genus level is indicated. Size of the points reflect the log-mean abundance of the sequence data. (C) qPCR results for the differential expression of bacterial genes related to glutamate bacterial degradation. Results are indicated in differential Cts count.

We performed a differential OTU analysis (employing DESeq2) of the core microbiomes in the control and fibromyalgia samples. We identified 32 OTUs distributed among 3 phyla (Actinobacteria, Bacteroidetes and Firmicutes) (Fig. 3B) whose abundance differed between the two groups, with an adjusted p-value of 0.05. In fibromyalgia patients, the Bacteroidetes and Firmicutes had OTUs both with increased and decreased abundance, and Actinobacteria levels were reduced in this group (Fig. 3B).

The number of OTUs with the unassigned genus in Bacteroidaceae and Lachnospiraceae families were decreased in fibromyalgia samples; there were also fewer Bifidobacteriaceae and Erysipelotichaceae OTUs in fibromyalgia patients. The Rikenellaceae family showed an increased abundance in fibromyalgia patients (Supplementary Table S2).

Finally, at the genus level, the abundance of Bacteroides OTUs was reduced in fibromyalgia patients, as were Bifidobacterium, Eubacterium and Clostridium OTUs. However, the abundances of the genera Dorea, Roseburia and Alistipes were increased in this group (Fig. 3B).

There were no significant differences between microbiome composition abundances in the two sexes. We did not observe any significant association between drug types (as summarized in Supplementary Table S1) and the relative microbiome abundance at the genus level (data not shown).

We validated the reduction of the abundance of bacterial species by qPCR technique. For that, we amplified a set of genes dedicated to the glutamate incorporation to bacterial cytoplasm and its transformation to GABA (gadC, glnA, glsA and glsB). We designed specific primers for amplifying genes from 5 bacterial families that we found to be diminished in fibromyalgia patients (Bacteroides, Bifidobacterium, Eubacterium, Lachnospiraceae and Ruminococcaceae) (Fig. 3C). We found that the gene encoding the transporter of glutamate into bacterial cytoplasm, represented by gadC, was diminished, as it was also the genes encoding enzymes involved in the transformation of glutamate to l-glutamine (glnA, glsA) and to GABA (gadB) (Supplementary Fig. S2B), in agreement with the taxonomic analysis of 16S rDNA gene.

Metabolomics analysis

The metabolomics analysis yielded 8543 different metabolic features defined by retention time and mass/charge. One sample was removed due to technical failure. The PCA analysis revealed that the metabolomics profiles differed between hospitals (Supplementary Fig. 3). This was expected because of the autoclaving performed in one of the hospitals. Thus, to avoid the bias caused by the chemicals released during the autoclaving procedure, the discriminating hospital features (p = 661), were removed from the study, as well as the features with >30% of missing values. Two hundred and twenty-eight features differed between the fibromyalgia and control groups (Fig. 4A). Of these 228, only 88 had tentative IDs in the METLIN database. Using MS/MS data and chemical standards, we found that the levels of 7 of these metabolites were significantly altered in the fibromyalgia samples (Supplementary Table S3): ornithine, L-arginine, Nε-Methyl-l-lysine, l-glutamate, l-glutamine, asymmetric dimethylarginine (ADMA) and platelet activating factor (PAF-16) (Fig. 4B). Another metabolic feature among the 228 altered in fibromyalgia was tentatively identified as l-threonine or DL-homoserine (Fig. 4B). We could not discriminate between these two metabolites as they are structurally similar and have the same molecular mass and fragmentation pattern in LC-MS. We also analysed the metabolites described in the literature, such as creatinine [31,55], platelet activating factor [11] and acetylcarnitine [25]. To infer alterations in the biological processes and metabolic and functional pathways associated with the differentially expressed metabolites, we used MetScape [44] and Ingenuity Pathway Analysis® (QIAGEN) (IPA). The analyses showed that cell signalling and inflammatory and hypersensitivity responses were the most relevant biological processes. The most represented metabolic pathways were arginine, nitric oxide (NO) and glutamate metabolism.

Fig. 4

Univariate metabolomics analysis. (A) Volcano plot of 1070 metabolic features detected in serum samples after background subtraction and removal of the features found in <30% of the data or differing between hospitals. (B) Volcano plot of the identified metabolites. Positive log2 FC indicates increased abundance in fibromyalgia patients. All p-values were adjusted using the Bonferroni method.

To study the potential dependencies between microbiome composition and the host metabolism and metabolome, we examined the correlations between the two datasets. We computed the Spearman's correlation coefficient for the full set of metabolomics features and microbiome variables. We did not see any clear association patterns between the two complete datasets (Supplementary Fig. S4). We also constructed a heatmap of the scaled correlations between the bacteria whose abundance was changed in fibromyalgia and the identified metabolites (Fig. 5A). Metabolites were grouped into two clusters, depending on the correlations. These were seen mainly with for genera Bifidobacterium and Dorea, which behaved in an opposite manner. The first cluster contained 4 metabolites (3-methyl-l-Lysine, PAF C-16, ADMA, l-Lysine). The second cluster was formed by 8 metabolites (glutamate, l-threonine/DL-homoserine, glutamine, Nε-methyl-l-Lysine, creatinine, ornithine, arginine and acetylcarnitine), although the metabolite acetylcarnitine behaved differently from the other metabolites in this cluster. Bifidobacterium, whose abundance was reduced in fibromyalgia patients, correlated negatively with the first metabolite cluster and positively with the second one. Dorea, with increased abundance in fibromyalgia patients, correlated positively with the first metabolite cluster and negatively with the second one.

Fig. 5

Heatmap of scaled correlations between the bacteria whose abundance was altered in fibromyalgia and the identified metabolites. The dendrograms were unsupervised. Red arrows mark the bacteria with increased abundance in fibromyalgia, green arrows, with decreased abundance, and “equals” symbol indicates the OTUs with both increased and decreased abundance (A). Omics correlations with indexes used in fibromyalgia diagnostics, as defined by ACR 2010 criteria. Only significant correlations (p-value < .05) are coloured. Positive correlations are indicated in red and negative correlations, in blue. Correlations between circulating miRNA levels (B), circulating cytokine levels (C), identified serum metabolites (D) and microbiome composition (at genus level) (E). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Finally, we checked, using Virtual Metabolic Human [65] database, whether the different metabolites were produced by the differentially abundant bacteria. We also wanted to study whether they were made by the genera for which we found most correlations (Fig. 5A). Thus, we limited the search to Bifidobacterium and Dorea genera. For glutamate, we identified the metabolites upstream and downstream of its production/degradation. For lysine, threonine, homoserine, glutamine, ornithine and arginine (and their modifications), we found that the metabolites themselves, their precursors and degradation products might had been produced by bacteria. No bacterial associations were found for creatinine, PAF C-16, ADMA and acetylcarnitine, consequently suggesting that their origin was exclusively human.

Serum factors and miRNA analysis for a subset of samples

A subset of the samples (n = 72; nC = 36 controls and nF = 36 fibromyalgia samples) was used to perform multiplex assays for different serum molecules, including miRNAs and cytokines. For the multiplex design, we used 70 molecules and 68 miRNAs that have been associated with fibromyalgia and/or chronic pain. The protein content assays and the miRNAs analysis did not show any differences between the fibromyalgia and the control groups. We observed statistically significant differences for ten serum proteins: PCSK9, mesothelin, BST2 (↑), procalcitonin, Axl, myoglobin, MIG, TNF-alpha, ICAM2 and IL-9 (↓) with fold changes ranging from 0.76 (lower level in patients) for IL-9 to 1.07 for BST2 (Supplementary Fig. S5A). However, the levels of only one miRNA differed significantly between the fibromyalgia patients and the control group, the hsa-miR-335-5p (Supplementary Fig. S5B). Predicted target genes were obtained using miRWalk 2.0 database [20]; they were selected if they mapped to at least 8 of the 12 database options. The enrichment of the miRNA targets was performed using ConsensuspathDB [42], selecting the targets with a p-value < .01. Notably, we identified several pathways related to signalling dedicated to gene regulation processes. The complete results are provided in SUPPLEMENTARY TABLE S4.

Correlations between omics data and clinical data

To determine which differences could be associated with the disease, we examined the correlations between different diagnostic indexes obtained for the fibromyalgia patients and the omics data (Fig. 5B, C, D and E). Notably, miRNA data constituted the omics dataset most correlated with pain indicators (Fig. 5B), followed by the results of serum protein profiling (Fig. 5C). Metabolomics also showed a considerable number of correlations with several pain indexes (Fig. 5D). The microbiome composition (at the genus level) (Fig. 5E) was the omics dataset with the weakest correlation with pain indicators.

We also considered possible effects of medication on the observed differences between the patient and control samples. We checked whether the samples clustered depending on the drug regimen followed. However, we did not find any clusters of samples (neither for serum factors nor for miRNAs) that could be associated with a specific drug or drug combination. We also checked whether any data correlated with distinct drug types; no such correlation was observed (data not shown).

Modelisation of microbiome, metabolomics, cytokine and miRNA datasets

We combined the four datasets of the 71 samples (nC = 36, nF = 35) that had all the data. This allowed us to discriminate between the control and fibromyalgia samples when a block sparse PLS-DA model was applied (block sPLS-DA) (Fig. 6A). The analysis of the individual contribution of each dataset to the differences showed that the most correlated datasets were the microbiome composition and the metabolomics data. We also found that the major contributor to the separation of the sample groups was the microbiome dataset, followed by serum metabolomics, proteins and, finally, miRNAs (Fig. 6B and C). In this analysis, we used only the metabolomics dataset containing the identified metabolites (n = 14). The sPLS-DA analysis using the whole unidentified metabolomics dataset (n = 1070) showed that using the metabolomics dataset improved the discrimination between the two sample groups, becoming the strongest factor distinguishing the patients from controls (Supplementary Fig. S6) although the microbiome showed slightly better predictive ability.

Fig. 6

Multi-omics integration. (A) sPLS-DA consensus plot for the combination of the 4 datasets, showing the nearly complete discrimination of the 71 samples (36 fibromyalgia and 35 control samples). (B) The individual contribution of each dataset to the sPLS-DA final model, in each case showing the score plots for the two first components, indicating the best separation capability for microbiome data, followed by cytokines, metabolomics and miRNAs. (C) ROC curves for each omics dataset, with the Area under the Curve (AUC) values.

Discussion

In this study, we applied an omics approach and identified a set of potential molecular markers (Table 2) for the diagnosis of fibromyalgia.

Table 2

Differences between fibromyalgia and healthy control groups observed using each omics technique (showing alterations in the fibromyalgia patients).

	Increased (↑)	Decreased (↓)
Microbiome	DoreaRoseburiaPapillibacterSubdoligranulum	BifidobacteriumEubacteriumLachnospiraceae (family)ClostridiumFirmicutes (phylum)
Metabolomics	l-glutaminel-threonine/DL-homoserinel-arginineADMAl-glutamateNε-methyl-l-lysineOrnithine	PAF-16
Cytokines	PCSK9MesothelinBST2	ProcalcitoninAxl-UFOMyoglobinMIGTNF-alphaICAM2IL-9
miRNAs	hsa-miR-335-5p

The gut microbiome analysis revealed two clusters (Fig. 2B), one cluster for fibromyalgia patients (modified 2010 ACR diagnostic criteria) and the other for individuals without any clinical manifestation of fibromyalgia. Both core microbiome and alpha-diversity analyses showed a reduction in the bacterial diversity in the fibromyalgia group. This is in agreement with the report of reduced microbiota diversity in other pain disorders, such as myalgic encephalomyelitis/chronic fatigue syndrome [27]. Interestingly, our fibromyalgia microbiome analysis showed a reduction in the abundance of several bacterial strains associated with healthy microbiome, such as those linked to SCFA production (Bifidobacterium, Eubacterium and Lachnospiraceae) [40,64,77,94], and/or to the reduction in Firmicutes phylum OTUs ([75]; Human Microbiome Project Consortium et al., 2012; [51]), suggesting dysbiosis events in fibromyalgia patients. Currently, there is no consensus on the use of the term “dysbiosis” or its meaning [35]. Thus, we would like to clarify that we refer to alterations in microbiome composition linked to disease (either causing the disease or appearing as its consequence). Dysbiosis events are also associated with the disruption of the intestinal barrier; this increases the interactions of bacteria with the immune system of the host, producing local inflammation [41]. This is supported not only by the large proportion of patients reporting abdominal pain (>90%) but also by the number of intestinal diseases considered co-morbidities of fibromyalgia. The maintenance of the intestinal barrier is associated with the production of SCFAs, including butyric acid and butyrate [77]. In fibromyalgia, we found a decrease in the abundance of several members of the Lachnospiraceae family, the bacteria involved in butyric acid production [61]. Butyrate, the conjugate base of butyric acid, is produced by a small number of bacteria, including several Eubacterium species [64], a genus also underrepresented in fibromyalgia patients. The reduction in the diversity of bacteria, especially of those engaged in the production of protective SCFAs, suggests that this process might be implicated in the development of fibromyalgia. If this is the case, the dysbiosis events, as understood here, should be persistent. Thus, we recognise that multiple time-point data should be acquired and studied; lack of this data is a limitation of our study. We would like to emphasise that this is a pilot study and that a follow-up analysis, which might reinforce our findings, is recommended.

We also found differences between neurotransmitter metabolisms in the patients and control individuals. We detected a significant increase in the serum levels of glutamate in fibromyalgia patients. Moreover, the abundance of bacteria from Bifidobacterium and Lactobacillus genera (involved in the transformation of glutamate into GABA; [4,8,105] was reduced in the fibromyalgia group. This might contribute to the elevated systemic levels of glutamate. The effect of GABA on the gut-brain axis, via the vagus nerve, has been described by several authors [8,16]. Glutamate affects the development of pain, via glutamatergic synapses [69], and stress can alter the regulation of this pathway [74]. Stress-related events have also been associated with microbiome modifications [8]. The 2010 modified ACR criteria for fibromyalgia diagnosis include several stress-associated symptoms. Whether such elevated systemic levels of glutamate affect the ENS and alter the CNS is still unclear. However, some authors have demonstrated the activation of glutamatergic neurons and glutamate-mediated neurotransmission in the ENS [13,46,50,84]. As a result of a reduction in bacterial diversity, the glutamate might enter the host bloodstream after the disruption of the intestinal barrier by the inflammation caused by the dysbiosis. Interestingly, several patients presented with symptoms associated with IBD as fibromyalgia co-morbidities (irritable bowel syndrome (46%), abdominal pain (13%) and the pain in the upper abdomen (45%), diarrhoea (20%), etc.). The role of microbiome in IBD pathogenesis has been broadly demonstrated [23,86]; a dysregulation of intestinal immune system caused by microbiome alterations may lead to disease [91], as demonstrated by patients presenting T-cell responses against commensal bacteria [73]. Specifically, a reduction in the abundance of Firmicutes phylum bacteria (observed in fibromyalgia patients) has been recurrently associated with IBD pathogenesis and progression [24,63]. These common alterations in microbiome composition could explain some of the most frequent co-morbidites reported by the patients in our study.

Furthermore, it has been shown that the blood-brain barrier increases its permeability after a decrease in the numbers of SCFA-producing bacteria. This alters the tight junction organisation, which can be recovered by colonisation with SCFA-producing bacteria and/or by the administration of these bacterial metabolites [7]. Cytokines can also modify the blood-brain barrier permeability [6,103]. Importantly, glutamate levels increase in the cerebrospinal fluid (CSF) of fibromyalgia patients [85]. These data suggest an important role of this neurotransmitter in the pathogenesis of fibromyalgia. The manner in which the peripheral levels of gut microbiome derived neurotransmitters can affect the brain function is still under debate [84], although several mechanisms have been proposed. Alterations in the blood-brain barrier permeability could modify the interchange of serum metabolites with the brain. Serum levels of 5-HTs are altered in germ-free mice [101,104]. Even though 5-HT itself is not known to cross the blood-brain barrier, its precursor can. The microbiome might alter the 5-HT precursor (e.g. tryptophan) levels, as has been proposed by several authors [67,88]. The same mechanism has been suggested for other gut microbiome neurotransmitters, such as dopamine and GABA [56,84,97].

It is essential to keep in mind the relationship between GABAergic pain inhibition and gender as fibromyalgia is 3 times more prevalent in women than in men [76]. Steroid 17β-estradiol (E2) suppresses the GABAergic inhibition in female rats via a sex-specific oestrogen receptor ERα, mGluR and endocannabinoid-dependent mechanism [92]. This suppression requires the activation of mGluR type I receptors by glutamate [36]. Therefore, in the presence of excess glutamate, as observed here in fibromyalgia patients, the pain inhibition by GABA might be suppressed in female patients by this E2-specific regulation. This might partly explain the increased prevalence of fibromyalgia in the female population.

The functional analysis of the metabolomics dataset showed that the most represented pathways were those dedicated to the metabolism of known neurotransmitters, such as glutamate and serine. Both arginine and ornithine levels, related to the widespread pain in fibromyalgia, increased in the sera of fibromyalgia patients. Consistently, IPA analysis identified several pathways related to arginine, such as arginine degradation (I and II) canonical pathways and proline biosynthesis from arginine. These two metabolites are required for the synthesis of nitric oxide (NO) [31]. NO plays an important role in both acute and chronic pain as it is a mediator of nociception [17]. However, NO contributes not only to nociception; it also mediates in analgesia and increases the effect of morphine on pain inhibition [17]. Here, we also observed a strengthening of this pathway in fibromyalgia patients (by using IPA). The role of NO in fibromyalgia pathogenesis has been studied but without reaching a consensus [72]. Notably, the levels of iNOS isoform increase in female fibromyalgia sufferers in comparison with healthy controls, while the levels of constitutive isoforms (nNOS and eNOS) do not change [59]. It is important to remember that our functional profiling was performed using the results obtained from the serum sample analysis. One of the limitations of this study is the metabolomics analysis, and specifically, the metabolite identification step. We could only identify a small subset of all the metabolic features observed. Thus, the results obtained here are constrained by the relatively small number of identified metabolites. An improved metabolite identification procedure could not only expand the list of potential metabolite biomarkers but also advance the identification of potentially affected biological pathways and functionalities.

Patients afflicted by chronic pain are likely to participate in many different long-term treatments, which could affect their microbiome composition. Differences in diets and lifestyles will also have some effect. Thus, it is difficult to be certain whether the detected alterations in the microbiota are the cause or consequence of fibromyalgia. No association between microbiome composition and drug type was found for fibromyalgia patients. However, it has been demonstrated that clinical drugs have an impact upon microbiome composition; this seems to be true for antibiotic, non-antibiotic [54] and psychotropic [19] drugs. The lack of associations shown here could have been caused by the small number of patients taking medication from a specific drug family and/or by the interactions between different drugs prescribed. Proton pump inhibitors (PPI), for example, have an anticommensal activity and were taken by nearly 30% of the patients. One study has reported a reduction in the abundance of Lachnospiraceae and Ruminococcaceae in PPI consumers [39], which is consistent with our observations for fibromyalgia patients. Another study obtained similar results and considered in its analysis the decrease in the abundance of Bifidobacterium genus in PPI consumers [38]. Both studies have reported a decrease in α-diversity after PPI administration, which is also consistent with our findings. It has been reported that psychotropics target a similar pattern of bacterial species irrespective of the degree of their chemical similarity. This suggests that the anticommensal activity of these drugs is a part of their mechanism of action rather than a secondary effect [19].

We did not observe any microbiome alterations that could be associated with antidepressant drugs, either for the tricyclic antidepressants (taken by 12% of patients) or for the selective serotonin reuptake inhibitors (SSRI), 54% of patients). The antiepileptic drugs (here taken by 29% of patients), such as lithium or valproate, do not have a significant anticommensal activity. However, lithium may increase the relative abundance of Ruminococcaceae and reduce the abundance of Bacteroides, while valproate alters the levels of SCFA [18]; there were also alterations found in fibromyalgia patients. Finally, while no antimicrobial activity has been reported for morphine [83], chronic use of opioids (prescribed to 45% of patients) has been associated with a reduction in Bacteroidaceae (which we also observed in fibromyalgia patients) and Ruminococcaceae [3]. Even though no associations between specific drugs and microbiome composition was found, probiotics could be useful in the treatment of fibromyalgia as they affect the microbiome composition [34]. Notably, several authors have used this approach to treat the chronic fatigue syndrome [82] and one pilot study has examined the effects of probiotics on fibromyalgia patients [80]. The authors have shown some improvements, mainly in depression symptoms and impulsive behaviour, in comparison with the placebo group [81].

Conclusions

To the best of our knowledge, this is the first study to report differences between the microbiome composition of fibromyalgia patients and healthy controls. We provided a list of these differences and reported the alterations in the levels of various molecules in the fibromyalgia sufferers, which might be useful as diagnostic biomarkers. We examined the functionality of these molecules and found that the most altered metabolic pathways were related to neurotransmitters, such as glutamate and nitric oxide. We checked possible interactions between the gut microbiome and serum metabolome; our analysis found several individual correlations between the two datasets. We also demonstrated that the combined microbiome and serum metabolome analyses could discriminate between the fibromyalgia patients and control individuals. Thus, we report a new set of molecules and bacteria that might improve the diagnosis process, compensating for the current lack of objective biomarkers. Our results should help to shed some new light on the pathogenesis of this disease, provide biomarkers within a biological framework and improve our knowledge of this relatively unknown disease.

Ethics approval and consent to participate

All donors signed the informed consent form, and the study was approved by the ethical committee (CEIC-PI2016037).

Availability of data and materials

Sequencing data was submitted to the ENA repository, under the Project Accession code PRJEB27227.

We used free, open-access software, except for SIMCA-P+ 12.0.1 (Umetrics AB, Umeå, Sweden) and Ingenuity Pathway Analysis® from QIAGEN, which require a license. Further information and requests for resources and reagents should be directed to and will be fulfilled by the Lead Contact, Juan M. Falcón-Pérez (jfalcon@cicbiogune.es). Sharing of patient data may be restricted due to anonymity considerations.

Funding

This project was funded by the Basque Government's Health Department (ref. 2015111149 (2016-2018)).

Author contributions

Experimental design: Marc Clos-Garcia, Naiara Andrés, Gorka Fernández-Eulate, Olga Maíz, Adolfo López de Munain and Juan Manuel Falcón-Pérez.

Patient recruitment and sample collection: Naiara Andrés, Gorka Fernández-Eulate, Alejandro Valero, Nerea Errazquin, María Cristina Gómez Vallejo, Leila Govillard, Ana María Callejo Orcasitas, Olga Maíz and Adolfo López de Munain.

Sample processing and data analysis: Marc Clos-Garcia, Leticia Abecia, José L. Lavín, Sebastiaan van Liempd, Diana Cabrera, Félix Royo, Esperanza Gónzalez, Ana M. Aransay, Michael R Tackett and Genesis Tejada.

Manuscript writing: Marc Clos-Garcia, Leticia Abecia, José L. Lavín, Naiara Andrés, Gorka Fernández-Eulate and Juan Manuel Falcón-Pérez.

Discussion and review of the manuscript: Marc Clos-Garcia, Leticia Abecia, José L. Lavín, Juan Anguita, Luis Bujanda, Ana María Callejo Orcasitas, Ana M. Aransay, Olga Maíz, Adolfo López de Munain and Juan Manuel Falcón-Pérez.

The following are the supplementary data related to this article.

Rarefaction curves computed for observed OTUs (upper panel) and Shannon (lower panel) indexes. The vertical line at 12,000 reads indicates the value selected for rarefying the samples (A). PCoA analysis of unweighted UNIFRAC distances of the non-rarefied and the rarefied datasets. Axis titles indicate the variance explained by each component. Control samples are coloured in blue and fibromyalgia in red (B).

Glutamate incorporation to bacterial cytoplasm and transformation pathways. Decreased genes are indicated with a red triangle.

PCA analysis of full metabolomics data, coloured by the hospital sample origin.

Scaled correlations between the full metabolomics dataset and the microbiome OTUs, at the genus level. Red colouring indicates positive correlations and blue, negative correlations.

Volcano plots for the cytokine (A) and miRNA (B) profiling. All p-values were adjusted using the Bonferroni method.

Multi-omics integration, employing DIABLO functionality in the mixOmics R package, using the full metabolomics feature dataset. (A) sPLS-DA consensus plot for the combination of the 4 datasets, showing the perfect discrimination between the two sets of samples (36 fibromyalgia and 35 control samples). (B) The individual contribution of each dataset to the final sPLS-DA model, showing the score plots for the two first components in each case. Using the metabolomics data gives the best separation capability, followed by the microbiome, cytokines and miRNAs. (C) ROC curves for each omics dataset, showing the AUC value for each set.

Table S1

Summary of the drug regime followed by patients (sheet 1) and statistical distribution of the drugs grouped by function and by active ingredient.

Table S2

Significant DESeq2 results for the comparison of OTUs in the control and fibromyalgia groups. OTU names (see Fig. 3) and random OTU classifiers obtained from GreenGenes database are shown.

Table S3

Full metabolomics features, after removing those not associated with any ID in the METLIN database. Retention time, imputed mass, fold change, p-value, tentative ID and validation indication are shown.

Table S4

hsa-miR-335-5p mRNA targets obtained from miRWalk 2.0 (Sheet 1) and pathway enrichment results obtained from ConsensusPathDB (Sheet 2). In blue, the mRNA targets selected for the posterior enrichment analysis in ConsensusPathDB.

Supplementary material

Declaration of Competing Interest

The authors declare no competing interests.