Gut bacteria and blood and urine chemicals linked to vestibular migraine

Demographic information [age, sex, body mass index (BMI)] and disease-related information (disease duration, attack frequency, clinical manifestations during attacks) were collected using Excel spreadsheets. Migraine patients were assessed using the Visual Analog Scale (VAS) for headache severity, and VM patients were assessed using the Dizziness Handicap Inventory (DHI) scale.

Fecal sample collection occurred within 24 h of the collection of demographic/disease information and plasma/urine samples. Participants were instructed to collect fecal samples in the hospital according to a standardized protocol. Prior to collection, participants were asked to empty their bladder. Approximately 3 g of fecal material was collected from the middle inner part of the stool using a uniform 5 ml sterile fecal collection device. The sample was immediately labeled with a marker pen, placed on ice, sealed with parafilm within 2–4 h, and stored at −80 °C until analysis.

Plasma sample collection occurred within 24 h of the collection of demographic/disease information and fecal/urine samples. After an overnight fast, 5 ml of peripheral blood was collected from each participant at the outpatient department or Neurology Ward IV into EDTA anticoagulant tubes. Samples were transported on ice. Plasma was separated by centrifugation at 3,000 rpm for 10 min at room temperature within 2–4 h of collection. The supernatant was aliquoted into 1.5 ml centrifuge tubes (0.2 ml per tube), labeled, sealed with parafilm, and stored at −80 °C until analysis.

Urine sample collection occurred within 24 h of the collection of demographic/disease information and fecal/plasma samples. Participants were instructed to collect mid-stream morning urine (approximately 20 ml) into sterile containers in the hospital according to a standardized protocol. Samples were transported on ice. Urine was centrifuged at 1,000 rpm for 5 min at 4 °C within 2–4 h of collection. The supernatant was filtered through a 0.22 μm membrane filter, aliquoted into 1.5 ml centrifuge tubes (1 ml per tube), labeled, sealed with parafilm, and stored at −80 °C until analysis.

Statistical analyses were performed using IBM SPSS Statistics software (version 25.0). Categorical data are presented as number and percentage [n (%)]. Group comparisons for categorical variables were conducted using the Chi-square test or Fisher's exact test, as appropriate. Continuous data are expressed as mean ± standard deviation (SD) if normally distributed, or as median (interquartile range; IQR) if non-normally distributed. Group comparisons for continuous variables were performed using analysis of variance (ANOVA). When the assumption of homogeneity of variances was violated (assessed by Levene's test), Welch's ANOVA with Tamhane's T2 post-hoc test was used. A P-value < 0.05 was considered statistically significant.

DNA from different samples was extracted using the E.Z.N.A.^® Stool DNA Kit (D4015, Omega, Inc., USA) according to manufacturer's instructions. The reagent which was designed to uncover DNA from trace amounts of sample has been shown to be effective for the preparation of DNA of most bacteria. Nuclear-free water was used for blank. The total DNA was eluted in 50 μl of Elution buffer and stored at −80 °C until measurement in the PCR. The 5′ ends of the primers were tagged with specific barcods per sample and sequencing universal primers. PCR amplification was performed in a total volume of 25 μl reaction mixture containing 25 ng of template DNA, 12.5 μl PCR Premix, 2.5 μl of each primer, and PCR-grade water to adjust the volume. The PCR conditions to amplify the prokaryotic 16S fragments consisted of an initial denaturation at 98 °C for 30 s; 32 cycles of denaturation at 98 °C for 10 s, annealing at 54 °C for 30 s, and extension at 72 °C for 45 s; and then final extension at 72 °C for 10 min. The PCR products were confirmed with 2% agarose gel electrophoresis. Throughout the DNA extraction process, ultrapure water, instead of a sample solution, was used to exclude the possibility of false-positive PCR results as a negative control. The PCR products were purifyied by AMPure XT beads (Beckman Coulter Genomics, Danvers, MA, USA) and quantified by Qubit (Invitrogen, USA). The amplicon pools were prepared for sequencing and the size and quantity of the amplicon library were assessed on Agilent 2100 Bioanalyzer (Agilent, USA) and with the Library Quantification Kit for Illumina (Kapa Biosciences, Woburn, MA, USA), respectively. The libraries were sequenced on NovaSeq PE250 platform. Samples were sequenced on an Illumina NovaSeq platform according to the manufacturer's recommendations. Paired-end reads were assigned to samples based on their unique barcode and truncated by cutting off the barcode and primer sequence. Paired-end reads were merged using FLASH. Quality filtering on the raw reads was performed under specific filtering conditions to obtain the high-quality clean tags according to the fqtrim (v0.94). Chimeric sequences were filtered using Vsearch software (v2.3.4). After dereplication using DADA2, we obtained feature table and feature sequence. Alpha diversity and beta diversity were calculated by normalized to the same sequences randomly. Then according to SILVA (release 132) classifier, feature abundance was normalized using relative abundance of each sample. Alpha diversity is applied in analyzing complexity of species diversity for a sample through 5 indices, including Chao1, Observed species, Goods coverage, Shannon, Simpson, and all these indices in our samples were calculated with QIIME2. Beta diversity was calculated by QIIME2, the graphs were drawn by R package. Blast was used for sequence alignment, and the feature sequences were annotated with SILVA database for each representative sequence.

