Abstract

Background:

Vestibular migraine (VM) pathophysiology remains unclear, with research often extrapolating from migraine (M) studies. The microbiota-gut-brain axis represents a novel avenue for exploring VM mechanisms and treatments. This study aimed to compare gut microbiota, plasma and urine metabolome alterations among VM patients, M patients, and healthy controls (Hcs).

Methods:

A cross-sectional study recruited 15 VM patients, 15 M patients, and 15 age-/gender-matched Hcs (April–September 2024) from a tertiary hospital. Final analysis included 10 VM, 15 M, and 15 Hc participants. All underwent fecal 16S rRNA gene sequencing, plasma and urine metabolomics. Demographic and clinical data were collected.

Results:

Gut microbiota alpha/beta diversity showed no significant inter-group differences. At the phylum level, Verrucomicrobiota and Chloroflexi were reduced in VM vs. M. Genus-level analysis revealed trends (e.g., decreased Akkermansia, increased Ruminococcus gnavus group in VM) and LEfSe (Linear Discriminant Analysis Effect Size) identified several genera (e.g., Latilactobacillus, Parvimonas) with higher relative abundance in VM. Plasma metabolomics identified numerous differential metabolites; bioinformatics implicated amino acid metabolism pathways, with pyruvate as a key differential metabolite. Urine metabolomics and bioinformatics strongly associated differential metabolites with tyrosine metabolism and norepinephrine (NE).

Conclusions:

In this exploratory multi-omics study, VM patients showed distinct gut microbiota composition compared to M and HC, particularly at the genus level. Plasma metabolomics revealed alterations in pyruvate and amino acid metabolism pathways, suggesting possible energy metabolism disturbances in VM. Urine metabolomics highlighted the tyrosine metabolism pathway, with norepinephrine emerging as a metabolite of interest. These preliminary findings point to potential involvement of gut dysbiosis, metabolic perturbations, and neurotransmitter-related pathways in VM pathophysiology, providing a foundation for hypothesis-driven research. Given the exploratory nature of this study, these observations require validation in larger cohorts.

1 Introduction

VM is one of the most prevalent vestibular disorders, with an estimated 1-year prevalence of 2.7% in the general population (1) and 1.07% in Asians (2). Among neurology outpatients, the proportions of definite vestibular migraine (dVM) and probable vestibular migraine (pVM) are 4.6% and 1.1%, respectively (3). VM predominantly affects middle-aged women (4–7). Its core symptoms involve various forms of vertigo (vestibular symptoms), often accompanied by a history of migraine or migraine-like features during attacks, including visual aura, photophobia, or phonophobia (8). Moderate to severe vestibular symptoms significantly impair patients' daily life, work, and study. Reports indicate that VM patients have a higher rate of work absenteeism or school absence (60%) compared to patients with other vertigo or balance disorders (51%, P < 0.001). Furthermore, the proportion experiencing falls in the past 5 years (52.2%) is significantly higher than in the general population (17.0%, P < 0.001) or in patients with dizziness (43.7%, P < 0.01) (1). Additionally, Akdal et al. (9) identified cognitive impairment in VM patients, with another study (10) suggesting that the severity of dizziness may influence the degree of cognitive impairment.

The pathophysiological mechanisms underlying VM remain unclear, and no specific treatment guidelines exist. Despite increasing clinical research and the release of relevant expert consensus in China (11), current understanding and therapeutic approaches for VM are largely extrapolated from migraine. Proposed theories for VM pathogenesis include: a lowered trigeminovascular activation threshold; cortical spreading depression affecting vestibular processing centers in the brain; inherited ion channel defects; labyrinthine vasospasm; sterile inflammation affecting inner ear structures; and dysfunction in central vestibular pathways modulated by the dorsal raphe nucleus and locus coeruleus, potentially influencing cognitive behavior (12–18).

The growing recognition of the microbiota-gut-brain axis has spurred interest in the relationship between gut microbiota and neuropsychiatric disorders. Numerous studies have revealed gut dysbiosis, increased intestinal permeability, elevated pro-inflammatory cytokines, reduced levels of beneficial short-chain fatty acids (SCFAs) with anti-inflammatory properties, and dysregulated metabolism (e.g., tryptophan) in migraine patients. These alterations may contribute to migraine pathophysiology, potentially via neuroinflammation. Interventions like probiotics have shown promise in alleviating some migraine symptoms.

