Male C57BL/6J mice, aged 6–8 weeks, were obtained from SiPeiFu Biotechnology Co., Ltd. (Beijing, China). Mice were housed in an animal incubator under a 12-h light/dark cycle with ad libitum access to food and water. To minimize variability in behavioral outcomes both within and between groups, mice underwent pre-training before the commencement of the experiment. Additionally, experimenters who conducted the behavioral tests were blinded to the treatment groups to ensure unbiased assessment.

Mice were anesthetized using vaporized isoflurane (2.5% for induction and 2.0% for maintenance; R510-22-10, RWD Life Science, Shenzhen, Guangdong, China) and subsequently positioned in a stereotaxic frame for precise surgical procedures. Collagenase IV (Coll IV; C5138, Sigma-Aldrich, St. Louis, MO, USA; 0.01 U/10 nL), dissolved in saline, was stereotactically injected into the right ventral posterior medial and ventral posterior lateral nuclei of the thalamus in mice designated for the model group (group B). The injection coordinates relative to the bregma were as follows: 0.82–2.30 mm in the anterior-posterior axis, 1.30–1.95 mm lateral to the midline, and 3.01–4.25 mm ventral to the skull surface. The injection was performed at a rate of 0.01 nL/sec. Mice in the control group (group A) received an equal volume of sterile saline. After injection, the glass micropipette was left in place for 10 min to ensure complete diffusion of Coll IV. The micropipette was carefully withdrawn, and the cranial opening was sealed with bone wax. The scalp was sutured using 4-0 surgical sutures.

Mice in the intervention group (group C) received intraperitoneal injection of the NAT10 inhibitor, Remodelin (5 mg/kg, S7641, Selleck, Shanghai, China), beginning on the first day post-stereotaxic injection. Treatment was administered twice a week for a total of six injections. The health status, physiological parameters, and pain-related behaviors of all mice were monitored at baseline (pre-injection) and on days 1, 3, 5, 7, 14, and 21 post-stereotaxic injection to validate the successful establishment of the model.

Pain behavior was assessed through thermal and cold stimulation tests, measuring paw withdrawal latency in mice using a hot and cold plate pain meter (Model 35100, Ugo Basile, Gemonio, Italy) [18]. During the tests, each mouse was placed in a separate resin glass box. For the thermal stimulation test, mice were positioned on a glass plate maintained at 26 $\pm$ 0.5 °C. A light beam was directed at the center of the hind paw, and the duration (in seconds) until the mouse displayed a rapid paw withdrawal or licking response was recorded as the paw withdrawal latency. To prevent tissue damage, a minimum interval of 5 min was maintained between each stimulation, and each mouse underwent a total of 5 trials. The average latency from these trials was calculated and used as the final result. Similarly, for the cold stimulation test, the mouse was placed on an aluminum plate maintained at 0 °C, and the duration until the first paw withdrawal, licking response, or jumping was recorded. Each cold stimulation trial was separated by a 10-minute interval, in which three repetitions were conducted for each mouse, and the average latency was calculated.

Upon completion of the pain behavior assessments, motor function evaluations were conducted, including the cliff avoidance response, forelimb grip strength, and righting reflex tests [19, 20]. For the cliff avoidance test, a mouse was placed at the edge of a table, with its forepaws gripping the edge and its hind limbs dangling. The ability of the mouse to promptly return its hind paws to the table was then recorded. In the forelimb grip strength test, the mouse was placed on a wire grid, and the ability of its forepaws to grasp the grid was assessed. The righting reflex test measured the time it took for a mouse to reorient itself when placed on its back on a flat surface. Mice that required more than 2 seconds to right themselves were excluded to minimize individual variability. Each trial was conducted for 30 seconds, with a normal righting response defined as occurring within 2 seconds.

Following the behavioral tests, thalamic brain tissue samples from mice in each group were prepared for paraffin sectioning. Anesthesia was induced with isoflurane (5% for induction and 2% for maintenance, R510-22-10, RWD Life Science), after which perfusion was performed using phosphate-buffered saline (PBS) (G5001, Servicebio Technology Co., Ltd., Wuhan, Hubei, China), followed by 4% paraformaldehyde (20210827, GUANG FU, Tianjin, China). The thalamus was rapidly dissected and stored at 4 °C. Tissues were fixed in 4% paraformaldehyde for 24 hours. After fixation, the samples were dehydrated in graded ethanol solutions (100092683, Sinopharm), cleared with xylene (10023418, Sinopharm, Shanghai, China), and then embedded in paraffin (69019961, Sinopharm). Sections of 20 µm thickness were cut for Nissl staining and immunohistochemistry (IHC). Following the experiment, the carcasses were then sealed in bags, frozen, and transferred to certified disposal facilities for safe and compliant disposal.

