Integrated Multi-Omics Analyses Identify Acupuncture-Associated Changes in Gut Microbiota, Short-Chain Fatty Acids, and TNC Gene Expression in Chronic Migraine Rats

Introduction

Migraine is one of the most common neurological disorders worldwide, affecting approximately 1 billion individuals in 2019.¹ Chronic migraine (CM), a severe subtype of migraine, features recurrent, moderate-to-severe headaches that occur more than 15 days per month for over 3 months.² The progression from episodic migraine to CM is closely associated with persistent sensitization of trigeminovascular and central pain-processing pathways, together with neuroimmune dysregulation.^3,4 Beyond these neural mechanisms, migraine is increasingly recognized as a complex neurological disorder with systemic biological involvement, in which peripheral immune, metabolic, and gastrointestinal alterations may interact with central pain-processing networks, and contribute to headache persistence.^5–8

In this context, gut microbiota and microbial metabolites have attracted growing attention in migraine research. The microbiota-gut-brain axis provides a biologically plausible framework for understanding how gut microorganisms and their metabolites may communicate with the central nervous system through immune, endocrine, metabolic, and neural pathways.⁹ Gastrointestinal symptoms and disorders are frequently reported in individuals with headache or migraine,^10–12 and observational studies have reported altered gut microbial composition in patients with migraine.¹³ Mendelian randomization analyses have further suggested that specific microbial taxa, including Bifidobacteriaceae, may be associated with migraine risk.¹⁴ Among gut-derived metabolites, short-chain fatty acids (SCFAs), including acetate, propionate, and butyrate, are of particular interest because they may influence blood-brain barrier integrity, microglial activity, inflammatory signaling, and neurotransmitter-related processes.^15–17 Nevertheless, whether changes in gut microbiota and SCFA profiles are linked to central molecular alterations in CM remains insufficiently understood.

Pharmacological therapies remain the mainstay of CM management, but medication overuse, incomplete response, and adverse effects continue to be clinical concerns.¹⁸ To investigate mechanisms relevant to CM pathophysiology, the nitroglycerin (NTG)-induced rodent model has been widely used because repeated NTG administration produces migraine-like hypersensitivity and neuroinflammatory responses, providing a suitable preclinical platform for evaluating candidate interventions.¹⁹

Acupuncture is used as a non-pharmacological intervention for migraine prevention and pain modulation. Our previous clinical and preclinical studies have shown that acupuncture reduces migraine attack frequency, alleviates migraine-like hyperalgesia, modulates brain function, suppresses central inflammatory responses, and regulates peripheral exosomal miRNAs.^20–24 Other studies have indicated that acupuncture may be associated with changes in peripheral metabolites, proteins, and gastrointestinal function.^25,26 Nevertheless, acupuncture studies remain heterogeneous in treatment protocols and sham-control designs, and the biological changes accompanying acupuncture treatment in CM, particularly those involving gut microbial composition, microbial metabolites, and central gene expression, remain incompletely characterized.

To address this gap, we combined full-length 16S rRNA sequencing, targeted SCFA profiling in fecal and plasma samples, and transcriptomic analysis of the trigeminal nucleus caudalis (TNC), a key region involved in trigeminovascular pain processing, in NTG-induced CM rats treated with true or sham acupuncture. This multi-omics approach was designed to characterize acupuncture-associated changes across gut microbial taxa, SCFA profiles, and TNC gene expression, thereby providing exploratory multi-omic evidence for potential microbiota–metabolite–central transcriptomic associations in this preclinical CM model.

Methods

Animals, Randomization, and Ethics

Male Sprague-Dawley rats, aged 6–8 weeks and weighing 200–220 g, were purchased from the Experimental Animal Center of Byrness Weil Biotech Ltd (Chongqing, China). The rats were housed under standard conditions (temperature: 22–25 °C; relative humidity: 45–65%; 12-h light/dark cycle; free access to food and water) at the Experimental Animal Center of Chengdu University of Traditional Chinese Medicine (TCM) (Chengdu, China). Twenty-four rats were randomly assigned to four groups: vehicle control group (VEH), nitroglycerin injection group (NTG), NTG plus true acupuncture group (NTG+TA), and NTG plus sham acupuncture group (NTG+SA), with six rats in each group. Randomization was performed before the experiment by an independent researcher using a computer-generated random number table in SPSS 26.0, with animals allocated at a 1:1:1:1 ratio. Allocation concealment was maintained using coded cage labels. The allocation list was kept by the independent researcher and was not disclosed to the experimenters until raw data collection and primary statistical analyses had been completed.

Blinding was implemented for outcome assessment and data analysis whenever feasible. Investigators responsible for behavioral testing, sample collection, ELISA, RT-qPCR, and data analysis were blinded to group allocation or worked with coded samples. Full-length 16S rRNA sequencing, transcriptome sequencing, and GC-MS-based SCFA quantification were performed by external service providers using coded samples only, without access to group allocation information. Because acupuncture required manual administration, the acupuncture operator could not be blinded but was not involved in outcome assessment or data analysis. Treatment identities were disclosed only after the primary analyses had been completed.

All experimental procedures were approved by the Animal Ethics Committee of Chengdu University of TCM (2023DL-047) and conducted in accordance with the National Institutes of Health Guidelines for the Care and Use of Laboratory Animals, adhering to the ARRIVE guidelines.

Sample Size and Power Analysis

