Gut microbial and metabolomic alterations associated with herpetic trigeminal neuralgia: an exploratory multi-omics study

2.1 Participant recruitment

This study recruited healthy controls (HC) without known diseases and patients diagnosed with HTN from the First Affiliated Hospital of Zhengzhou University (Henan, China) between October 2025 and January 2026. The clinical diagnosis of HTN was strictly formulated in accordance with the ICHD-3 diagnosis criteria. Fecal and blood samples were collected from all enrolled participants, and baseline demographic and clinical information, including age, sex, and body mass index (BMI), was systematically recorded. For HTN patients, clinical pain severity was evaluated via the Visual Analogue Scale (VAS), with scores documented on the day of sample collection. Additionally, disease duration and the distribution of affected trigeminal nerve branches were comprehensively recorded. A total of 37 participants were enrolled in this study, consisting of 18 HC and 19 HTN patients.

2.2 Inclusion and exclusion criteria

The inclusion criteria were as follows for HTN patients: (1) Consistent with the diagnostic guidelines of the ICHD-3 [9], enrolled patients were diagnosed with two subtypes of herpes zoster-related trigeminal neuralgia: (a) Painful trigeminal neuropathy attributed to herpes zoster (ICHD-3 13.1.2.1, Unilateral facial pain of less than three months’ duration but more than 1 month’s in the distribution(s) of one or more branches of the trigeminal nerve, caused by and associated with other symptoms and/or clinical signs of acute herpes zoster); (b) Trigeminal postherpetic neuralgia (ICHD-3 13.1.2.1): Unilateral facial pain persisting or recurring for at least three months in the distribution(s) of one or more branches of the trigeminal nerve, with variable sensory changes, caused by herpes zoster; (2) Aged between 34 and 77 years; (3) Moderate to severe pain with a VAS score ≥4; and (4) Disease duration ≥1 month.

The inclusion criteria for HC were: (1) Age-matched with HTN patients (34–77 years old) with consistent gender distribution; (2) No personal history of herpes zoster, trigeminal neuralgia, or any other neuropathic pain disorder; and (3) Free from severe systemic diseases, including diabetes mellitus, autoimmune disorders, malignant tumors, and active infectious diseases.

The exclusion criteria applicable to all participants were as follows: (1) A confirmed diagnosis of irritable bowel syndrome (IBS) or inflammatory bowel disease (IBD); (2) Clinical presentation with diarrhea symptoms; (3) A definitive diagnosis of any psychiatric disorders; and (4) A history of active smoking or alcohol dependence.

2.3 Sample collection and sequencing

Fecal samples were immediately frozen at −80 ℃ after collection and processed for metagenomic and metabolomic sequencing within one month. Whole blood samples were processed promptly to isolate serum, which was immediately aliquoted and preserved at −80 ℃ for subsequent metabolomic analysis. All metabolomic sequencing procedures were completed within one month following sample preservation to ensure sample stability and data reliability.

2.4 Metagenomic sequencing and untargeted metabolomics (LC–MS/MS) analysis

Metagenomic sequencing was performed by Shanghai Biotree Biomedical Technology Co., Ltd. All sequencing procedures were implemented on the Illumina NovaSeq 6000 platform (San Diego, CA, USA). Briefly, total microbial DNA was extracted from qualified fecal samples and subjected to rigorous quality control. Raw sequencing data were generated in FASTA with Quality scores (FASTQ) format, and the acquired reads were subsequently corrected and assembled into contigs. After dereplication, gene prediction was conducted to construct a nonredundant amino acid sequence dataset. Multiple public databases were adopted for the functional annotation of amino acid sequences, yielding functional abundance profiles at distinct hierarchical levels, which supported the subsequent analyses of microbial functional composition, diversity, and differential enrichment. Furthermore, high-quality datasets and nonredundant amino acid sequences were used for taxonomic annotation, thereby characterizing microbial community profiles at the species and subspecies levels.