Paraffin-embedded tissue sections were baked at 64 °C for 1 h, followed by dewaxing in xylene and rehydration in graded ethanol. The sections were thrice rinsed with PBS for 5 minutes each. They were then stained with Nissl staining solution (G1430, Solarbio, Beijing, China) at 50–60 °C for 30–60 minutes, followed by differentiation in 95% ethanol. The sections were then dehydrated in absolute ethanol, cleared in xylene, and finally mounted with neutral gum (10004160, Sinopharm) for long-term preservation. Once dried, the stained thalamic tissue was examined under a microscope (Leica DM500, Leica Microsystems, Wetzlar, Germany), and the number of neurons with Nissl bodies was quantified using ImageJ software (V1.53k; NIH, Bethesda, MD, USA).

Following the same dewaxing and hydration steps as in Nissl staining, antigen retrieval was performed using citrate buffer (G1202, Servicebio). Endogenous peroxidase activity was blocked by incubating the sections in 3% H₂O₂ (H1009, Sigma-Aldrich) for 20 minutes at room temperature. After blocking with 5% bovine serum albumin (BSA; A1933, Sigma-Aldrich), sections were incubated overnight with the following primary antibodies: rabbit anti-NAT10 (DF12427, 1:50; Affinity, Shanghai, China), rabbit anti-xanthine dehydrogenase (XDH, DF8111, 1:50, Affinity), rabbit anti-nucleoside phosphorylase (NP, HA500062, 1:150; Huabio, Woburn, MA, USA), rabbit anti-hypoxanthine phosphoribosyltransferase 1 (HPRT1 15059-1-AP, 1:50; Proteintech, Rosemont, IL, USA), and rabbit anti-guanine deaminase (GDA, DF2701, 1:50, Affinity). The next day, sections were treated with a reaction-enhancing solution for 20 minutes, followed by incubation with a universal secondary antibody for another 20 minutes. Staining was visualized using the 3,3^′-diaminobenzidine (DAB) substrate (G1212, Servicebio), and sections were counterstained with hematoxylin (G1004, Servicebio) for at least 15 minutes to achieve a blue background. The sections were dehydrated and cleared following the same protocol used for Nissl staining, with three 5-minute rinses in PBS between each step. Mean optical density was analyzed using ImageJ software.

Protein extraction from thalamic homogenates was performed using a commercial whole-protein extraction kit (EX1100, Solarbio). After quantifying protein concentration, samples were heated at 99 °C for 5 minutes. Subsequently, 30 µg of protein per sample was loaded onto 4–20% precast polyacrylamide gel electrophoresis (PAGE) gels (M42015, GenScript, Nanjing, Jiangsu, China) for electrophoresis. Proteins were then transferred onto polyvinylidene fluoride (PVDF) membranes (1620177, Bio-Rad, Hercules, CA, USA) via wet transfer. The membranes were blocked for 1 hour with 5% skim milk in Tris-buffered saline containing 0.1% Tween-20 (TBST). Primary antibodies, diluted in blocking buffer, were incubated overnight at 4 °C, including rabbit anti-XDH (DF8111, 1:1000; Affinity), anti-NP (HA500062, 1:1000; Huabio), anti-HPRT1 (15059-1-AP, 1:2000; Proteintech), anti-GDA (DF2701, 1:1000; Affinity), and anti-beta actin (AF7018, 1:10,000; Affinity). After washing with TBST, membranes were incubated for 2 hours at room temperature with a goat anti-rabbit IgG secondary antibody (S0001, 1:10,000; Affinity). Protein bands were visualized using ECL (1705061, Bio-Rad) and quantified with ImageJ software.

Protein extracts from mouse thalamic tissues collected 21 days post-modeling were analyzed for inflammatory cytokine levels using commercial ELISA kits (EM0005 for IL-10, EM0004 for IL-6, and EM0029 for IL-1$\beta$hua, Huabio) following the manufacturer’s instructions. Samples were diluted 5-fold before analysis. Briefly, 50 µL of standards or diluted samples were added to the wells, followed by 100 µL of enzyme-conjugated detection antibody, and incubated at 37 °C for 60 minutes. After five automated wash cycles, 50 µL of each of substrates A and B was added, followed by 15-minute incubation at 37 °C in the dark. The reaction was stopped with 50 µL of stop solution, and absorbance was measured at 450 nm. Cytokine concentrations were calculated using a four-parameter logistic standard curve.