Additionally, taxonomic classification and subsequent analyses were performed using the SILVA and NT-16S databases. Based on the annotation results of Amplicon Sequence Variants (ASVs) and the ASV abundance table for each sample, species abundance tables at the kingdom, phylum, class, order, family, genus, and species levels were obtained. Compositional analysis and differential abundance analysis were conducted on the abundance tables at different taxonomic levels to identify differentially abundant taxa. This included the Kruskal–Wallis test, with P-values adjusted for multiple comparisons using the Benjamini–Hochberg method for False Discovery Rate (FDR) correction, applying a significance threshold of q < 0.05. For the LEfSe analysis, no results remained statistically significant after multiple testing correction. Therefore, in line with the exploratory nature of this study, results with an LDA score > 3 and an uncorrected P < 0.05 were considered statistically significant for exploratory purposes and are presented in a bar chart. Diagrams were implemented using the R package (v3.5.2).

100 μl of sample was taken, mixed with 400 μl of extraction solution [MeOH:ACN, 1:1 (v/v)], the extraction solution contains deuterated internal standards, the mixed solution was vortexed for 30 s, sonicated for 10 min in 4 °C water bath, and incubated for 1 h at −40 °C to precipitate proteins. Then the samples were centrifuged at 12,000 rpm (RCF = 13,800 ( × g), R = 8.6cm) for 15 min at 4 °C. The supernatant was transferred to a fresh glass vial for analysis. The quality control (QC) sample was prepared by mixing an equal aliquot of the supernatant of samples. For polar metabolites, LC-MS/MS analyses were performed using an UHPLC system (Vanquish, Thermo Fisher Scientific) with a Waters ACQUITY UPLC BEH Amide (2.1 mm × 50 mm, 1.7 μm) coupled to Orbitrap Exploris 120 mass spectrometer (Orbitrap MS, Thermo). The mobile phase consisted of 25 mmol/L ammonium acetate and 25 ammonia hydroxide in water (pH = 9.75) (A) and acetonitrile (B). The auto-sampler temperature was 4 °C, and the injection volume was 2 μl. The Orbitrap Exploris 120 mass spectrometer was used for its ability to acquire MS/MS spectra on information-dependent acquisition (IDA) mode in the control of the acquisition software (Xcalibur, Thermo). In this mode, the acquisition software continuously evaluates the full scan MS spectrum. The ESI source conditions were set as following: sheath gas flow rate as 50 Arb, Aux gas flow rate as 15 Arb, capillary temperature 320 °C, full MS resolution as 60,000, MS/MS resolution as 15,000, collision energy: SNCE 20/30/40, spray voltage as 3.8 kV (positive) or −3.4 kV (negative), respectively.