Untargeted metabolomic analysis based on Liquid Chromatography-Tandem Mass Spectrometry (LC–MS/MS) was also performed by Shanghai Biotree Biomedical Technology Co., Ltd. Metabolite separation was implemented using an ultra-high-performance liquid chromatography system (Vanquish, Thermo Fisher Scientific, Waltham, MA, USA) equipped with a Waters ACQUITY Premier HSS T3 column (1.8 μm, 2.1 × 100 mm). The mobile phases consisted of aqueous phase A (ultrapure water with 0.1% formic acid) and organic phase B (acetonitrile containing 0.1% formic acid). For data preprocessing, metabolite features with missing values exceeding 50% were excluded from downstream analyses. Remaining features with a missing value proportion of ≤50% were imputed using the k-nearest neighbor (KNN) algorithm (k = 10), where missing values were estimated based on the weighted average of measurements from adjacent samples within the same group [10].

Metabolite annotation was accomplished by matching acquired MS/MS spectra against the proprietary database of Biotree Biotech Co., Ltd (Shanghai, China), with supplementary verification from public databases, including the Human Metabolome Database (HMDB) and Kyoto Encyclopedia of Genes and Genomes (KEGG). Metabolite identification confidence was categorized following the standardized criteria proposed by the Metabolomics Standards Initiative (MSI). Specifically, level 1 (Confident Identification) required consistent matching of both retention time and MS/MS fragment profiles with authentic reference standards. Level 2 (Putative Identification) was defined by spectral matching against public databases (e.g., HMDB, METLIN (A metabolite mass spectral database)). Level 3 (Tentative Annotation) relied solely on accurate mass matching. All metabolites identified in this study satisfied Level 2 identification standard [11].

2.5 Data analysis

For microbiome data analysis, alpha and beta diversity analyses were conducted using the vegan and picante packages in R. Phylogenetic trees were constructed and visualized with the ggtree, ape, ggstar, and ggplot2 packages. Linear discriminant analysis effect size (LEfSe) was performed using the microeco R package to identify differential microbial features between groups. Differential metabolite analysis was implemented using the limma R package. Linear models were fitted with lmFit, and intergroup contrast matrices were constructed. Heatmaps of differential metabolites were generated using ggplot2. Volcano plots were constructed with log₂ fold change (log₂FC) on the x-axis and −log₁₀ (p-value) on the y-axis, while p-values were corrected using the Benjamini-Hochberg (BH) method with false discovery rate (FDR) <0.05. For metagenomic sequencing data, the LEfSe algorithm was applied, and the Kruskal-Wallis test was used to identify differentially abundant bacterial taxa between groups (p < 0.05). Gene Ontology (GO) functional annotation of metagenomic sequences was performed using Blast2GO [12], with the non-redundant (NR) protein database as the reference. Genes were annotated into three core GO categories, including biological process (BP), molecular function (MF), and cellular component (CC), to establish gene-function correspondence [13].

The GO annotation results were then integrated with the transcripts per million (TPM) abundance matrix. For each GO term, TPM values of all assigned genes were summed to construct a GO term-level abundance matrix, thereby minimizing the influence of sequencing depth differences on intergroup comparisons. Metabolic KEGG enrichment analysis of metabolites was conducted using the clusterProfiler R package (version 4.16.0), with functional annotation reference to the KEGG database. The Wilcoxon rank-sum test was used to evaluate the differential alterations between the HC and HTN groups. Multiple testing correction was performed using the BH method. GO terms and KEGG pathways with adjusted p < 0.05 were defined as statistically significantly enriched.