The conserved regions of the 16S ribosomal deoxyribonucleic acid (16S rDNA)/Internal Transcribed Spacer 2 (ITS2) sequences were targeted for primer design. Total microbial DNA was extracted from mouse feces and used as a template for one-step polymerase chain reaction (PCR) with Phusion enzyme (F531L, Thermo Fisher Scientific, Waltham, MA, USA). Amplification was performed for 35 cycles. Sequencing adaptors and sample-specific barcodes were added to the amplified products via PCR. The PCR products were analyzed using 1.5% agarose gel electrophoresis, and the target fragments were extracted. The AxyPrep MAG PCR Clean-up kit (MAG- PCR-CL-50, Axygen Scientific, Union City, CA, USA) was utilized for the recovery. Purified PCR products were quantified and pooled using the Promega QuantiFluor fluorescence system with the Quant-iT PicoGreen dsDNA Assay kit (E2670, Promega Corporation, Madison, WI, USA). Finally, Illumina sequencing (PE250) was performed according to standard procedures.

Total RNA was isolated from thalamic tissue using the TRIzol reagent (15596026, Invitrogen, Carlsbad, CA, USA) following the manufacturer’s protocol. RNA concentration and purity were evaluated using a NanoDrop ND-1000 spectrophotometer (NanoDrop, Wilmington, DE, USA). The RNA integrity was assessed via a 2100 Bioanalyzer (Agilent, Santa Clara, CA, USA) and verified by agarose electrophoresis. RNA samples with a concentration greater than 50 ng/µL, a RIN above 7.0, an OD260/280 ratio exceeding 1.8, and a total RNA amount over 1 µg were selected for subsequent analysis. Oligo(dT) magnetic beads were utilized to capture polyadenylated mRNA, which was then fragmented under high temperatures. First-strand cDNA was synthesized with reverse transcriptase, followed by second-strand synthesis via RNase H digestion. The resulting double-stranded DNA was then end-repaired, followed by A-tailing and ligation with T-tailed adaptors. The resulting fragments (300 bp $\pm$ 50 bp) were purified and size-selected with magnetic beads using the Hieff NGS Ultima Dual-mode mRNA Library Prep kit for Illumina (12608ES, Yeasen, Ltd., Shanghai, China). Sequencing was conducted on an Illumina NovaSeq 6000 sequencing system (Illumina, San Diego, CA, USA) the paired-end (PE150) mode.

Blood samples (20 µL) were collected from each group of mice, and metabolites were extracted using a 50% methanol buffer (A452-4, Thermo Fisher Scientific). Briefly, 120 µL of precooled 50% methanol was added to 20 µL of the sample, followed by 1-minute vortexing and 10-minute incubation at room temperature. The samples were stored overnight at –20 °C. After centrifugation at 4000 g for 20 min, the supernatants were transferred to a 96-well plate. Quality control (QC) samples were prepared by mixing an equal volume (10 µL) of diluent from each sample. Samples were stored at –80 °C before liquid chromatography–mass spectrometry (LC-MS) analysis.

Chromatographic separations were conducted using an ultra-performance liquid chromatography (UPLC) system (SCIEX, Macclesfield, UK) with an ACQUITY UPLC T3 column (100 mm $\times$ 2.1 mm, 1.8 µm; Waters, Cheshire, UK). The column oven was maintained at 35 °C, with a flow rate of 0.4 mL/min. The mobile phase included solvent A (water with 0.1% formic acid; A117-50, Thermo Fisher Scientific) and solvent B (acetonitrile, A998-4, Thermo Fisher Scientific) with 0.1% formic acid. The gradient elution profile was as follows: 0–0.5 min, 5% B; 0.5–7 min, 5% to 100% B; 7–8 min, 100% B; 8–8.1 min, 100% to 5% B; 8.1–10 min, 5% B. The injection volume was 4 µL per sample.

Mass spectrometry analysis was conducted using a TripleTOF 5600+ system (SCIEX, Framingham, MA, USA) in both positive and negative ion modes. The ion source pressure was set to 30 Pounds per square inch (PSI), while gas 1 and gas 2 pressures were maintained at 60 PSI. The source temperature was maintained at 650 °C, and the ionization voltage was +5000 V for positive mode and –4500 V for negative mode. Data acquisition was performed in Information Dependent Acquisition (IDA) mode, with accuracy corrections implemented for every 20 samples. QC samples were analyzed after every 10 samples to monitor and mitigate batch effects.

To enhance the reliability of the 16S rRNA sequencing results, data processing included error correction, de-duplication, and chimera filtering to improve the accuracy of identifying amplicon sequence variants (ASVs). A Venn diagram was plotted to illustrate the common and unique ASVs among different sample groups, providing a visual comparison of similarities and specificities across groups. Differences in species diversity between samples were depicted using rarefaction curves, which allowed for the visualization of species richness and saturation.