However, no studies to date have investigated the association between VM and gut microbiota. This observational study aims to compare the gut microbiota (via fecal 16S rRNA sequencing), plasma metabolome, and urine metabolome among VM patients, migraine (M) patients, and healthy controls (Hcs). We seek to identify differential gut bacteria and metabolites, and utilize bioinformatics to explore potential microbe-metabolite pathways. This foundational work aims to inform future research into the pathophysiology and therapeutic strategies for VM.

2 Materials and methods

This study was approved by the Ethics Committee of the First Affiliated Hospital of Harbin Medical University (Approval No. IRB-AF/SC-12/03.0). All participants provided written informed consent.

2.1 Participants

From April 2024 to September 2024, 15 patients with VM and 15 patients with migraine were recruited from the Dizziness and Vertigo Clinic, Headache Clinic, and Neurology Outpatient Department, as well as inpatients from Fourth Neurology Ward of the First Affiliated Hospital of Harbin Medical University. Additionally, 15 healthy controls were recruited from the family members of patients attending the aforementioned departments. The three groups were matched for age and sex. The inclusion Criteria are as follows: (1) Age ≥ 18 years; (2) Possession of a stable contact method and willingness to cooperate; (3) VM group meeting the diagnostic criteria for dVM or pVM as defined by the 2022 Bárány Society criteria, M group meeting the diagnostic criteria for migraine with aura or migraine without aura according to the International Classification of Headache Disorders, 3rd edition (ICHD-3) published by the International Headache Society (IHS) in 2018 and Hc Group should be absence of any current or past symptoms or history of headache or dizziness. Exclusion criteria applied to all groups included: (1) History of gastrointestinal diseases; (2) Use of antibiotics (including medications containing bismuth subsalicylate), acid-suppressing drugs, or microecological modulators (including probiotics/prebiotics/synbiotics) within 1 month prior to enrollment; (3) Comorbid psychiatric disorders, other chronic diseases (including cardiovascular diseases, hypertension, diabetes mellitus, malignancy, thyroid disorders, etc.), or any other major systemic illnesses; (4) Chronic use of oral medications. Both definite vestibular migraine (dVM) and probable vestibular migraine (pVM) patients meeting the diagnostic criteria were enrolled and are collectively referred to as the “VM group” throughout the manuscript.

This study was designed as an exploratory pilot investigation to generate hypotheses regarding multi-omic alterations in VM. Given the absence of prior multi-omic data in VM to inform effect size estimates, a formal a priori sample size calculation was not performed. We acknowledge that this sample size may limit statistical power; therefore, all findings, particularly those not surviving multiple testing correction, should be interpreted as preliminary and hypothesis-generating.

2.2 Data and sample collection2.2.1 Demographic and disease-related information

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.

2.2.2 Fecal sample collection and processing

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.

2.2.3 Plasma sample collection and processing

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.

2.2.4 Urine sample collection and processing

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.

2.3 Data processing and analysis2.3.1 Analysis of demographic and disease-related data

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.

2.3.2 Fecal 16S rRNA gene sequencing and data analysis

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).

2.3.3 Plasma and urine metabolomics and data analysis

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.

2.3.4 Correlation analyses of differential gut genera, metabolites, and clinical features

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.

3 Results3.1 Demographic data