Given the exploratory design and the limited sample size of this study, all machine-learning analyses, including least absolute shrinkage and selection operator (LASSO), random forest, and WGCNA, were performed for hypothesis generation rather than confirmatory validation. WGCNA was specifically applied to explore potential modular correlations among differential serum and fecal metabolites, and the sample size adopted in this study complied with the empirical reference range recommended by the official guideline [14]. Differential serum metabolites and the top 150 fecal metabolites ranked by variance were incorporated for network construction. The optimal soft-thresholding power was determined as 11 based on the scale-free topology criterion (Scale-Free Topology Model Fit R² ≥ 0.8), and a signed weighted co-expression network was subsequently established. Module partitioning was conducted using the blockwiseModules function, with the minimum module size set to 10 and the module merging threshold (mergeCutHeight) set to 0.25.

For LASSO regression analysis, all 37 enrolled participants were included to construct a binary diagnostic model for distinguishing HTN patients from healthy controls. The full dataset was randomly partitioned into a training set (n = 26) for model construction and an internal test set (n = 11) for preliminary validation.

Module eigengenes (MEs) were subsequently calculated and correlated with clinical phenotype, with HTN status encoded as a binary variable (0/1), using Pearson correlation analysis. Phenotype-related key modules were screened according to module–trait correlation coefficients (r) and corresponding p-values. For each key module, module membership (MM; representing the correlation of individual metabolites with MEs) and gene significance (GS; representing metabolites-phenotype correlation) were quantified. GS–MM scatter plots were generated to verify the biological validity of core modules. Furthermore, metabolites within individual modules were ranked by absolute GS values to screen core metabolites strongly correlated with the HTN phenotypes.

Random forest regression models were constructed for exploratory analysis to characterize potential associations between differential microorganisms and metabolites. The relative abundances of the 47 differential microbial taxa were defined as independent variables, while metabolite abundances were defined as dependent variables. Separate random forest models were built for each metabolite using 500 trees (ntree = 500) with a random seed set to 1234.

Model performance was evaluated using the coefficient of determination (R²), root mean square error (RMSE), and mean absolute error (MAE). Leave-one-out cross-validation (LOOCV) was performed to compute cross-validated mean squared error (CV-MSE), and the relative error was calculated as the ratio of CV-RMSE to the observed mean value for evaluating model generalizability. Metabolites models with R² > 0.5 and relative error <0.5 were regarded as statistically reliable. Feature importance scores (%IncMSE) were extracted to quantify the explanatory power of individual microbial taxa. The top 20 metabolites with the highest cumulative %IncMSE scores were selected for visualization, and stacked bar plots were generated to delineate the relative contributions of different microbial taxa to metabolite variation.

Spearman’s rank correlation analysis was conducted to explore the relationships between the relative abundances of microbial genera and VAS pain scores. Pain scores were defined as independent variables, and microbial relative abundances were set as dependent variables. Spearman correlation coefficients (ρ) and corresponding p-values were calculated for each microbial genus, with the BH method applied for multiple testing correction. Scatter plots with fitted linear regression trend lines and 95% confidence intervals (CI) were generated for visualization.

All R-based analyses and visualizations were implemented in RStudio (version 2025.05.1).

Baseline clinical data were analyzed using SPSS 30.0 (IBM, Armonk, NY, USA). Specifically, age and BMI differences between the two groups were compared using an independent-samples t-test, while gender distribution was assessed using the chi-square (χ²) test.

3.1 HTN is associated with the altered gut microbiota composition

Baseline demographic and clinical characteristics of all enrolled participants are presented in Table 1 . No significant intergroup differences in age, gender, or BMI were observed between the HTN and HC groups (all p > 0.05), indicating well-matched baseline conditions for subsequent comparative analysis.

Fecal metagenomic sequencing was performed to profile gut microbial compositional and functional alterations in HTN patients. Significant reductions in microbial αdiversity were detected in the HTN group, as reflected by decreased Gini-Simpson index and Shannon index (p < 0.05; Fig. 1a).

The Shannon index comprehensively reflects microbial community richness and evenness, therefore, its reduction indicated simultaneous declines in species richness and uniform distribution in the gut microbiota of HTN patients. The Gini-Simpson index primarily characterizes the diversity of dominant taxa; its decrease suggested simplified community structure and reduced diversity of dominant microbial populations in HTN patients. Although the Observed Species index did not reach statistical significance, a consistent declining trend in species richness was observed in the HTN group.