The sample size of six rats per group was determined by considering the 3Rs principle, experimental feasibility, and commonly used designs in controlled rodent multi-omics studies.^27,28 Although no formal a priori sample-size calculation was performed, post-hoc power analysis using G*Power version 3.1.9.7 showed large effect sizes for representative behavioral and inflammatory outcomes, including mechanical pain threshold and TNC IL-1β. The corresponding Cohen’s d values were 7.02 and 2.25, respectively, with achieved powers >0.90 at α = 0.05. For RNA-seq and 16S rRNA sequencing, six biological replicates per group were used, which is consistent with commonly used controlled rodent omics designs^27–30 and published recommendations for RNA-seq differential expression analysis.³¹ SCFA analyses were exploratory, and plasma SCFA profiling was interpreted with caution because only four samples per group passed quality control.

NTG-Induced CM Model

The CM model was established by repeated subcutaneous injections of NTG, as described in a previously published study.³² Briefly, the diluted NTG was subcutaneously injected into rats of the NTG, NTG+SA, as well as NTG+TA groups at a dose of 10 mg/kg every other day for 9 days (ie., days 1, 3, 5, 7, and 9). Rats in the VEH group received the same volume of 0.9% saline according to the same schedule. To maintain blinding during model induction, an independent researcher prepared coded NTG or saline injection solutions according to the randomization list. The investigator administering the injections was unaware of group allocation and treatment identity.

Acupuncture Treatment

The rats were immobilized and the skin over the acupoints was disinfected with 75% alcohol. Acupuncture was administered once daily for 9 consecutive days, with each session lasting 20 min. On modeling days, acupuncture was performed 1 h before NTG or saline injection. The acupuncture stimulation was given manually using four needles, each size at 0.18 mm×13 mm (Suzhou Medical Supplies Co. Ltd, China) with an insertion depth of 5 mm. After needle insertion, the needles were twisted at a small amplitude with a frequency of 120 times/min once every 5 minutes in the TA group but not in the SA group. Rats in the NTG and VEH groups were immobilized without needle insertion.

For the TA group, bilateral GB8 (Shuaigu, located on the head, 6 mm directly above the ear apex) and GB34 (Yanglingquan, located on the lower leg, 3 mm anterior and inferior to the fibular head) were selected. The selection was based on the traditional Chinese medicine principle of combining local and distal acupoints, supported by neuroanatomical rationales. GB8 is located in the temporal region, which is primarily innervated by the auriculotemporal nerve, a branch of the mandibular division (V3) of the trigeminal nerve. Acupuncture stimulation at GB8 may activate trigeminal somatosensory afferents that project to the TNC, where they converge with nociceptive inputs from the trigeminovascular system,^33,34 thereby modulating TNC neuronal activity involved in migraine pain.³⁵ GB34 lies adjacent to the common peroneal nerve. Acupuncture at GB34 attenuates systemic neuroinflammation,³⁶ modulates autonomic function,³⁷ and activates descending pain inhibitory pathways.³⁸ Previous studies have found that this acupoint combination effectively alleviates hyperalgesia and neuroinflammation in migraine models.^35,39 For the SA group, two non-acupuncture points (one located approximately 10 mm above the iliac crest, and the other located at the posterior superior iliac spine) were applied.³⁹

Behavioral Tests

Paw withdrawal threshold (g) and tail-flick latency (s) were measured to assess the development of mechanical allodynia as well as thermal hyperalgesia 2 hours before and after NTG injection, following previous protocol.³² The paw withdrawal threshold test and tail flick test were conducted on rats on days 1, 3, 5, 7, and 9,⁴⁰ with a 30-minute interval between each test. For each test, the final value was calculated as the average of three repeated measurements, with a 5 min interval between measurements.

Paw Withdrawal Mechanical Threshold Test

Von Frey apparatus (Model 2390, IITC Life Science Inc., USA) was used to measure the mechanical threshold of the rat hind paw. Individual rat was housed in a plastic cage with a wire-mesh floor for at least 10 minutes to acclimate to the testing environment. Electronic von Frey tips were applied vertically, gradually increasing the strength, to the midplantar surface of the hind paw, while avoiding the toe pads, through the mesh floor. The withdrawal threshold was automatically recorded. A positive response was defined as either an immediate withdrawal of the hind paw or paw licking. Voluntary movements associated with locomotion were not classified as withdrawal responses.⁴¹

Tail-Flick Test

Thermal nociceptive responses in the rats were evaluated using the tail-flick test. The rat’s tail (about 20 mm from the tip) was exposed to infrared radiant heat emitted by the tail-flick instrument (SW-200, Chengdu Taimeng Software Co., Ltd., China) at a power level of 30%. The latency of withdrawal, indicated by tail flicking, was measured and recorded. A cutoff time of 10 seconds was implemented to prevent tissue damage.

Enzyme-Linked Immunosorbent Assay (ELISA)

On day 10, rats were anesthetized with 1% pentobarbital sodium (40 mg/kg) via intraperitoneal injection and then euthanized by decapitation. Samples of plasma, feces, and brain were collected for further experiments. The trigeminal nucleus caudalis (TNC) region in brain tissue was extracted, and the levels of TNF-α, IL-1β, as well as IL-6 in TNC tissues were quantified via ELISA kits (Jingmei, Jiangsu, China), following the manufacturer’s instructions.

Full-Length 16S rRNA Sequencing Analysis

Following the manufacturer’s instructions, the TGuide S96 Magnetic Soil/Stool DNA Kit (Tiangen Biotech Co., Ltd., China) was used to extract the total genomic DNA from fecal samples. Amplification of the 16S rRNA gene’s V1–V9 hypervariable regions was carried out via primers (27F:AGRGTTTGATYNTGGCTCAG;1492R:TASGGHT ACCTTGTTASGACTT). Following measurement, the PacBio Sequel II platform (Beijing Biomarker Technologies Co., Ltd., Beijing, China) was used to pool and sequence the amplicons at normalized equimolar concentrations.