Five participants in the VM group (designated as VM1, VM8, VM9, VM12, and VM14) withdrew from the study and revoked their informed consent. Consequently, the final VM group comprised 10 patients, all female, with a mean age of 41.5 ± 14.3 years. The M group included all 15 enrolled patients (1 male, 14 females), with a mean age of 30.9 ± 8.0 years. The Hc group had a mean age of 29.1 ± 6.3 years and included 2 males and 13 females. There were no significant differences in sex distribution (Fisher's exact test, P = 0.772) or mean age (P > 0.05) among the three groups (Table 1).

CharacteristicVM group (N = 10)M group (N = 15)Hc group (N = 15)F/PPAge41.5 ± 14.3 (31.3, 51.7)30.9 ± 8.0 (26.5, 35.4)29.1 ± 6.3 (25.6, 32.5)5.740.007Gender0.1590.772Male0 (0%)1 (6.7%)2 (13.3%)Female10 (100%)14 (93.3%)13 (86.7%)

Age and gender differences between groups.

One-way ANOVA showed a significant overall difference in age among the three groups (P = 0.007). However, due to unequal variances (Levene's test, P < 0.05), Tamhane's T2 post-hoc test was applied for pairwise comparisons, which revealed no significant differences between any two groups after correction for multiple comparisons: VM vs. M, P = 0.152; VM vs. Hc, P = 0.072; M vs. Hc, P = 0.862. Therefore, the three groups were considered comparable in terms of age distribution.

3.2 16S rRNA results3.2.1 ASVs analysis

The Venn diagram constructed based on ASVs (Supplementary Figure S1) 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.

3.2.2 Alpha diversity analysis

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 (Supplementary Figures S2A–E, Supplementary Table S1). Meanwhile, as shown in Supplementary Figure S2F, Goods_coverage values approached 1 in all three groups, indicating that the sequencing results adequately reflected the true microbial composition of the samples.

3.2.3 Beta diversity analysis

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).

3.2.4 Changes in gut microbiota compositional abundance

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).

Taxonomic classification of the top 30 species in terms of abundance at phylum (A) and genus levels (B).

3.2.5 Differential abundance analysis of gut microbiota composition

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.

The results of LEfSe difference analysis. (A) The results of LEfSe difference analysis across three groups (Hc, VM, M). (B) The results of LEfSe difference analysis between VM and Hc group. (C) The results of LEfSe difference analysis between VM and M group.

3.3 Plasma metabolomics results3.3.1 Raw data preprocessing results

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 Supplementary Figure S6. 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 Table S2).

3.3.2 Orthogonal partial least squares-discriminant analysis

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).

3.3.3 Screening of differential metabolites and kegg enrichment analysis

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).

3.3.4 KEGG annotation and enrichment analysis of differential metabolites

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).

KEGG annotation analysis of differential metabolites. (A) KEGG annotation analysis of differential metabolites between the VM and Hc groups. (B) KEGG annotation analysis of differential metabolites between the VM and M groups.

ComparisonDifferential metabolitesHc-M-VMPyruvate, 2-Oxobutyric acidHc-VMPyruvate, 2-Oxobutyric acid, L-glutamateM-VMPyruvate, 2-Oxobutyric acid, 2-Aminobutyric acid, Phenylalanine

Differential metabolites involved in the 2-Oxocarboxylic acid metabolism from KEGG annotation analysis of inter-group comparisons.

Bold entries belong to the top 20 significant differential metabolites mapped to the pathways in specific inter-group comparison.

3.3.5 KEGG topological analysis

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).

KEGG topological analysis of differential metabolites. (A) KEGG topological analysis of differential metabolites among the three groups. (B) KEGG Topological analysis of differential metabolites between the VM and Hc groups. (C) KEGG topological analysis of differential metabolites between the VM and M groups. The x-axis (Impact) represents the topological influence factor, with circle size proportional to pathway enrichment level. The y-axis (-ln P) displays the negative natural logarithm of P-values (base e). Color intensity denotes statistical significance, where -ln P > 3 (equivalent to P < 0.05) indicate significance thresholds, with darker hues reflecting higher confidence.

3.4 Urine metabolomics results3.4.1 Raw data preprocessing results

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 Supplementary Figure S10. 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 Table S3).

3.4.2 Orthogonal partial least squares-discriminant analysis

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).

3.4.3 Screening of differential metabolites and kegg enrichment analysis

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).

3.4.4 KEGG annotation and enrichment analysis of differential metabolites

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 (Supplementary Figure S13). Among the top 20 most significant differential metabolites associated with these pathways, Norepinephrine featured prominently.

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.

KEGG annotation analysis of differential metabolites. (A) KEGG annotation analysis of differential metabolites between the VM and Hc groups. (B) KEGG annotation analysis of differential metabolites between the VM and M groups.

3.4.5 KEGG topological analysis

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.

KEGG topological analysis of differential metabolites. (A) KEGG Topological Analysis of Differential Metabolites among the three groups. (B) KEGG Topological Analysis of Differential Metabolites between the VM and M groups. The x-axis (Impact) represents the topological influence factor, with circle size proportional to pathway enrichment level. The y-axis (-ln P) displays the negative natural logarithm of P-values (base e). Color