The raw data were converted to the mzXML format using ProteoWizard and processed with an in-house program. which was developed using R and based on XCMS, for feature detection, extraction, alignment, and integration. The R package and the BiotreeDB (V3.0) were applied in metabolite identification. In this study, features were detected and metabolites were left after relative standard deviation de-noising. Then, the missing values were filled up by half of the minimum value. Also, internal standard normalization method was employed in this data analysis. The final dataset containing the information of feature number, sample name and normalized feature area was imported to SIMCA18.0.1 software package (Sartorius Stedim Data Analytics AB, Umea, Sweden) for multivariate analysis. Data was scaled and logarithmic transformed to minimize the impact of both noise and high variance of the variables. After these transformations, PCA (principal component analysis, PCA), an unsupervised analysis that reduces the dimension of the data, was carried out to visualize the distribution and the grouping of the samples. 95% confidence interval in the PCA score plot was used as the threshold to identify potential outliers in the dataset. In order to visualize group separation and find significantly changed metabolites, supervised orthogonal projections to latent structures-discriminate. analysis (OPLS-DA) was applied. Then, a 7-fold cross validation was performed to calculate the value of R2 and Q2. R2 indicates how well the variation of a variable is explained and Q2 means how well a variable could be predicted. To check the robustness and predictive ability of the OPLS-DA model, 200 times permutations was further conducted. Afterward, the R2 and Q2 intercept values were obtained. Here, the intercept value of Q2 represents the robustness of the model, the risk of overfitting and the reliability of the model, which will be the smaller the better. Furthermore, the value of variable importance in the projection (VIP) of the first principal component in OPLS-DA analysis was obtained. It summarizes the contribution of each variable to the model. Since no metabolites remained statistically significant after multiple testing correction, and in keeping with the exploratory aim of the study, metabolites with a VIP score > 1 and an uncorrected P < 0.05 were selected as meaningful for further analysis. In addition, commercial databases including KEGG (http://www.genome.jp/kegg/↗) was used for pathway enrichment analysis.

For further exploration, correlation analyses were performed between clinical features and microbiome/metabolome results, as well as between microbiome and metabolome results themselves. Clinical measures were attack frequency [categorized as 1 (once every few months or years) or 2 (several times per month)] and DHI scores (total and subscale scores: DHI-P, DHI-F, DHI-E). Spearman correlation analysis and associated statistical graphs were performed/generated using GraphPad Prism 10.

The Venn diagram constructed based on ASVs () illustrates compositional differences in gut microbiota among the three groups. Compared with the Hc group, the VM group shared 929 ASVs and exhibited 909 unique ASVs, while the Hc group exhibited 2,129 unique ASVs. Compared with the M group, the VM group shared 903 ASVs and exhibited 935 unique ASVs, while the M group exhibited 2,101 unique ASVs. Between the M group and the Hc group, 1,158 ASVs were shared; the M group exhibited 1,846 unique ASVs, and the Hc group exhibited 1,900 unique ASVs. Supplementary Figure S1

There were no significant differences in the Chao1 index, observed species, Shannon index, Simpson index, or Pielou's evenness (pielou-e) index among the three groups, suggesting no marked differences in microbial community richness or diversity. However, the Chao1 index and observed species were significantly higher in the M group compared to the VM group, indicating greater microbial community richness in the M group (–,). Meanwhile, as shown in, Goods_coverage values approached 1 in all three groups, indicating that the sequencing results adequately reflected the true microbial composition of the samples. Supplementary Figures S2A E Supplementary Table S1 Supplementary Figure S2F

Principal Coordinates Analysis (PCoA) scatter plots revealed noticeable differences in the clustering distribution of microbial communities among the three groups. Similarity analysis (ANOSIM) further confirmed significant separation between the groups (M vs. VM vs. Hc, R = 0.4175, P = 0.001) (Supplementary Figure S3).

We focused on compositional changes at the phylum and genus levels. Analysis of the gut microbiota composition selected the top 30 most abundant taxa, and their relative abundance distribution across groups is presented as stacked bar charts (Figure 1). At the phylum level, Firmicutes, Actinobacteriota, Bacteroidota, Proteobacteria, and Verrucomicrobiota dominated the abundance in all three groups. Within the VM group, the relative abundance of Verrucomicrobiota demonstrated a decreasing trend, while Firmicutes and others showed an increasing trend, as detailed in Figure 1A. At the genus level, compared to the other two groups, genera such as Akkermansia exhibited a marked decreasing trend in relative abundance, whereas the Ruminococcus gnavus group and others showed a marked increasing trend (Figure 1B).

Following differential abundance analysis with FDR correction, significant reductions in the abundance of Verrucomicrobiota (q = 0.041) and Chloroflexi (q = 0.022) were observed in the VM group at the phylum level compared to the other two groups. Further pairwise comparisons revealed that the relative abundance of Chloroflexi was significantly lower in the VM group compared to the M group (q = 0.01), while Verrucomicrobiota showed a decreasing trend approaching significance (q = 0.055). Neither phylum exhibited a significant difference in abundance when compared to the Hc group (Supplementary Figure S4).

At the genus level, the relative abundances of all 44 differentially abundant genera among the three groups were below 0.5%. Among these, genera including Latilactobacillus, JG30-KF-AS9_unclassified, Parvimonas, Olsenella, Cellulomonas, Porphyromonas, and Peptostreptococcus exhibited higher relative abundance in the VM group. Further pairwise comparisons within these genera revealed that compared to the Hc group, Latilactobacillus was significantly elevated in the VM group (q = 0.017); compared to the M group, JG30-KF-AS9_unclassified was significantly elevated (q = 0.021), while Olsenella showed an increasing trend approaching significance (q = 0.0502). Notably, the relative abundance of multiple genera was zero (“0”) in one or two groups, which may reflect limitations in sequencing depth or associations with factors related to the disease's pathophysiology (Supplementary Figure S5).

Further LEfSe analysis revealed no taxa that remained statistically significant after correction for multiple comparisons. However, for exploratory purposes, we also report findings from the uncorrected analysis. At the genus level, Coprococcus was specifically enriched in the VM group (LDA = 3.91, P = 0.034) (Figure 2A). This enrichment was also significant when comparing the VM group to the Hc group (LDA = 3.87, P = 0.046) (Figure 2B). However, no significant enrichment of Coprococcus was observed in the VM group compared to the M group (LDA < 3, P > 0.05) (Figure 2C); instead, specific enrichment was found for the Ruminococcus gnavus group (LDA = 4.41, P = 0.020), Flavonifractor (LDA = 3.75, P = 0.046), and Parasutterella (LDA = 3.71, P = 0.027) in the VM group relative to the M group. However, no specific findings at the phylum level.

Following preprocessing, a total of 17,311 features were retained, of which 1,157 were secondarily identified metabolites. The classification of these secondarily identified metabolites is shown in. Further annotation and differential analysis revealed that compared to the Hc group, the VM group exhibited 54 up-regulated and 46 down-regulated metabolites; compared to the M group, the VM group exhibited 115 up-regulated and 102 down-regulated metabolites (). Supplementary Figure S6 Supplementary Table S2

OPLS-DA was employed to filter out orthogonal variables unrelated to the categorical variables and separately analyze the predictive (non-orthogonal) and orthogonal components. This approach provides more reliable information regarding the group differences in metabolites and their correlation with the experimental groups. The OPLS-DA score plot demonstrated clear separation between the three sample groups and good intra-group reproducibility (). Figure S7

Since no metabolites remained statistically significant after multiple testing correction among the three groups or between pairwise groups, the differential metabolites discussed below were screened based on an uncorrected P < 0.05 and VIP score > 1. It is important to note that these results carry a potential false positive risk; therefore, the subsequent analyses should be interpreted strictly as exploratory. Heatmaps were generated to visualize the top 20 differential metabolites identified in the comparisons among the three groups, between the VM and Hc groups, and between the VM and M groups (Supplementary Figure S8).

KEGG annotation and enrichment analysis of all differential metabolites among the three groups revealed that the pathways with the largest proportions were 2-Oxocarboxylic acid metabolism and Biosynthesis of amino acids (Supplementary Figure S9). Among the top 20 most significant differential metabolites associated with these pathways were pyruvate and 2-ketobutyric acid. Further pairwise comparisons showed that, compared to the Hc group (Figure 3A) and the M group (Figure 3B), 2-Oxocarboxylic acid metabolism remained the most prominent biological pathway. Similarly, pyruvate and 2-ketobutyric acid were consistently among the top 20 most significant differential metabolites involved in this pathway in both comparisons (Table 2).

Further KEGG topological analysis of differential metabolites among the three groups (Figure 4A) indicated that the Alanine, aspartate and glutamate metabolism pathway was significantly enriched (Rich factor = 0.177, P = 0.007). Pyruvate was among the top 20 most significant differential metabolites associated with this pathway (Supplementary Table S4). Pairwise comparisons revealed that analysis of differential metabolites between the VM and Hc groups (Figure 4B) also showed significant enrichment of the Alanine, aspartate and glutamate metabolism pathway (Rich factor = 0.177, P = 0.002), with pyruvate again being among the associated top 20 most significant metabolites (Supplementary Table S4). However, when compared to the M group (Figure 4C), this pathway was not significantly enriched (Rich factor = 0, P = 0.199). Instead, only the Ubiquinone and other terpenoid–quinone biosynthesis pathway showed a statistically significant difference (Rich factor = 0.038, P = 0.042), although its enrichment level was not prominent. Furthermore, this latter pathway did not include any of the top 20 most significant differential metabolites (Supplementary Table S4).

Following preprocessing, a total of 26,044 features were retained, of which 2,608 were secondarily identified metabolites. The classification of these secondarily identified metabolites is shown in. Further annotation and differential analysis revealed that compared to the Hc group, the VM group exhibited 11 up-regulated and 40 down-regulated metabolites; compared to the M group, the VM group exhibited 4 up-regulated and 161 down-regulated metabolites (). Supplementary Figure S10 Supplementary Table S3

OPLS-DA was employed to filter out orthogonal variables unrelated to the categorical variables and separately analyze the predictive (non-orthogonal) and orthogonal components. This approach provides more reliable information regarding the group differences in metabolites and their correlation with the experimental groups. The OPLS-DA score plot demonstrated clear separation between the three sample groups and good intra-group reproducibility (). Supplementary Figure S11

Same as the plasma results, since no metabolites remained statistically significant after multiple testing correction among the three groups or between pairwise groups, the differential metabolites discussed below were screened based on an uncorrected P < 0.05 and VIP score > 1. It is important to note that these results carry a potential false positive risk; therefore, the subsequent analyses should be interpreted strictly as exploratory. Heatmaps were generated to visualize the top 20 differential metabolites identified in the comparisons among the three groups, between the VM and Hc groups, and between the VM and M groups (Supplementary Figure S12).

KEGG annotation and enrichment analysis of all differential metabolites among the three groups identified Tyrosine metabolism and ABC transporters as the pathways with the largest proportions (). Among the top 20 most significant differential metabolites associated with these pathways, Norepinephrine featured prominently. Supplementary Figure S13

Further pairwise comparisons revealed that, compared to the Hc group (Figure 5A), the pathways with the largest proportions were Starch and sucrose metabolism, Galactose metabolism, Biosynthesis of amino acids, ABC transporters, and Taste transduction. Norepinephrine (NE) was among the top 20 most significant differential metabolites associated with these pathways. In contrast, when compared to the M group (Figure 5B), Tyrosine metabolism and Biosynthesis of cofactors were the most prominent pathways; however, none of their associated differential metabolites ranked within the top 20 most significant.

Further KEGG topological analysis of differential metabolites among the three groups (Figure 6A) indicated significant enrichment of the Tyrosine metabolism pathway (Rich factor = 0.153, P = 0.004), with NE featuring among the top 20 most significant associated metabolites. Analysis of differential metabolites between the VM and Hc groups yielded no statistically significantly enriched pathways. However, comparison with the M group (Figure 6B) showed higher enrichment levels for the Tyrosine metabolism pathway (Rich factor = 0.158, P = 0.003) and the Citrate cycle (TCA cycle) pathway (Rich factor = 0.100, P = 0.030). Neither pathway included any of the top 20 most significant differential metabolites.

To further explore the gut microbiota-disease link, we analyzed correlations between the relative abundance of specific genera and clinical measures in the VM group. The genera included were those with high relative abundance and significant intergroup differences (Latilactobacillus, JG30-KF-AS9_unclassified, Parvimonas, Olsenella, Cellulomonas, Porphyromonas, Peptostreptococcus) and Coprococcus, which was specifically enriched in the VM group per LEfSe analysis. Parvimonas correlated positively with attack frequency (r = 0.7307, P = 0.0250). JG30-KF-AS9_unclassified (r = 0.8845, P = 0.0011) and Olsenella (r = 0.6838, P = 0.0330) correlated positively with DHI-P scores. Peptostreptococcus correlated positively with DHI-F scores (r = 0.6831, P = 0.0350). JG30-KF-AS9_unclassified (r = 0.7964, P = 0.0084) and Porphyromonas (r = 0.7771, P = 0.0120) correlated positively with total DHI scores. No significant correlations were found between any genera and DHI-E scores (Supplementary Figure S14A–D).

We next analyzed correlations between plasma pyruvate, urinary NE, and clinical features (attack frequency, DHI scores) in the VM group. Plasma pyruvate showed a positive correlation with DHI-P scores (r = 0.7152, P = 0.0225; Supplementary Figure S15) but not with other clinical measures. Urinary NE showed no significant correlations with any clinical features analyzed.

Correlation analyses were performed between the relative abundance of the differential genera and the levels of plasma pyruvate and urinary NE. Plasma pyruvate correlated positively with the relative abundance of Cellulomonas (r = 0.6727, P = 0.0390; Supplementary Figure S16). No other significant correlations were observed between the metabolites and the other genera. Urinary NE showed no significant correlations with any of the differential genera.