After removing low-quality sequences from raw data, high-quality reads were clustered into operational taxonomic units (OTUs) based on a 97% identification using USEARCH (version 10.0). Alpha-diversity indexes, encompassing Chao1, ACE, Shannon, as well as Simpson, were employed to assess species richness as well as evenness. Beta-diversity analysis was conducted to assess the complexity of species composition among different groups. In our research, principal coordinate analysis (PCoA) and the Bray-Curtis dissimilarity metric were applied to visualize the differences between the groups. Subsequently, the composition of the gut microbiota was further analyzed. Characteristic taxa were identified using Linear discriminant analysis effect size (LEfSe) analysis⁴² based on Linear Discriminant Analysis (LDA) with a minimum score of 3.

Transcriptome Sequencing and Differential Expression Analysis

The TRIzol Reagent (Life Technologies, California, USA) was used to homogenize the TNC area in a rat brain in order to isolate total RNA. The NanoDrop 2000 device (Thermo Fisher Scientific, Wilmington, USA) was used to measure the concentration and purity of RNA. The RNA Nano 6000 Assay Kit was used on an Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA) to assess the integrity of the RNA. Following the manufacturer’s instructions, 1 μg of RNA per sample was used to create the library using the Hieff NGS Ultima Dual-mode mRNA Library Prep Kit for Illumina (Yeasen Biotechnology, Shanghai Co., Ltd., China). Each sample’s sequences were labeled using index codes. To produce 150 bp paired-end reads, the libraries were sequenced using an Illumina NovaSeq platform.