Principal coordinate analysis (PCoA) based on metagenomic data was further conducted to evaluate intergroup differences in β-diversity (Fig. 1b). The overall gut microbial community structure differed significantly between HTN patients and HC (p = 0.01), implying that HTN development is closely associated with global alterations in gut microbiota composition.

Partial overlap of confidence ellipses between the two groups in the two-dimensional PCoA plot, nevertheless, a clear intergroup separation trend was identified. HC samples exhibited tight clustering and were predominantly distributed on the left axis, indicating stable gut microbial composition with lower individual variability. In contrast, HTN samples were broadly dispersed across multiple quadrants, demonstrating increased interindividual heterogeneity and disrupted gut microbiota homeostasis.

LEfSe analysis was further applied to screen characteristic gut microbial taxa that distinguished the two groups. The cladogram results revealed that differential microbial taxa were predominantly enriched in the HC group (Fig. 1c). Using an linear discriminant analysis (LDA) score >3 and an adjusted p-value (padj) < 0.05 as screening criteria, a total of 89 characteristic microbial taxa were identified, and the top 30 ranked by LDA score were visualized (Fig. 1d).

Fig. 1.Gut microbial diversity and differential taxa between healthy controls (HC) and patients with herpetic trigeminal neuralgia (HTN). (a) Boxplot of microbial alpha diversity. *indicates significant differences between groups. (b) Microbial beta diversity between HC and HTN groups. (c) Phylogenetic tree of LEfSe differential analysis. (d) LDA score plot of LEfSe differential analysis. *indicates a statistically significant difference (*, p < 0.05); **indicates a highly statistically significant difference (**, p < 0.01). PCoA: Principal coordinate analysis; PERMANOVA: Permutational Multivariate Analysis of Variance; LDA: Linear Discriminant Analysis.

At the phylum level, Proteobacteria and Firmicutes were identified as the primary differential taxa between the two groups. Proteobacteria were markedly enriched in the HTN group, whereas Firmicutes predominated in HC group.

At the genus level, the HC group showed significant enrichment of Faecalibacterium, Blautia, Roseburia, Dialister, and Anaerobutyricum. In comparison, the HTN group exhibited increased abundance of Mediterraneibacter and Escherichia, which contain multiple opportunistic pathogenic species.

At the species level, significantly enriched species in the HC group included [Eubacterium] rectale, Oscillospiraceaebacterium, Faecalibacterium prausnitzii, Bifidobacterium adolescentis, Prevotella copri, and Blautia obeum. By contrast, [Ruminococcus] gnavus and Escherichia coli were significantly enriched in the HTN group.

GO functional enrichment analysis was further performed to explore microbial functional variations (Fig. 2). In the Biological Process category, multiple pH homeostasis-related pathways were significantly downregulated in the HTN group, including cellular pH reduction, intracellular pH reduction, and lithium ion response, suggesting impaired capacity to maintain intracellular acid-base balance. Moreover, both positive and negative regulatory pathways of endoribonuclease activity were significantly suppressed, indicating impaired microbial RNA degradation and regulatory functions. Conversely, upregulated pathways in HTN patients mainly involved transcriptional regulation (DNA-binding transcription factor activity), alcohol stress responses (ethanol and alcohol response), and quinolinate metabolism and biosynthesis.

Fig. 2.Functional differences of gut microbiota between HTN patients and HC. The Wilcoxon rank-sum test was performed on the abundance matrix obtained from GO functional annotation to evaluate functional differences for each GO term between groups. The bar chart displays the top 20 most significantly enriched pathways in each category. Blue represents downregulated functions in the HTN group, while red indicates upregulated functions in the HTN group. HTN: Herpetic trigeminal neuralgia; HC: healthy controls; UV: Ultraviolet; GTPase: guanosine triphosphatase; ATP: adenosine triphosphate; NAD: nicotinamide adenine dinucleotide; NADP: nicotinamide adenine dinucleotide phosphate.