Alpha and beta diversity analyses were conducted to evaluate the diversity within and between the samples, respectively. Alpha diversity, which reflects species richness and evenness, was assessed using the Chao1, Observed, Shannon, and Simpson indices. The Chao1 and Observed indices provide estimates of the number of species present in the community, while the Shannon index, based on information entropy, quantifies species diversity, and higher values indicate greater diversity. The Simpson index ranges from 0 to 1, and values approaching 1 indicate higher diversity. The differences in group-wise diversity were statistically tested using the Kruskal-Wallis test (p 0.05). Beta diversity was assessed using the principal coordinate analysis (PCoA) based on the weighted UniFrac distance to examine the differences in microbial community structure among the groups.

According to the species richness and taxonomic annotation, the relative abundance of the top 30 taxa at the genus level was compared across the groups. To further explore the differences and similarities between samples, the Bray-Curtis distance was used for hierarchical clustering. This clustering method enabled differentiation between taxa of high and low abundance, providing insights into the similarities and variations in the genus-level composition among the three groups, which were visualized using color gradients and similarity measures.

Subsequently, linear discriminant analysis (LDA) combined with linear discriminant analysis effect size (LEfSe) was applied to compare microbial compositions between groups (p 0.05 and $|$LDA score$|$ $\ge$3). The Kruskal-Wallis test was also employed to compare species differences among the three groups at the genus level, enabling the identification of differentially abundant taxa for further analysis (p 0.05).

The principal component analysis (PCA) was employed to investigate differences in metabolomic profiles across the three groups. Additionally, orthogonal partial least squares discriminant analysis (OPLS-DA) was utilized for pairwise comparisons (group A vs. group B and group B vs. group C) using both positive and negative ion modes. The variable importance in projection (VIP) values for each metabolite were calculated. The t-test was employed to analyze group 1 (model vs. control) and group 2 (intervention vs. model), identifying metabolites with p 0.05 and VIP $>$1 as differentially expressed metabolites. To identify shared differential metabolites affected by inhibitors, the intersection of upregulated and downregulated metabolites from group 1 was compared with downregulated and upregulated metabolites from group 2, respectively. Metabolites from both positive and negative ion modes were intersected, combined, and deduplicated to define the shared differentially expressed metabolites. Finally, the shared differentially expressed metabolites were mapped to relevant signaling pathways using the keggList function from the KEGGREST package (v1.38.0; Bioconductor, Fred Hutchinson Cancer Research Center, Seattle, WA, USA) [21] to identify associated Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways and their involvement in specific biological processes.

Differential gene expression between group 1 (model vs. control) and group 2 (intervention vs. model) was analyzed using the DESeq2 R package (v1.38.0; Bioconductor) [22]. The criteria for identifying differentially expressed genes (DEGs) were p 0.05 and $|$Log₂ Fold-Change (FC)$|$ $>$0.5. Shared DEGs were identified by intersecting the upregulated and downregulated genes from group 1 with the downregulated and upregulated genes from group 2, respectively. These shared DEGs were then subjected to the Gene Ontology (GO) and KEGG pathway enrichment analyses using the ClusterProfiler R package (ver. 4.7.1; Bioconductor) [23], with the statistical significance defined as p_adjust 0.05, to explore potential biological functions and pathways.

Multi-omics association analysis enables the identification of potential association patterns across different omics layers, facilitating the systematic elucidation of regulatory networks involved in complex biological processes and the discovery of key regulatory nodes. In this study, we integrated three omics layers: microbiome (differentially abundant microbial taxa), metabolome (shared differentially expressed metabolites), and transcriptome (shared DEGs). This study concentrated on molecules that were significantly altered across multiple layers, based on the hypothesis that coordinated changes may reflect functionally relevant regulatory crosstalk.