Transcriptome data were analyzed using BMKCloud (https://www.biocloud.net). First, low-quality reads and sequencing adapters were eliminated from the raw data. HISAT2 was used to align clean reads to the rat genome sequences, and StringTie was used for transcript assembly and quantification.⁴³ Fragments per kilobase of transcript per million fragments (FPKM) normalization was used to determine the levels of gene expression. With the following cut-off criteria, edgeR was utilized to find candidate differentially expressed genes (DEGs) in a group comparison: |log2FC| > 0.585, p value < 0.05. KOBAS conducted a pathway analysis using the Kyoto Encyclopedia of Genes and Genomes (KEGG).⁴⁴

RT-qPCR Analysis

Total RNA was reverse-transcribed into cDNA using MonScript™ RTIII All-in-One Mix with dsDNase (Monad, China; Cat. No. MR05101M) according to the manufacturer’s instructions. Quantitative real-time PCR was performed using SYBR^® Green Pro Taq HS Premix (Agbio, China; Cat. No. AG11701) on a Q225 real-time PCR system (Monad, China). Each 10 μL reaction contained 5 μL of 2× SYBR Green Pro Taq HS Premix, 1 μL of forward and reverse primer mix, 2 μL of cDNA template, and nuclease-free water to a final volume of 10 μL. The PCR amplification program was as follows: 95°C for 5 min, followed by 40 cycles of 95°C for 10s and 60°C for 30s. Melting curve analysis was performed from 65°C to 95°C with increments of 0.5°C every 5s. Relative gene expression levels were normalized to GAPDH and calculated using the 2^−ΔΔCt method. Each sample was analyzed in technical triplicate, with six biological samples included in each group. Primer sequences are provided in Supplementary Table 1.

Measurement of Fecal and Plasma Short-Chain Fatty Acids (SCFAs)

Blood samples were attained by cardiac puncture at the site of maximal heart pulsation. Plasma was collected from the supernatant after centrifugation. Gas chromatography-mass spectrometry (GC-MS) was leveraged to quantify SCFA levels, including isobutyric acid, propionic acid, acetic acid, butyric acid, isovaleric acid, valeric acid, and caproic acid, in fecal as well as plasma samples. The samples were prepared according to the protocol previously published.⁴⁵

Using helium as the carrier gas at a flow rate of 1 mL/min, the GC analysis was carried out using a TRACE 1310 GC (Thermo Fisher Scientific, USA) fitted with a capillary column (Agilent HP-INNOWAX: 30 m × 0.25 mm ID × 0.25 μm). The protocol served as the basis for programming the injection method, rate, and volume in addition to the column temperature.⁴⁵ Metabolites were detected by mass spectrometry using the electron impact ionization mode on an ISQ LT (Thermo Fisher Scientific, USA). With an electron energy of 70 eV, the selected ion monitoring (SIM) mode was employed.

Exploratory Multi-Omics Correlation Analysis

Potential relationships between the species-level relative abundance of gut microbiota and the abundance of fecal SCFA metabolites as well as between the species-level gut microbiota and the level of mRNA expression in the TNC brain region were investigated using Pearson correlation analysis. Correlation analyses were performed using coded datasets before group identities were disclosed. Because the multi-omics correlation analysis was performed for exploratory purposes and the sample size was limited, P values were not adjusted for multiple comparisons. Therefore, unadjusted P values are reported for the correlation analyses. Correlations with an absolute Pearson correlation coefficient |r| > 0.8 and P < 0.05 were retained for exploratory interpretation.^46,47 These criteria were used to identify potential cross-omics associations that may warrant further validation, rather than to draw confirmatory conclusions or infer causal relationships.

Statistical Analysis

Data are presented as mean ± SD unless otherwise indicated. Analyses were performed using GraphPad Prism 8.0. Normality and homogeneity of variance were assessed using the Shapiro–Wilk and Brown–Forsythe tests, respectively. For ELISA, SCFA, and RT-qPCR data, differences among groups were analyzed by one-way ANOVA with Tukey’s post hoc test or, when assumptions were not met, by the Kruskal–Wallis test with Dunn’s post hoc test. Behavioral data were analyzed using two-way repeated-measures ANOVA, with group and time as factors. Mauchly’s test was used to assess sphericity, and the Greenhouse–Geisser correction was applied when required. Post hoc comparisons were adjusted using Bonferroni correction. SCFA analyses were considered exploratory metabolic analyses, no individual SCFA was prespecified as a primary SCFA endpoint. Because only four plasma samples per group passed quality control, plasma SCFA results were interpreted with caution as supportive exploratory findings rather than confirmatory evidence. Effect estimates were reported as mean differences (MD) or mean rank differences with 95% CIs. Adjusted P < 0.05 was considered statistically significant.

Results

TA Alleviated Hyperalgesia and Neuroinflammation in CM Rats

NTG-induced migraine is a reliable animal model for investigating the pathogenesis and mechanisms of migraine.⁴⁸ We used this model to validate acupuncture efficacy on the modulation of mechanical allodynia and thermal hyperalgesia by conducting Von Frey test and tail flick test. The workflow is schematically illustrated in Figure 1A.

Figure 1 Experimental design and effects of acupuncture on hyperalgesia and neuroinflammation in NTG-induced CM rats. (A) Schematic illustration of the NTG-induced CM rat model and acupuncture treatment regimen; (B) Basal and post-treatment hind paw mechanical threshold; (C) Basal and post-treatment tail-flick latency; (D) Levels of IL-1β, IL-6, and TNF-α in TNC tissues quantified by ELISA. Data are presented as mean ± SD. In (B and C) data points represent group means at each time point, and statistical indicators denote group comparisons at each time point; in (D) each dot represents an individual rat. For (B and C) * indicates comparisons vs. VEH, # indicates comparisons vs. NTG, and $ indicates comparisons vs. NTG+TA. **adjusted P<0.01, ***adjusted P<0.001, ****adjusted P<0.0001; #adjusted P<0.05, ##adjusted P<0.01, ###adjusted P<0.001, ####adjusted P<0.0001; $adjusted P<0.05, $$adjusted P<0.01, $$$adjusted P< 0.001, $$$$adjusted P< 0.0001. For (D) brackets indicate the compared groups, and ampersands above brackets indicate adjusted P values for the bracketed comparisons. && adjusted P<0.01, &&& adjusted P<0.001, &&&& adjusted P<0.0001. (A) Created in BioRender. Liu, L. (2026) https://BioRender.com/e981rc3.

On day 1, a baseline assessment (basal response) was established prior to treatment and NTG injection. Post-treatment responses were evaluated after treatment followed by NTG injection, measuring hind paw mechanical threshold and tail-flick latency. No significant differences were seen among the four groups (adjusted P > 0.05, Supplementary Table 2). However, on day 5, a notable difference emerged between the NTG and VEH groups in both basal and post-treatment evaluations, with a time-dependent increase observed in the following days (adjusted P < 0.01, Supplementary Table 2). Compared to the NTG group, the NTG+TA group exhibited a significantly increased mechanical threshold (Figure 1B) and prolonged tail-flick latency in both basal and post-treatment evaluations from day 5 to 9 (Figure 1C) (adjusted P < 0.05, Supplementary Table 2). Conversely, these changes were not observed in the NTG+SA group, which displayed a similar trend to the NTG group.

At the Day 9 post-treatment assessment, the NTG group displayed a substantial reduction in mechanical threshold (MD = –15.39, 95% CI [–17.17, –13.61], adjusted P < 0.0001) and tail-flick latency (MD = –2.51, 95% CI [–2.85, –2.17], adjusted P < 0.0001) vs. VEH. The NTG+TA group showed significant improvements compared to the NTG group in mechanical threshold (MD = 13.78, 95% CI [9.00, 18.56], adjusted P = 0.0004) and tail-flick latency (MD = 1.74, 95% CI [0.79, 2.69], adjusted P = 0.0034). No significant differences were observed between the NTG+SA and NTG groups (adjusted P > 0.05, Supplementary Table 2).

Furthermore, we measured the TNC region’s levels of the proinflammatory cytokines TNF-α, IL-1β, and IL-6. When compared to the VEH control, NTG caused noticeably higher values of these variables (TNF-α: MD = 24.58, 95% CI [12.22, 36.93], adjusted P = 0.0001; IL-1β: MD = 3.17, 95% CI [1.09, 5.24], adjusted P = 0.002; IL-6: MD = 12.17, 95% CI [6.57, 17.77], adjusted P < 0.0001) (Figure 1D). Interestingly, TA treatment decreased their levels (TNF-α: MD = −22.39, 95% CI [−34.75, −10.04], adjusted P = 0.0003; IL-1β: MD = −2.89, 95% CI [−4.97, −0.82], adjusted P = 0.0045; IL-6: MD = −8.11, 95% CI [−13.72, −2.51], adjusted P = 0.0032), although SA treatment did not result in any discernible downregulation (adjusted P>0.05, Supplementary Table 3). This suggests that acupuncture may have an anti-neuroinflammatory impact in this CM model.

TA Treatment Was Associated with Changes in Gut Microbiota Composition in CM Rats

Full-length 16S rRNA gene sequencing was used to characterize gut microbiota profiles associated with TA treatment in CM rats. Fecal samples from 24 rats (n=6 per group) were analyzed, and the rarefaction curve as well as Shannon Index curve remained balanced as the number of reads increased (Supplementary Figure 1A), indicating sufficient sequencing data and reasonable sequencing depth.

Alpha-diversity indexes, encompassing Chao1 and ACE, showed no remarkable difference (P > 0.05) among the four groups (Supplementary Figure 1B), suggesting that TA treatment was not associated with obvious changes in the richness of the bacterial community. Only the Simpson index showed a significant decrease in the NTG+TA group compared to the NTG group (P < 0.05), implying that TA treatment was associated with altered community evenness. Beta-diversity analysis revealed a clear separation of the samples (scattered points, Figure 2A) between the VEH and NTG groups, as illustrated in the PCoA plot based on Bray-Curtis distance and PERMANOVA test (P=0.001, R²=0.2925). The samples in the NTG+TA group were clustered closer to those in the VEH group, while those in the NTG+SA group were closer to the NTG group, suggesting that TA treatment was associated with a partial shift of the gut microbiota profile toward that of VEH rats.

Figure 2 TA-associated changes in gut microbiota composition in NTG-induced CM rats. (A) PCoA plot displaying the grouped discrimination among the groups based on the relative abundance of OTUs, with individual samples color-coded for VEH group (grey square), NTG group (Orange diamond), NTG+TA (turquoise triangle), and NTG+SA (dark red, dot); (B) Component proportion of bacterial phylum and species in each group; (C) Characteristic taxa identified in the NTG+TA vs. NTG group by LEfSe analysis; (D) Relative abundance of selected characteristic taxa including Lachnoclostridium, Ruminococcus bromii, Oscillospiraceae, Lactobacillus murinus, Alistipes massiliensis, Parabacteroides goldsteinii, and Ligilactobacillus murinus depicting differences among groups. Bars represent mean ± SD, and each dot represents an individual fecal sample from one rat.

Next, we analyzed proportional changes in gut microbiota composition at different taxonomic levels between groups, focusing specifically on the phylum and species levels (Figure 2B). Additional taxa are depicted in Supplementary Figure 1C. At the phylum level, Bacteroidota and Firmicutes remained dominant across all groups, with a combined relative abundance exceeding 60%. In the NTG+TA group, Verrucomicrobiota showed a notable decrease. At the species level, the NTG group showed decreased relative abundance of Lactobacillus murinus and increased relative abundances of unclassified Muribaculaceae, Lactobacillus acidophilus, Lactobacillus johnsonii, and Akkermansia muciniphila compared with the VEH group. Although Lactobacillus acidophilus and Lactobacillus johnsonii have been reported as probiotic-associated taxa in some contexts, their enrichment in the NTG group occurred together with hyperalgesia and increased inflammatory responses, suggesting that these changes may reflect context-dependent, CM-associated microbial remodeling rather than a protective effect. TA treatment was associated with a microbial profile closer to that of the VEH group and a higher relative abundance of Lactobacillus murinus.

LEfSe analysis identified bacterial taxa that significantly differentiated the NTG+TA group from the NTG group, as illustrated in the histogram (Figure 2C). Five OTUs in the NTG+TA group and 17 OTUs in the NTG group met the LDA significance threshold (LDA>3). Taxa enriched in the NTG group included Lachnoclostridium, Ruminococcus, and Oscillospiraceae, whereas taxa enriched in the NTG+TA group included Lactobacillus murinus, Alistipes massiliensis, Parabacteroides goldsteinii, and Ligilactobacillus murinus. The relative abundances of these taxa across groups are displayed in Figure 2D. The distribution of other taxa among the additional groups is shown in Supplementary Figure 1D.

Exploratory Transcriptomic Profiling of the TNC Following TA Treatment in CM Rats

Gene expression profiles of migraine patients exhibit distinctions from those of the general population. To investigate whether TA treatment was associated with changes in gene expression, we conducted transcriptome profiling in the TNC regions of CM rats. Each sample yielded approximately 40 million clean reads, with about 96.65% of reads mapped to the rat genome following quality control analysis. A Venn diagram visually represented the overlapped and distinct genes among the groups (Figure 3A). Applying a discovery threshold (|fold change| > 1.5 and unadjusted p value < 0.05), a total of 1482 candidate DEGs were identified. Compared to the VEH group, the NTG group showed 187 downregulated and 72 upregulated candidate DEGs (Figure 3B). TA treatment was associated with 642 upregulated and 501 downregulated candidate DEGs compared to the NTG group (Figures 3B), including S100a4 and Vamp8, which were both upregulated (S100a4: log2FC = 0.72, P = 0.0039; Vamp8: log2FC = 0.64, P = 0.0063). We further validated the expression of S100a4 and Vamp8 by RT-qPCR. Consistent with the transcriptomic findings, RT-qPCR analysis showed that, compared with the VEH group, the NTG group exhibited significantly decreased mRNA levels of both S100a4 and Vamp8. Conversely, compared with the NTG group, the NTG+TA group showed significantly increased expression of both genes (all adjusted P < 0.0001, Figure 3D, Supplementary Table 4). Additionally, when compared to the SA group, the NTG+TA group showed 344 upregulated and 327 downregulated candidate DEGs (Figures 3B). SA treatment, on the other hand, was associated with 55 upregulated and 52 downregulated candidate DEGs compared to the NTG group (Figure 3B). The volcano plot displayed upregulated genes in red and downregulated genes in blue.

Figure 3 TNC transcriptomic profiles in NTG-induced CM rats following TA treatment. (A) Venn diagram showing overlapping and unique candidate DEGs among pairwise comparisons; (B) Volcano plots displaying exploratory candidate DEGs identified by pairwise group comparisons. Upregulated (red), downregulated (blue); (C) Selected KEGG enriched pathway presented for comparisons between the VEH and NTG groups, and the NTG and NTG+TA groups; (D) RT-qPCR validation of the expression levels of selected candidate DEGs, including S100a4 and Vamp8, in TNC tissues; Each dot represents an individual rat. Brackets indicate the compared groups. ***adjusted P < 0.001, ****adjusted P < 0.0001.

KEGG pathway analysis further explored the potential biological functions associated with these candidate DEGs. The candidate DEGs between the VEH and NTG groups revealed enrichment in pathways such as neuroactive ligand-receptor interaction, long-term depression, cGMP-PKG signaling pathway, and platelet activation (Figure 3C, Supplementary Table 5). In contrast, the candidate DEGs between the NTG and NTG+TA groups were mainly enriched in pathways, including serotonergic synapse, complement and coagulation cascades, chemokine signaling pathway, cAMP signaling pathway, TGF-beta signaling pathway, NOD-like receptor signaling pathway, among others (Figure 3C, Supplementary Table 5).

In addition, 28 immune-related candidate DEGs identified in the comparison of NTG+TA and NTG groups (Table 1) were used for subsequent integrative analysis to explore potential associations between gene expression and gut microbiota.

Table 1 28 Immune-Related Candidate DEGs Identified in the NTG+TA Group Compared to the NTG Group

Exploratory Analysis of Fecal and Plasma SCFA Profiles Following TA Treatment in CM Rats

The fecal SCFA levels, including isobutyric acid, propionic acid, acetic acid, butyric acid, isovaleric acid, valeric acid, as well as caproic acid, are shown in Figure 4A and Supplementary Table 6. Our findings revealed higher levels of acetic acid in the NTG vs. VEH group, with a subsequent decrease observed only after TA treatment (MD = −1028, 95% CI [−1815, −241.4], adjusted P = 0.0078). Similarly, propionic acid levels were elevated in the NTG group compared to the VEH group but decreased after acupuncture treatment (adjusted P>0.05, Supplementary Table 6). Due to quality control issues, only four plasma samples per group qualified for the final analysis. Nevertheless, the distribution of SCFAs in plasma closely aligned with their level in fecal samples across all groups (Figure 4B, Supplementary Table 6).

Figure 4 Fecal and plasma SCFA profiles in NTG-induced CM rats following acupuncture treatment. (A) Fecal concentrations of isobutyric acid, propionic acid, acetic acid, butyric acid, isovaleric acid, valeric acid, and caproic acid in the VEH, NTG, NTG+TA, and NTG+SA groups; (B) Plasma concentrations of isobutyric acid, propionic acid, acetic acid, butyric acid, isovaleric acid, valeric acid, and caproic acid in the VEH, NTG, NTG+TA, and NTG+SA groups. n=6 per group for fecal SCFAs and n=4 per group for plasma SCFAs. Data are presented as mean ± SD, and each dot represents an individual fecal or plasma sample from one rat. *adjusted P < 0.05, **adjusted P < 0.01, ***adjusted P < 0.001. Only significant bracketed comparisons are annotated in the figure; full adjusted P values and effect estimates are provided in Supplementary Table 6.

Exploratory Correlation Analysis Among Gut Microbiota, Fecal SCFAs, and Immune-Related Genes

Given the well-established reciprocal influence between gut microbiota and metabolic profiles,^49,50 we conducted a correlation analysis between the gut microbiota and fecal SCFAs. The heatmap (Figure 5A, Supplementary Table 7) shows that the identified bacterial taxa at the species level were associated with SCFAs to varying degrees. Notably, a remarkable negative correlation was seen between acetic acid and Parabacteroides goldsteinii (P=0.0029). Propionic acid exhibited a negative correlation with Bacteroides uniformis and Ruminococcus sp WAL 17306 (P<0.05).

Figure 5 Exploratory correlation analysis of gut microbiota with fecal SCFAs and TNC candidate genes within the NTG+TA group. Correlation heatmaps showing nominal exploratory Pearson correlations between gut microbial taxa and (A) fecal SCFAs and (B) TNC candidate genes. *unadjusted P < 0.05, **unadjusted P < 0.01. P values in the correlation analyses were not adjusted for multiple comparisons. (C) Representative linear fitting plots showing associations among gut microbial taxa, fecal SCFAs, and TNC candidate genes.

Similarly, the exploratory pairwise correlations between bacterial taxa and the 28 immune-related candidate genes listed in Table 1 are shown in Figure 5B and Supplementary Table 8. Parabacteroides goldsteinii showed a positive correlation with S100a4, Vamp8, and Ccl19, while displaying a negative correlation with Foxo3, Prkacb, and Map2k4 (Figure 5C, P<0.05). Lactobacillus murinus showed a negative correlation with Adcy1 (Figure 5C, P<0.05).

Collectively, these findings suggest potential exploratory associations among specific bacterial taxa, fecal SCFAs, and TNC transcriptomic alterations, thereby providing hypotheses for future mechanistic studies on acupuncture intervention in migraine.

Discussion

Acupuncture has been reported to modulate pain-related, inflammatory, and neuroimmune processes through multi-level and multi-target biological responses. However, the biological basis underlying these effects remains incompletely understood, particularly with respect to how peripheral microbial and metabolic alterations may be associated with central neuroinflammatory responses.^20–24 In the present study, we used an integrative multi-omics approach to characterize acupuncture-associated changes across gut microbiota, SCFA profiles, and TNC transcriptomic responses in an NTG-induced CM rat model, providing exploratory evidence for the possible involvement of the microbiota-gut-brain axis (Figure 6). We observed that TA, but not SA, attenuated neuroinflammatory responses and hyperalgesia in NTG-induced CM rats. Notably, TA and SA were associated with distinct SCFA levels, gut microbiota compositions, and gene expression patterns in the TNC of CM rats. Our findings suggest that the alleviation of migraine-like hyperalgesia by TA in this preclinical model was associated with reduced inflammatory markers, as evidenced by decreased levels of pro-inflammatory cytokines and alterations in neuroinflammation-related signaling pathways in brain tissue. In addition, TA was associated with changes in gut microbiota composition and SCFA profiles. Moreover, alterations in gut microbiota were correlated with changes in the brain transcriptome, particularly the expression of the genes S100a4 and Vamp8, supporting a possible association between TA-associated microbiota remodeling and TNC transcriptomic changes along the microbiota-gut-brain axis.

Figure 6 Integrated multi-omics overview of TA-associated alterations in the microbiota-gut-brain axis in NTG-induced CM rats. The schematic summarizes TA-associated alterations in TNC neuroinflammation, TNC transcriptomic signatures, fecal gut microbiota and SCFA profiles, and exploratory plasma SCFA profiles. Upward and downward arrows indicate increases and decreases, respectively, in the NTG+TA group compared with the NTG group. Solid arrows indicate statistically significant differences, whereas dashed arrows indicate non-significant trends. Dashed curved lines indicate exploratory associations between omics layers. Integrated correlation analysis showed that Parabacteroides goldsteinii was positively correlated with S100a4 and Vamp8 and negatively correlated with fecal acetic acid. *unadjusted P < 0.05, **unadjusted P < 0.01. P values in the correlation analyses were not adjusted for multiple comparisons. Created in BioRender. Liu, L. (2026) https://BioRender.com/rfm9dwp.

Repeated NTG administration was used to induce CM-like pain hypersensitivity and neuroinflammatory responses in rats. Although this model is useful for investigating migraine-related mechanisms and candidate interventions, it does not fully recapitulate the clinical heterogeneity of human CM; therefore, the findings require further validation before clinical extrapolation. NTG has been shown to decrease mechanical thresholds in mouse hind paws⁵¹ and increase the levels of IL-1β, IL-6, and TNF-α.^52–54 Dysregulated inflammation can contribute to enhanced pain sensitivity and hyperalgesia. Acupuncture has been reported to reduce neuroinflammation by diminishing the release of trigeminal-activated neuropeptides, inhibiting dural immune cell activation, and downregulating inflammatory mediators.⁵⁵ Consistent with these findings, our results showed that TA significantly improved mechanical and thermal hyperalgesia, and reduced IL-1β, IL-6, and TNF-α levels in CM rats.

Exploratory transcriptome profiling of the TNC region suggested that TA-associated candidate genes were linked to pathways involved in neurotransmission and neuroimmune regulation, including serotonergic synapse, cAMP signaling, TGF-β signaling, chemokine signaling, and NOD-like receptor signaling.^56–58 Serotonergic synapse and cAMP signaling may be involved in pain-related neurotransmission and neuronal excitability,^59,60 whereas TGF-β, chemokine, and NOD-like receptor signaling are closely associated with immune and inflammatory regulation.^61–63 However, these pathway results should be interpreted as exploratory functional annotations of TA-associated transcriptomic changes, and further validation is needed to clarify their biological relevance. Although substance P signaling, mainly encoded by Tac1 and mediated through neurokinin receptors such as Tacr1, has recently been reconsidered in migraine research,⁶⁴ Tac1/Tacr1-related signals were not identified as prominent transcriptomic features in the present dataset. Therefore, substance P-related pathways were not interpreted as part of the current acupuncture-associated multi-omic findings and should be examined directly in future studies.

Increasing evidence has linked gut microbial dysbiosis to migraine pathophysiology. In our study, Chao1 and ACE indices did not differ significantly among groups, whereas the Simpson index differed between the NTG and NTG+TA groups, suggesting that TA-associated microbiota changes may involve diversity/evenness-related community structure rather than overall microbial richness. Firmicutes and Bacteroidota are predominant bacterial phyla in the gut and have been widely studied in relation to host metabolism, immune regulation, and inflammatory conditions. Alterations in the Firmicutes/Bacteroidetes ratio have been reported in several inflammatory and metabolic disorders, although its biological interpretation remains context-dependent.⁶⁵ In the present study, TA treatment increased the relative abundance of Firmicutes and decreased that of Bacteroidota, shifting the overall bacterial composition toward that of the VEH group. However, because 16S rRNA sequencing provides relative rather than absolute abundance information, these compositional changes should be interpreted cautiously. LEfSe analysis further identified several taxa with differential relative abundances between the NTG and NTG+TA groups. The NTG group had higher levels of Lachnoclostridium, Ruminococcus, and Oscillospiraceae. Lachnoclostridium, a genus within the Lachnospiraceae family, has been reported as the most abundant group in migraine-afflicted mice.⁶⁶ Also, a Mendelian randomization study has linked Lachnoclostridium to neurodegenerative disease.⁶⁷ In contrast, the NTG+TA group was enriched in Parabacteroides goldsteinii, Ligilactobacillus murinus, Lactobacillus murinus, and Alistipes massiliensis. Among these taxa, Lactobacillus murinus and Alistipes massiliensis have been associated with anti-inflammatory and immunomodulatory effects in various experimental settings.^68–70 In our study, TA treatment increased the relative abundance of Lactobacillus murinus, which was decreased in the NTG group. These findings suggest a possible association between TA-induced microbiota remodeling and reduced inflammatory responses, although causality remains to be determined. Interestingly, Lactobacillus acidophilus and Lactobacillus johnsonii, often regarded as probiotic-associated bacteria,⁷¹ were enriched in the NTG group. This apparent paradox may reflect the context- and strain-dependent effects of Lactobacillus species,^72,73 a compensatory response to NTG-induced inflammatory or metabolic disturbance, or proportional changes inherent to relative-abundance-based 16S rRNA sequencing.⁷⁴ Thus, the enrichment of these taxa should be interpreted as part of the dysbiotic microbial pattern in CM rats. Instead, TA increased Lactobacillus murinus and shifted the overall microbial profile toward the VEH group, suggesting a possible microbiota-remodeling association rather than a general increase in all probiotic-associated bacteria.

SCFAs, as microbiota-derived metabolites, can influence gut-brain communication through immune, endocrine, neural, and metabolic pathways, including modulation of serotonin-related signaling.⁷⁵ Elevated serotonin levels in the peripheral and CNS are implicated in the pathogenesis of migraine.⁵⁵ In the present study, exploratory KEGG analysis based on nominal-threshold candidate DEGs suggested involvement of the serotonergic synapse pathway, and TA was accompanied by altered SCFA profiles, particularly lower fecal acetic acid. Acetic acid has previously been reported to be increased in CM rats subjected to repeated inflammatory stimulation.⁷⁶ In our correlation analysis, fecal acetic acid was negatively correlated with Parabacteroides goldsteinii, a taxon enriched in the NTG+TA group. Conversely, taxa enriched in the NTG group, including unclassified Muribaculaceae and Lactobacillus acidophilus, have been associated with SCFA production in previous studies.^77,78 These findings suggest a possible association among TA treatment, gut microbial composition, SCFA alterations, and serotonergic pathway annotations.

We further explored how acupuncture’s modulation of the gut microbiota might influence neuroinflammatory pathways in the brain, which are increasingly recognized as contributing factors to neurological disorders through transcriptomic and epigenetic regulation.^79–83 Transcriptomic data from the TNC region identified a set of candidate DEGs, including S100a4 and Vamp8, between NTG and NTG+TA groups. RT-qPCR validation further supported the TA-associated changes in S100a4 and Vamp8 expression, strengthening the reliability of the transcriptomic findings. Notably, these genes exhibited significant correlations with specific gut microbial taxa, with Parabacteroides goldsteinii showing a positive correlation with both S100a4 and Vamp8. Nevertheless, the observed associations among the gut microbiota, SCFAs, and gene expression remain correlative and warrant further mechanistic investigation. VAMP8 is involved in synaptic vesicle fusion with the presynaptic membrane and has been implicated in neuroinflammation relief.⁸⁴ It also mediates the exocytosis of MUC2 mucin from colonic goblet cells, maintaining intestinal homeostasis,⁸⁵ with its loss exacerbating colitis in mice.⁸⁶ S100A4, a Ca²⁺-binding protein, has a crucial role in the function of different neural cells, encompassing astrocytes, microglia, infiltrating cells, as well as neurons. It contributes to neuronal plasticity, survival, as well as the CNS development.⁸⁷ S100A4 has been implicated in inflammatory signaling, including TGF-β/NF-κB-related processes, in a context-dependent manner.^88,89 Interestingly, Parabacteroides goldsteinii has been shown to enhance gut recovery from injury by reducing inflammatory cytokine release via inhibition of the MyD88/NF-κB signaling pathway.⁹⁰ In our data, Parabacteroides goldsteinii was enriched in the NTG+TA group and was positively correlated with S100a4 and Vamp8 expression. Together, these findings suggest a possible association between TA-associated microbiota remodeling and TA-associated immune‑related transcriptomic changes in the TNC.

Although this study characterized acupuncture-associated changes in gut microbiota, SCFA profiles, inflammatory responses, and TNC transcriptomics, several limitations should be noted. First, 16S rRNA sequencing provides compositional information but does not directly assess microbial function, strain-level variation, or absolute bacterial abundance. Second, the multi-omics integration was exploratory and correlation-based. Because of the modest sample size and the use of unadjusted P values in the cross-omics correlation analyses, these findings should be interpreted as hypothesis-generating associations rather than causal evidence. Behavioral–molecular correlation analyses were also not performed because behavioral outcomes were assessed repeatedly over time, whereas omics endpoints were measured at euthanasia. Third, no formal a priori sample-size calculation was performed for all omics endpoints, and plasma SCFA profiling was limited by the reduced number of samples passing quality control; therefore, plasma SCFA results were considered supportive exploratory data. Fourth, potential cage effects, environmental influences, and inter-individual microbiota variability could not be fully excluded. Finally, only male rats were included, and the NTG-induced model recapitulates selected migraine-like features but does not fully reproduce the clinical heterogeneity of human chronic migraine. Future studies should include larger and sex-balanced cohorts, prospective sample-size calculations, metagenomic sequencing, absolute microbial quantification, and mechanistic validation such as germ-free animal models or fecal microbiota transplantation.

Conclusions

Our findings suggest that in the NTG-induced CM rat model, TA reduces migraine-like hyperalgesia and neuroinflammatory responses, accompanied by changes in gut microbiota composition, SCFA profiles, and TNC gene expression. These findings provide preliminary associative evidence linking acupuncture-associated microbial, metabolic, and transcriptomic changes in this preclinical model, and warrant further mechanistic and clinical validation.