In the Cellular Component category, the aminodeoxychorismate synthase complex, pellicle, and inner membrane pellicle complex were significantly enriched in the HTN group. For Molecular Function, ribonuclease inhibitor activity exhibited the most prominent downregulation. Additional suppressed functions included nitrite transmembrane transporter activity, oxidoreductase activity, and hydrogen dehydrogenase nicotinamide adenine dinucleotide phosphate, oxidized form (NADP⁺) activity. In contrast, upregulated molecular functions in HTN group comprised kinase regulator activity and protein kinase regulator activity, as well as multiple metabolic enzyme activities, including acetylesterase, short-chain carboxylesterase, mannitol-1-phosphate 5-dehydrogenase, mannosyltransferase, and glutamate-cysteine ligase activities. These functional alterations were primarily associated with oxidative stress defense and polysaccharide biosynthesis and metabolism.

Collectively, the gut microbiota of HTN patients exhibited systematic downregulation of core homeostasis functions, including pH stabilization, nitrite metabolism, redox balance, and ribonuclease inhibition. Meanwhile, compensatory upregulation was observed in pathways related to guanosine triphosphatase (GTPase)/kinase-mediated signaling, polysaccharide metabolism, bacterial membrane structure remodeling, and oxidative stress responses. These findings indicate that gut microbial dysfunction may participate in HTN pathological progression through multiple biological pathways, although causal inference cannot be drawn from this observational study design.

3.2 HTN is associated with altered fecal metabolite profiles in patients

Volcano plot analysis of fecal metabolites revealed extensive metabolic perturbations in HTN patients (Fig. 3a). Compared with HC, the HTN group presented a greater number of significantly upregulated metabolites with larger fold changes, indicating profound disruption of the intestinal metabolic microenvironment in HTN patients.

To screen robust metabolite biomarkers, differentially abundant fecal metabolites identified via the limma package were subjected to LASSO regression analysis for dimensionality reduction and feature selection. A panel of 14 candidate fecal metabolites were finally retained (Fig. 3b). LASSO regression with binomial deviance cross-validation was further applied to determine core predictive variables and establish a predictive model. The area under the receiver operating characteristic (ROC) curve (AUC) for the model in the test set was 0.833, indicating favorable discriminatory power (Supplementary Fig. 1a,b).

Fig. 3.Identification and evaluation of fecal differential metabolites. (a) Differential metabolites between the HC and HTN groups were systematically screened using the limma linear model based on fecal metabolomic data. Raw data were transformed via log₂(x+1) to satisfy the normality assumption. A design matrix without an intercept term was constructed, with the contrast set as HTN versus HC. Linear fitting was performed using the lmFit function, and intergroup comparisons were calculated via contrasts fit. Empirical Bayesian smoothing was implemented with the eBayes function to improve the stability of statistical testing. The Benjamini-Hochberg (BH) method was applied for false discovery rate (FDR) correction in multiple testing. Metabolites with (p < 0.05 and |log₂FC| > 1) were defined as significantly differentially abundant metabolites. (b) All samples were randomly divided into a training set (n = 26) and a test set (n = 11) at a ratio of 7:3. Binomial logistic regression (family = “binomial”) was adopted as the modeling framework. Metabolites with non-zero coefficients were extracted as candidate diagnostic features, and their expression patterns were visualized by heatmap. HTN: Herpetic trigeminal neuralgia; HC: healthy controls.