To investigate these associations, Pearson and Spearman correlation analyses were performed ($|$R$|$ $>$0.4, p 0.05) among these molecules. Significant correlations involving differential flora, metabolites, and genes were further analyzed using least absolute shrinkage and selection operator (LASSO) penalized Cox regression to identify key microorganism–metabolite and metabolite–gene pairs. A combined interaction network of hub microorganisms, hub metabolites, and hub genes was then constructed using the Cytoscape R package (ver. 3.8.1; The Cytoscape Consortium, San Diego, CA, USA; https://cytoscape.org) [24]. Additionally, Pearson and Mantel’s tests were conducted to identify metabolites significantly correlated with hub genes or hub microorganisms ($|$R$|$ $>$0.4).

Finally, a metabolite–response–enzyme–gene regulatory network was generated using the MetScape (https://apps.cytoscape.org/apps/metscape) plugin in Cytoscape to map biochemical pathways and interactions. Differential expression of hub genes, hub metabolites, and hub microorganisms was analyzed across groups A, B, and C to further investigate their functional roles.

To investigate the influences of xanthine on the thalamus, three groups were established: control group (group D), model group (group E), and febuxostat (FBX) intervention group (group F). Mice in the FBX group received 10 mg/kg of febuxostat (HY-13418, MedChemExpress, Shanghai, China) via intraperitoneal injection once daily for five consecutive days, beginning on the first day post-stereotaxic injection [24]. The treatment protocols for groups D and E followed the same procedures as those for groups A and B, respectively.

Thalamus tissue RNA from groups D, E, and F was collected to assess the expression levels of key genes. RNA extraction was conducted following the previously described procedure. Subsequently, 1 µg of total RNA was reverse transcribed into cDNA using the SureScript First-Strand cDNA Synthesis kit (G3130, Servicebio). Primers for key genes were synthesized by Tsingke Biotech Co., Ltd. (Beijing, China), and further details are presented in Supplementary Table 1. Quantitative reverse transcription PCR (RT-qPCR) was performed on the CFX96™ Real-Time PCR Detection system (Bio-Rad) using 2$\times$ Universal Blue SYBR Green qPCR Master Mix (G3320, Servicebio). Relative expression levels were calculated using the 2^{-Δ⁢Δ⁢CT} method [25], and GAPDH served as the housekeeping gene to normalize cycle threshold (CT) values.

RNA extraction was performed according to previously described protocols. The purified RNA samples were applied to Hybond-N+ membranes (HNB001, Solarbio), briefly air-dried, and subjected to UV crosslinking at 254 nm for 60 s. Membranes were then blocked with 5% non-fat milk for 1.5 h and incubated overnight at 4 °C with rabbit anti-ac4C antibody (ab252215, 1:250; Abcam, Cambridge, UK). After washing, membranes were incubated for 1.5 h with horseradish peroxidase (HRP)-conjugated goat anti-rabbit IgG antibody (ab6721, 1:10,000; Abcam). Signals were detected using an enhanced chemiluminescence kit (34580, Thermo Fisher Scientific).

The statistical analysis was conducted using R (v4.2.2; R Foundation for Statistical Computing, Vienna, Austria). Statistical significance was determined using the Kruskal-Wallis and Wilcoxon tests, and p 0.05 was considered statistically significant.

Compared with the control group, mice subjected to TH experienced a significant halt in body weight gain, as well as a reduction in food and water intake. However, administration of a NAT10 inhibitor partially reversed these effects (Fig. 1A–C). Both heat and cold pain hypersensitivity were observed in mice with TH. Remodeling treatment significantly reduced hypersensitivity induced by Coll IV microinjection. After the intervention, the latency period for response to cold and heat stimuli was significantly prolonged in group C, approaching a period similar to that observed in group A (Fig. 1D,E). Additionally, results from the grip strength test, plane righting reflex, and cliff avoidance test indicated that in group C, grip strength and righting reflex scores fell between those of groups A and B (Fig. 1F–H). Compared with group A, group B showed a significant increase in ac4C-modified mRNA level, whereas group C exhibited a reduction in ac4C level relative to group B (Fig. 1I). These findings suggest that inhibition of NAT10 has a therapeutic effect on thalamic hemorrhage.

Fig. 1.

Effect of the NAT10 inhibitor in the TH model. (A) Mouse weight measurements. (B) Monitoring of water intake. (C) Monitoring of feed intake in mice. (D) Response to heat pain stimulation. (E) Response to cold pain stimulation. (F) Forelimb grip strength assessment. (G) Plane righting test results. (H) Cliff avoidance test outcomes on day 21 post-intervention. (I) Analysis of ac4C expression level by dot blot in the mouse brain on day 21 post-intervention. (J) Nissl body observation and (K) NAT10 IHC staining in the mouse brain on day 21 post-intervention. (L) Visualization of bleeding points in the mouse brain on day 3 post-intervention. (J,K) Scale bars: 1000 μm in the upper panel and 100 μm in the lower panel, respectively. (L) Scale bar: 1 mm (each division). For panels (A–H), n = 6; for panel (K), n = 9; for panels (I) and (J), n = 3. Statistical analysis: (A–E) Repeated ANOVA with Bonferroni post-hoc test; (F–K) One-way ANOVA with Bonferroni post-hoc test. *p 0.05, **p 0.01, ***p 0.001, and ****p 0.0001 compared with group A. #p 0.05, ##p 0.01, ###p 0.001, and ####p 0.0001 compared with group B. In each panel, the small letter A indicates the control group (blue), B indicates the TH model group (red), and C indicates the NAT10 inhibitor intervention group (green). NAT10, N-acetyltransferase 10; TH, thalamic hemorrhage; ac4C, N4-acetylcytidine; IHC, immunohistochemistry. Group A, control group; group B, model group; group C, intervention group.

On day 3 post-intervention, thalamic tissues from all three groups were examined. Visible hemorrhages were observed in both group B and group C, with more severe bleeding in group B (Fig. 1L). On day 21, neuronal changes in the thalamus were assessed using Nissl staining. Neurons in group A maintained normal morphology, whereas group B exhibited extensive areas of infarction. Notably, the infarct area in group C was significantly reduced compared with that in group B (Fig. 1J). IHC revealed that NAT10 protein expression was upregulated in group B, while significantly suppressed in group C (Fig. 1K). These histological findings further support the therapeutic potential of NAT10 inhibitors in treating TH.

The multi-omics association analysis was carried out to gain a comprehensive understanding of microbial and molecular interactions across the different groups. Firstly, the gut microbiota was analyzed using a Venn diagram to assess the commonality and uniqueness of microbial species. A total of 788 core ASVs were shared among all three groups, while groups A, B, and C had 1534, 1819, and 1353 unique ASVs, respectively (Fig. 2A).

Fig. 2.

Analysis of gut microbiota. (A) The Venn diagram showed the distribution of ASVs among three groups of samples in the experiment, and clearly presented the shared and unique microbial taxa. (B) Assessment of microbial richness across all samples, quantifying the diversity of gut microbiota. (C) PCoA depicting sample clustering based on microbial community composition, revealing relationships among the samples. (D) The heatmap of microbial community composition (at the genus level) showed the relative abundance of major microbial genera among the three groups of samples. The color gradient reflected the level of abundance, and clearly presented the differences in intestinal microbial community composition among different treatment groups. (E) Identification of taxa with significant differences in richness between groups, highlighting those that were differentially enriched. (F) Overview of differentially abundant genera, highlighting those with significant variations in abundance across samples. (G) Summary of the top 30 differential microbial taxa, which were ranked by their abundance, showcasing the most significant alterations in gut flora. Group A: control group, Group B: TH model group, Group C: NAT10 inhibitor intervention group. ASVs, amplicon sequence variants; PCoA, principal coordinate analysis; LDA, linear discriminant analysis.

Next, species richness was evaluated both between and within the groups. As illustrated in Fig. 2B, there were notable variations in richness across the three groups. However, Alpha diversity indices, Chao1, Observed, Shannon, and Simpson, exhibited no significant differences in microbial diversity among the three groups (Supplementary Fig. 1). In contrast, PCoA revealed that samples from the three groups clustered separately, indicating subtle beta diversity differences (Fig. 2C).

The microbial community composition was examined at the genus level to identify specific differences among groups A, B, and C (Fig. 2D). The stacked bar plot revealed notable variations in the relative abundance rates of different genera across the three groups. Among them, Muribaculaceae_unclassified consistently exhibited a relatively high abundance in all the groups.

Cluster analysis provided additional insights into the community composition, displaying clear distinctions among the three groups. Notably, while group A displayed a unique microbial profile, groups B and C demonstrated greater similarity in their species composition. This suggests that the microbial communities in groups B and C may share more overlapping characteristics, which may be attributed to the shared treatment conditions.

Following the analysis of microbial diversity and genus-level composition across the three groups, LEfSe was employed to identify taxa with significant differential abundance (Fig. 2E). The LEfSe revealed 27 distinct taxa with differential expression among the three groups. Specifically, group A exhibited 11 distinct bacterial taxa, including p_Proteobacteria, g_Escherichia/Shigella, and s_Escherichia/Shigella_unclassified. Group B was characterized by g_Rikenellaceae_RC9_gut_group, s_Rikenellaceae_RC9_gut_group_unclassified, and g_Klebsiella. In contrast, group C had 13 significantly enriched genera, such as s_Clostridia_UCG_014_unclassified, o_Clostridia_UCG_014, and f_Clostridia_UCG_014_unclassified (Fig. 2F). Additionally, genus-level comparisons across the three groups identified a total of 40 bacterial taxa with differential expression, and the top 30 taxa were visualized in a bar chart (Fig. 2G). These findings highlight distinct microbial profiles and potential biomarkers across the three groups, reflecting how treatment or intervention may influence gut microbial composition.

To further explore the metabolic variations between groups, PCA was conducted to assess the variations in metabolomic profiles. The PCA results, based on principal components PC1 and PC2, revealed a low degree of separation among the three groups, suggesting subtle differences in their metabolomes (Supplementary Fig. 2A,B). However, OPLS-DA, performed under both positive and negative ion modes, successfully distinguished among the three groups, highlighting distinct metabolic differences (Supplementary Fig. 2C–F).

Next, to identify metabolites shared across different experimental conditions, the differentially expressed metabolites from group 1 (model vs. control) and group 2 (intervention vs. model) were intersected. This analysis yielded 20 shared differentially expressed metabolites in the positive ion mode and 22 in the negative ion mode, which were subsequently combined and refined to identify a total of 42 shared differentially expressed metabolites (Fig. 3A–D).

Fig. 3.

Analysis of metabolites in thalamic tissue. (A–D) Venn diagrams illustrating the identification of differentially expressed metabolites across the three groups. (E,F) Heatmaps represent the patterns of differentially expressed metabolites. Positive ionization patterns are illustrated in (A,B,E), while negative ionization patterns are presented in (C,D,F). Group designations are as follows: m, model; c, control; i, intervention. Group A, control group; group B, model group; Group C, intervention group.

Heatmaps representing the differentially expressed metabolites for both ionization modes were generated, exhibiting clear differences in metabolic profiles (Fig. 3E,F). In addition, the KEGG pathway analysis of these shared differentially expressed metabolites identified 29 enriched pathways, including arginine and proline metabolism, alpha-linolenic acid metabolism, and glycerophospholipid metabolism (Supplementary Table 2). These pathways provide insight into the potential biological mechanisms underlying the observed metabolic shifts.

Differential expression analysis revealed 499 DEGs in group 1 (model vs. control), comprising 315 upregulated and 184 downregulated DEGs (Fig. 4A,B). Similarly, group 2 (intervention vs. model) included a total of 664 DEGs, with 345 upregulated and 319 downregulated DEGs (Fig. 4C,D). By intersecting the upregulated and downregulated genes from group 1 with the corresponding downregulated and upregulated genes from group 2, 46 shared DEGs affected by the inhibitors were identified (Fig. 4E,F). The shared DEGs, metabolites, and bacterial genera identified above were subsequently compiled in Supplementary Table 3.

Fig. 4.

Transcriptome analysis of thalamic tissue. (A,C) Volcano plots displaying DEGs for (A) group 1 (TH model vs. control) and (C) group 2 (intervention vs. TH model). (B,D) Heatmaps illustrating the expression profiles of DEGs in (B) group 1 and (D) group 2. (E,F) Venn diagrams showing the shared DEGs between (E) group 1 and (F) group 2. (G) GO analysis and (H) KEGG pathway analysis of all the shared DEGs, providing insights into their functional roles and associated pathways. Different colors correspond to distinct pathways, and the size of each dot reflects the number of genes enriched in that pathway. DEGs, differentially expressed genes; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes. Control group (group A, c); model group (group B, m); intervention group (group C, i).

Functional and pathway enrichment analyses were subsequently conducted to investigate the potential roles of these shared DEGs. A total of 66 GO biological processes were enriched, including the neuropeptide signaling pathway, regulation of tumor necrosis factor production, and positive regulation of responses to biotic stimuli. In terms of cellular components, the shared DEGs were mainly enriched in the serine-type peptidase complex, serine-type endopeptidase complex, and Golgi lumen. Molecular function analysis indicated involvement in neuropeptide receptor binding, hormone activity, and phospholipase A1 activity (Fig. 4G). Additionally, the shared DEGs were enriched in eight KEGG signaling pathways, including the adipocytokine signaling pathway and cAMP signaling pathway (Fig. 4H).

Correlation analysis revealed significant relationships among 23 differentially abundant bacterial species, 45 shared DEGs, and 31 shared differentially expressed metabolites (Supplementary Fig. 3A,B). These significantly correlated bacteria, metabolites, and DEGs were integrated into a joint network (Fig. 5A), identifying 35 hub genes, 30 hub metabolites, and 28 hub microorganisms. Notable relationships were identified, such as between g_Acinetobacter and Helianyl octanoate, as well as between Retinyl beta-glucuronide and Gm14002.

Fig. 5.

Multi-Omics association analysis. (A) Joint network representation incorporating differential bacteria, shared differentially expressed metabolites, and shared DEGs with significant correlations, constructed using Least Absolute Shrinkage and Selection Operator penalized Cox analysis. (B) Pearson and Mantel analyses illustrating the relationships among 35 hub genes, 30 hub metabolites, and 28 hub microorganisms. (C) Metabolite-response-enzyme-gene regulatory network centered on xanthine, which was constructed using MetScape, highlighting the interconnections among metabolites, enzymes, and genes. Green circles represent key genes, blue squares denote enzymes, gray diamonds, colored in orange and purple, distinguish different types of metabolic reaction pathways, and red hexagons indicate metabolites. XDH, xanthine dehydrogenase; NP, Purine nucleoside phosphorylase; HPRT1, hypoxanthine phosphoribosyltransferase 1; GDA, guanine deaminase; PC, Phosphatidylcholine; TG, Triglyceride; PS, Phosphatidylserine; CL, Cardiolipin.

Subsequent Pearson and Mantel analyses indicated a strongly positive correlation between the metabolite Xanthine and hub microorganisms (Mantel’s R $\ge$0.4) (Fig. 5B). This prompted the construction of a metabolite-enzyme-gene regulatory network centered around Xanthine to further elucidate its molecular regulatory mechanisms. Four key genes, including XDH, NP, HPRT1, and GDA were identified (Fig. 5C). For instance, XDH was found to interact with Xanthine, catalyzing the R02103 reaction that produces Xanthine.

Furthermore, the expression levels of hub genes, metabolites, and microorganisms across groups A, B, and C were visualized in Supplementary Fig. 4A–C. Among them, the Pomc and Rn7sk genes, the g__Bacteroides and g__Clostridia_UCG-014_unclassified microorganisms, and the LysoPC (16:0) metabolite were highly expressed in all the three groups.

As previously described, Xanthine exhibited a significantly positive correlation with hub microorganisms, and its expression levels were upregulated in group B, while they were downregulated following intervention with a NAT10 inhibitor. To investigate Xanthine’s role in TH, FBX, a xanthine oxidase (XO) inhibitor, was administered to the TH mouse model. Behavioral tests indicated that after 3 days, group F (FBX-treated) exhibited a longer latency period for cold and thermal stimuli compared with group E (Fig. 6A,B). In motor function tests, including grip strength, plane righting, and cliff avoidance, group F demonstrated intermediate performance compared with group E (Fig. 6C–E). These findings suggest that FBX provides a therapeutic benefit in alleviating clinical symptoms of TH.

Fig. 6.

FBX alleviates clinical symptoms of TH in mice. (A) Effects of heat pain stimulation. (B) Effects of cold pain stimulation. (C) Assessment of forelimb grip strength. (D) Results of the plane righting test. (E) Outcomes of the cliff avoidance test. mRNA levels of (F) GDA, (G) HPRT1, (H) NP, and (I) XDH in the mouse brain on day 21 post-intervention. (J) Analysis of ac4C expression level by dot blot in the mouse brain on day 21 post-intervention. (K) IHC (scale bars: 100 μm) and (L) WB of GDA, HPRT1, NP, and XDH in the mouse brain on day 21 post-intervention. (M) ELISA of IL-6, IL-10, and IL-1$\beta$ in the mouse brain on day 21 post-intervention. For panels (A–E) and (J), n = 6; (F–I), n = 3; (K), n = 9; (L,M), n = 3. Statistical analysis: (A,B) Repeated measures ANOVA with Bonferroni post-hoc test; (C–M) One-way ANOVA with Bonferroni post-hoc test. *p 0.05, **p 0.01, ***p 0.001, and ****p 0.0001 compared with group D; #p 0.05, ##p 0.01, ###p 0.001, and ####p 0.0001 compared with group E. In each panel, the small letter D indicates the control group (blue), E indicates the TH model group (red), and F indicates the FBX intervention group (green). FBX, febuxostat; WB, western blot; ELISA, enzyme-linked immunosorbent assay.

Gene expression analysis via RT-qPCR in groups D, E, and F revealed significant differences in the levels of key genes, including XDH, NP, HPRT1, and GDA between groups D and E. After FBX intervention, XDH and HPRT1 expression levels in group F were significantly lower than those in group E, while NP expression level was significantly higher (Fig. 6F–I). Compared to group D, group E exhibited a significant increase in ac4C-modified mRNA level, while a decrease in ac4C expression level was found in group F relative to group E (Fig. 6J). IHC further confirmed these differences. Group E exhibited significantly altered expression levels of all four key genes compared with group D. This was largely restored in group F, although only XDH exhibited a statistically significant difference (Fig. 6K). Western blot analysis showed that XDH, HPRT1, and GDA exhibited consistent expression patterns across the three experimental groups, aligning with the IHC results (Fig. 6L, Supplementary Fig. 5). ELISA of inflammatory cytokines (IL-10, IL-6, and IL-1$\beta$) in mouse brain tissue revealed elevated levels in group E compared with group D, which were significantly reduced following FBX treatment (Fig. 6M). Thus, FBX may alleviate the clinical symptoms of TH in mice by modulating the expression levels of these key genes.