Notably, several candidate metabolites corresponded to clinically administered drugs, including Ganciclovir, Paracetamol, and Metformin. Given the study design, the accumulation of drug-related metabolites likely reflected individual medication exposure histories rather than inherent HTN-specific metabolic alterations. Accordingly, these compounds were interpreted as medication-related signatures and excluded from disease-related biomarker interpretation. Additionally, multiple xanthine derivatives, such as 2′,7-Dihydroxy-4′-methoxy-8-prenylflavan and 8-[(2-Hydroxyethyl) amino]-1,3,7-trimethylxanthine, were significantly enriched in the HTN patients, implying that aberrant purine pathway activation may contribute to HTN progression. In contrast, 9,10,13-Trihydroxystearic acid and 4-(4-Fluoroanilino)-2H-chromen-2-one were predominant in HC, and their depletion in HTN patients suggested impaired intestinal lipid metabolic homeostasis under disease conditions.

3.3 HTN is associated with altered serum metabolite profiles in patients

Untargeted metabolomic profiling was further performed on serum samples to characterize systematic metabolic alterations in HTN patients. In contrast to the fecal metabolome, the serum metabolome exhibited relatively mild perturbations, with fewer significantly differentially abundant metabolites identified between the two groups (Fig. 4a).

Among the differential serum metabolites identified (Fig. 4b), Metoprolol, 2,2′-(4-(2-Hydroxyethylamino)-3-nitrophenylazanediyl)diethanol, dimethylidenebutanedioylcarnitine, Trigonelline, 3-Pyridylacetic acid, and Aurintricarboxylic acid were significantly elevated in HTN patients, whereas 3,3,4,4,5,5,6,6,6-Nonafluorohexyl 2-methylacrylate, Trifluoroacetic acid, and Euparin were significantly decreased. Notably, Metoprolol—a widely prescribed β-adrenergic receptor blocker—was detected at elevated levels in patient serum. This enrichment was attributed to clinical medication exposure rather than inherent HTN-associated metabolic dysregulation, thereby representing a medication-related signature instead of a disease-specific biomarker. Of note, both Trigonelline and 3-Pyridylacetic acid are key intermediates involved in nicotinate metabolism. Additionally, Euparin, a naturally coumarin derivative with documented anti-inflammatory and neuroprotective properties, was significantly downregulated in HTN patients, suggesting its potential protective role in HTN pathogenesis.

Subsequent KEGG pathway enrichment analysis was performed based on all differential serum metabolites (p < 0.05). The results indicated significant enrichment of the nicotinate and nicotinamide metabolism pathway and the retrograde endocannabinoid signaling pathway in the HTN group (Fig. 4c), implying that dysregulation of these two metabolic signaling axes may contribute to HTN progression.

Fig. 4.Identification of serum differential metabolites. (a) The limma linear model was similarly applied to analyze differential metabolites between the HC and HTN groups. With the screening thresholds set at p < 0.05 and |log2FC| > 1, a total of nine characteristic differential metabolites were identified. (b) Heatmap of serum differential metabolites. (c) Metabolite pathway Sankey diagram-enrichment Bubble Plot. HTN: Herpetic trigeminal neuralgia; HC: healthy controls; PC: Phosphatidylcholine.

3.4 Correlations between differential gut microbiota and metabolites

To systematically explore the potential regulatory associations between differential gut microbial taxa and fecal metabolites in HTN patients, random forest regression models were constructed. Differential microbial species and genera were incorporated as independent variables to predict the abundance variation of candidate fecal metabolites, thereby elucidating microbe–metabolite interaction patterns (Fig. 5).

Fig. 5.Contribution of characteristic microorganisms to changes in fecal signature metabolites. The prediction model was constructed using the randomForest package, with the number of decision trees set to 500 (ntree = 500) and variable importance evaluation enabled (importance = TRUE). All target variables were ranked in descending order according to their total feature importance scores, and the top 20 variables with the highest contribution were screened to plot a stacked bar chart. The length and corresponding value of each bar represent the influence level of alterations in characteristic microorganisms on changes in fecal signature metabolites.

The model analysis identified 2,6-Xylidine as the fecal metabolite most strongly modulated by gut microbiota. Lachnospira, Clostridium innocuum, Lachnospira eligens, Simiaoa sunii, and Roseburia hominis were recognized as the potential microbial contributors influencing its abundance.

Among medication-related metabolites, Metformin exhibited prominent correlations with gut microbiota composition. Specifically, Haemophilus parainfluenzae, Haemophilus, Veillonella parvula, Eubacterium, and Erysipelatoclostridium showed strong interactive relationships with Metformin, indicating that these intestinal metabolites of Metformin may be significantly modulated by these microbial taxa.

Furthermore, 2-Hydroxy-2-(propan-2-yl) butanedioylcarnitine exhibited feature importance scores >7 for Enterocloster bolteae, Eggerthella, and Eggerthella lenta, suggesting substantial microbial contributions to the regulation of this carnitine-related metabolite. In addition, 9,10,13-Trihydroxystearic acid was strongly associated with Wujia and Wujia chipingensis (importance score >7), suggesting that the abundance of this hydroxylated fatty acid is influenced by these newly identified microbial taxa.

Collectively, the exploratory random forest regression analysis revealed correlations between differential gut microbiota and fecal metabolites in the HTN cohort, providing novel clues for understanding the microbe-mediated metabolic dysregulation underlying HTN pathological progression.

WGCNA was further conducted to integrate fecal and serum metabolite profiles, aiming to systematically characterize coordinated cross-compartment metabolic perturbation patterns underlying HTN pathogenesis and reveal coordinated metabolic alteration patterns between the two compartments in the context of HTN (Fig. 6).

Fig. 6.Weighted gene co-expression network analysis (WGCNA) was performed for integrated modular analysis of serum and fecal metabolites. A total of 200 metabolites were included, consisting of 50 serum differential metabolites and 150 fecal differential metabolites with the largest variance. The analyzed samples comprised the HC group (n = 18) and the HTN group (n = 19).

A total of four metabolite co-expression modules were identified in the constructed network, including the turquoise module (86 metabolites), blue module (51 metabolites), brown module (44 metabolites), and grey module (19 unclassified metabolites). Pearson correlation analysis between MEs and the HTN phenotype status was subsequently performed to evaluate the phenotypic relevance of each metabolic module.

The blue module (r = 0.73, p = 3.8 × 10⁻⁷) and the brown module (r = 0.73, p = 3.3 × 10⁻⁷) were both found to be significantly positively correlated with the HTN phenotype, indicating that the metabolites comprising these two modules were globally upregulated in HTN patients. In contrast, although the turquoise module contained the largest number of metabolites, its correlation with the HTN phenotype did not reach statistical significance (r = 0.14, p = 0.19). Furthermore, the distribution of GS-HTN values within the turquoise module were predominantly negative, suggesting that the majority of metabolites in this module were downregulated in both feces and serum samples from HTN patients. These observations imply that the progression of HTN is accompanied by a widespread suppression of systemic metabolic activity.

Based on the GS–MM scatter plots, metabolites with currently high MM and high GS values were prioritized as candidate key metabolites within each module.

In the blue module, the fecal-derived metabolites 3-{[(1Z)-4-(4-Hydroxyphenyl) phthalazin-1(2H)-ylidene] amino} phenol and 2-(4-Pyridinyl)-4-quinolinecarbohydrazide both exhibited MM values >0.8 and GS values >0.3, indicating strong intramodular connectivity and robust association with the HTN phenotype. These metabolites were therefore designated as potential core metabolites of the blue module (Supplementary Fig. 2a, Supplementary Table 1).

In the brown module, multiple fecal-derived metabolites, including 5,7-Dihydroxy-2-(4-methoxyphenyl)-6-(3-methylbut-2-enyl)-2,3-dihydrochromen-4-one, 1,9-Dideoxyforskolin, Homoorientin, N-Acetyl-β-alanine, 7-[2,6-dimethyl-8-(2-methylbutanoyloxy)]heptanoic acid, Benzoic acid, Benzyl 6-O-(methylbutanoyl)-β-D-glucopyranoside, Isogravacridonechlorine, and Austalide L, along with serum-derived metabolites including Glycodeoxycholic acid, Glycochenodeoxycholic acid, and Glycohyodeoxycholic acid, all exhibited MM values >0.5 and GS values >0.3. These results suggested that the above metabolites were strongly correlated with both the brown module and the HTN phenotype. Notably, the clustered enrichment of conjugated bile acid metabolites implied that the bile acid metabolic axis may play an important role in the pathological progression of HTN (Supplementary Fig. 2b, Supplementary Table 1).

Hub metabolites within each module were subsequently identified based on connectivity values derived from eigengene analysis (kME (connectivity to the module eigengene)). The hub metabolite of the blue module was 3-{[(1Z)-4-(4-Hydroxyphenyl)phthalazin-1(2H)-ylidene]amino}phenol (kME = 0.83), whereas that of the brown module was 5,7-Dihydroxy-2-(4-methoxyphenyl)-6-(3-methylbut-2-enyl)-2,3-dihydrochromen-4-one (kME = 0.81). In the turquoise module, Cauloside C was designated as the hub metabolite (kME = 0.77). All three hub metabolites exhibited kME values >0.76, indicating a high degree of consistency with the co-expression patterns of metabolites within their respective modules (Supplementary Table 1).

To further evaluate the stability of module classification, bootstrap analysis was conducted by repeatedly reconstructing the co-expression network 100 times, each time using a randomly selected 80% subset of the samples. Module preservation analysis yielded a stability score of 0.840 ± 0.000 within the current dataset, indicating that the identified module structures were reasonably reproducible and relatively stable. These findings thus provide a preliminary basis for subsequent functional analyses (Supplementary Fig. 2c).

3.5 Association analysis of characteristic microorganisms and metabolites with HTN disease-related clinical features

To validate the associations between the identified candidate features and clinical pain characteristics, the VAS scores were used as a quantitative indicator of pain severity (Fig. 7).

Fig. 7.Correlations between metabolite/microbial markers and pain scores. (a) Correlation fitting curves between each metabolite marker and pain score using individually matched patient samples, with p < 0.05 set as the threshold for statistical significance. (b) Correlation fitting curves between each microbial marker and pain score in patient samples. VAS: Visual Analogue Scale.

Most characteristic metabolites were found to be negatively correlated with VAS scores (Fig. 7a). Among the medication-related metabolites, Ganciclovir (r = −0.446) and 8-Cyclopentyl-1,3-dimethylxanthine (r = −0.546) both exhibited significant negative correlations with VAS scores, with the latter showing a stronger association. Additionally, 2′,7-Dihydroxy-4′-methoxy-8-prenylflavan-2′,7-diglucoside also demonstrated a relatively strong negative correlation (r = −0.54), suggesting that the intestinal abundance of this flavonoid compound may be closely associated with pain alleviation. Notably, 4-(4-Fluoroanilino)-2H-chromen-2-one was the only characteristic metabolite that showed a positive correlation with VAS scores (r = 0.475), implying that its abundance increased with greater pain severity.

Correlation analysis between characteristic fecal microorganisms and VAS scores further revealed distinct association patterns among different microbial genera (Fig. 7b). Haemophilus, Roseburia, Simiaoa, Lachnospira, Faecalibacterium, and Anaerobutyricum were found to be significantly positively correlated with VAS scores, suggesting that the abundances of these genera increased with higher pain scores and may reflect reactive or adaptive responses during disease progression. In contrast, Eggerthella, Erysipelatoclostridium, Enterocloster, and Mediterraneibacter were significantly negatively correlated with VAS scores, indicating that a decrease in these genera was associated with aggravated pain severity. These microorganisms may, therefore, be involved in relieving HTN-related pain, potentially through anti-inflammatory activity or neuromodulatory metabolic pathways.

Among all microbial taxa analyzed, Haemophilus, Eggerthella, and Erysipelatoclostridium exhibited the strongest correlations with VAS scores, suggesting their potential utility as candidate microbial features indicative of pain severity in HTN patients.