Eight rats were randomly selected from each group. After anesthetizing the rats, the abdominal cavity was surgically opened, and blood was collected from the abdominal aorta. A volume of 100 μL of the collected blood was transferred into a 1.5 mL centrifuge tube, followed by the addition of 400 μL of an extraction solution composed of acetonitrile and methanol in a 1:1 ratio. This solution contained four internal standards, including L-2-chlorophenylalanine and its three isotopes at a concentration of 0.02 mg/mL. The mixture was subjected to vortex mixing for 30 s and subsequently underwent low-temperature ultrasonic extraction at 5°C and 40 KHz for 30 min. The samples were then allowed to stand at −20°C for 30 min. Centrifugation was conducted at 4°C for 15 min at 13,000 g. The supernatant was carefully removed and evaporated to dryness under a stream of nitrogen gas. The residue was reconstituted in 100 µL of a re-solution consisting of acetonitrile and water in a 1:1 ratio. This reconstituted solution was subjected to low-temperature ultrasonic extraction for 5 min at 5°C and 40 kHz, followed by centrifugation at 13,000 g for 10 min at 4°C. The resulting supernatant was transferred to a vial equipped with an internal cannula for subsequent analysis.

The mid-rectal feces were collected from rats after blood sampling, and 100 mg of solid sample was taken in a 2 mL centrifuge tube, to which a 6 mm diameter grinding bead was added. Metabolite extraction was conducted using 800 μL of an extraction solution composed of methanol and water in a 4:1 volume ratio, which included four internal standards such as L-2-chlorophenylalanine and its three isotopes at a concentration of 0.02 mg/mL. The sample solutions were subjected to grinding using a frozen tissue grinder for 6 min at −10°C and 50 Hz, followed by low-temperature ultrasonic extraction for 30 min at 5°C and 40 kHz. Subsequently, the samples were allowed to stand at −20°C for 30 min, centrifuged at 4°C and 13,000 g for 15 min, and the resulting supernatant was transferred into a vial equipped with an internal cannula for further analysis.

A 20ul equivalent volume of each sample metabolite was combined to create a quality control (QC) sample. During instrumental analysis, a QC sample was interspersed every 5 to 10 samples to assess the reproducibility of the entire analytical procedure.

The LC-MS analysis was conducted utilizing an ultra-high performance liquid chromatography tandem Fourier transform mass spectrometry system, specifically the UHPLC-Q Exactive HF-X model. A 3 μL aliquot of the sample was separated on an HSS T3 column (dimensions: 100 mm × 2.1 mm i.d., 1.8 µm particle size) prior to mass spectrometric detection. The mobile phase A comprised a solution of 0.1% formic acid in a water/acetonitrile mixture (95/5, v/v), while mobile phase B consisted of an acetonitrile/isopropanol/water mixture (47.5/47.5/5, v/v/v) with 0.1% formic acid. The chromatographic separation was carried out at a flow rate of 0.40 mL/min and a column temperature of 40°C. Mass spectrometric data acquisition was executed in both positive and negative ion scanning modes, covering a mass-to-charge ratio (m/z) range of 70–1050. The ion spray voltages were set at 3500 V for positive ions and −3000 V for negative ions. Additional parameters included a sheath gas flow rate of 50 arbitrary units, an auxiliary heated gas flow rate of 13 arbitrary units, an ion source heating temperature of 450°C, and cyclic collision energies of 20, 40, and 60 V.

Following on-boarding, LC-MS raw data were processed using Progenesis QI (Waters Corporation, Milford, United States) for baseline filtering, peak identification, integration, retention time correction, and alignment, resulting in a data matrix comprising retention time, mass-to-charge ratio, and peak intensities. Concurrently, MS and MSMS data were matched with public databases HMDB (http://www.hmdb.ca/), Metlin (https://metlin.scripps.edu/), and the Majorbio Database (MJDB) from Majorbio Biotechnology Co., Ltd. (Shanghai, China).

The data matrix, following the library search, was uploaded to the Majorbio cloud platform (https://cloud.majorbio.com) for further analysis. Initial pre-processing of the data matrix involved applying the 80% rule to address missing values; specifically, variables with more than 80% non-zero values in at least one set of samples were retained. Subsequently, missing values were imputed using the smallest values present in the original matrix. To mitigate errors arising from sample preparation and instrument instability, the response intensities of the sample mass spectrometry peaks were normalized using the sum-normalization method, resulting in normalized data matrices. Concurrently, variables exhibiting a relative standard deviation (RSD) greater than 30% in quality control (QC) samples were excluded. The remaining data were then subjected to log10 transformation to produce the final data matrix for subsequent analysis.

Genomic DNA was extracted from microbial communities utilizing the E.Z.N.A.® Soil DNA Kit (Omega Bio-tek, Norcross, GA, United States) following the manufacturer’s protocol. The integrity of the extracted DNA was assessed via agarose gel electrophoresis using a 1% agarose gel. DNA concentration and purity were quantified with a NanoDrop 2000 spectrophotometer (Thermo Scientific, United States).

The extracted DNA was employed as a template for the polymerase chain reaction (PCR) amplification of the V3-V4 variable region of the 16S rRNA gene, utilizing the primer pairs 338F (5′-ACTCCTACGGGGAGGCAGCAG-3′) and 806R (5′-GGACTACHVGGGTWTCTAAT-3′) (Liu et al., 2016). These primers incorporate a Barcode sequence. The PCR reaction mixture consisted of 4 μL of 5× TransStart FastPfu buffer, 2 μL of 2.5 mM deoxynucleotide triphosphates (dNTPs), 0.8 μL of each primer (5 μM), 0.4 μL of FastPfu polymerase, 10 ng of template DNA, and double-distilled water (ddH2O) to achieve a final volume of 20 μL. The amplification protocol included an initial denaturation at 95°C for 3 min, followed by 27 cycles of denaturation at 95°C for 30 s, annealing at 55°C for 30 s, and extension at 72°C for 30 s, with a final extension at 72°C for 10 min. The PCR products were subsequently stored at 4°C. The amplification was conducted using an ABI GeneAmp® 9,700 thermal cycler. The PCR products were resolved on a 2% agarose gel, purified using a DNA gel recovery and purification kit (PCR Clean-Up Kit, China), and quantified using a Qubit 4.0 fluorometer (Thermo Fisher Scientific, United States).

The purified PCR products were utilized for library construction employing the NEXTFLEX Rapid DNA-Seq Kit, following a series of steps: (1) junction linkage, (2) elimination of self-ligated junction fragments through magnetic bead screening, (3) enrichment of the library template via PCR amplification, and (4) recovery of PCR products using magnetic beads to obtain the final library. Sequencing was conducted on the Illumina Nextseq2000 platform in accordance with the standard protocols provided by Majorbio Bio-Pharm Technology Co. Ltd. (Shanghai, China).

The paired-end raw sequencing reads were subjected to quality control using the fastp software (https://github.com/OpenGene/fastp, version 0.19.6) and FLASH software (http://www.cbcb.umd.edu/software/flash, version 1.2.11) for read merging. The quality control process involved the following steps: (1) Bases with a quality score of 20 or lower at the ends of the reads were filtered out. A sliding window of 10 bp was applied, and if the average quality score within this window was below 20, the bases from the window onwards were truncated. Reads with a quality score of 50 bp or lower after quality control were discarded, as were reads containing ambiguous ‘N' bases. (2) Based on the overlap between paired-end reads, pairs of reads were merged into a single sequence, requiring a minimum overlap length of 10 bp. (3) The maximum allowable mismatch ratio in the overlap region of a merged sequence was set at 0.2, and sequences not meeting this criterion were excluded. (4) Samples were differentiated by barcodes and primers located at the beginning and end of the sequences, with sequence orientation adjusted accordingly. The barcode allowed for zero mismatches, while the maximum number of allowable primer mismatches was set at two. Chimeric sequences were eliminated, and Operational Taxonomic Units (OTUs) were clustered at a 97% similarity threshold utilizing UPARSE (http://drive5.com/uparse/, version 11.0.667). Sequences identified as having chloroplast or mitochondrial origins were excluded, as recommended in cases of contamination. To mitigate the effects of sequencing depth on alpha and beta diversity analyses, all samples underwent rarefaction to 20,000 sequences per sample, a recommended practice, resulting in an average Good’s coverage of 99.09% following normalization. OTU taxonomic annotation was conducted using the Silva 16S rRNA database (version 138.2) with the RDP classifier (http://rdp.cme.msu.edu/, version 2.11) at a 70% confidence threshold. Subsequent community composition analysis was performed across multiple taxonomic levels. Functional prediction of the 16S data was carried out using PICRUSt2 (version 2.2.0).

Following the acquisition of the OTUs abundance array, the gut microbiota of all rats was subjected to analysis for both Alpha and Beta diversity. The richness and diversity of the rat gut microbiota composition were evaluated using three principal ecological indices: ACE, Chao, and Shannon. Notably, the ACE and Chao indices exhibited a positive correlation with microbial species richness, while the Shannon index was positively correlated with microbial species diversity. Our findings revealed that the ACE, Chao, and Shannon indices were significantly elevated in the M group compared to the WKY group (p < 0.001), indicating a greater diversity and abundance of gut microbiota in the M group. Conversely, the diversity and abundance of the gut microbiota were markedly diminished in the EA and MA groups relative to the M group (p < 0.001), as illustrated in Figure 3.

In this study, principal coordinate analysis (PCoA) and non-metric multidimensional scaling (NMDS) were conducted utilizing the Bray-Curtis distance metric. These analytical methods effectively assessed alterations in gut microbiota structure across samples, with inter-sample distances elucidating the similarities and differences in community composition. The findings reveal that the majority of samples from the WKY and M groups were distinctly separated, indicating substantial differences in gut microbiota composition between these two groups. Conversely, the relatively small distances among samples from the WKY, MA, and EA groups suggest a high degree of similarity in microbial composition within these groups (ANOSIM, R = 0.55, p = 0.001). These results indicate that the gut microbiota in the EA and MA groups underwent modifications post-treatment, aligning more closely with the microbiota composition observed in the WKY group, as illustrated in Figure 4.

As illustrated in Figure 5, the bar charts depict the average abundance of highly prevalent abundance at both the phylum and genus levels for samples from the WKY, M, MA, and EA groups. The findings indicate that the predominant phyla in rat fecal samples from the WKY, MA, and EA groups were Firmicutes and Actinobacteriota, whereas the M group was primarily dominated by Firmicutes and Patescibacteria. At the genus level, the samples across all groups were chiefly dominated by Blautia, Lactobacillus, Allobaculum, UCG-005, and Romboutsia.

LEfSe analyses were conducted to identify key species that exhibited significant differences between groups, and a Linear Discriminant Analysis (LDA) score was calculated to distinguish signature microorganisms within each group that significantly influenced variations in colony structure. The findings indicated that microorganisms such as Collinsella, Anaerostipes, and Fusicatenibacter exerted a substantial influence on the structural differences of the bacterial flora in the normotensive control group. In group M, UCG-005 and Roseburia were identified as the signature microorganisms affecting the flora structure in this cohort of rats. In the MA group, the dominant bacteria that significantly impacted the differences in colony structure included Blautia, Allobaculum, Bifidobacterium, and Ruminococcus, among others. Similarly, in the EA group, Blautia, Romboutsia, and Bifidobacterium were the dominant bacteria influencing the structural differences in the flora, as illustrated in Figure 6.

In comparison to the WKY group, the M group exhibited a significantly reduced abundance of Blautia and Lactobacillus, alongside a significantly increased abundance of UCG-005, Bacteroides, and Helicobacter. Conversely, the EA and MA groups demonstrated a significantly higher abundance of Blautia, Bifidobacterium, and Romboutsia, as well as a significantly lower abundance of UCG-005, Bacteroides, and Helicobacter, relative to the M group. These findings suggest that both EA and MA have the potential to ameliorate gut flora disturbances in SHRs and restore balance to the gut microbiota. For detailed data, refer to Table 1 and Supplementary Tables S1-S3.

A metabolomic analysis of serum and fecal samples from four groups of rats was conducted utilizing two modes of mass spectrometry, anionic and cationic, to identify various metabolites at distinct ionic termini. Principal Component Analysis (PCA) was employed to elucidate the most accurate differences in metabolism among the groups. Figures 7A,B present the PCA score plots for the serum samples, where the motifs associated with the WKY group are distinctly separated from those of the other groups, indicating a greater metabolic divergence. In contrast, the closer proximity of the MA and EA groups suggests that the impacts of MA and EA on the serum metabolic profiles of rats are more similar. Correspondingly, the metabolic outcomes for the fecal samples, depicted in Figures 7C,D, mirror those observed in the serum samples. The model developed for the analysis of serum and fecal metabolism was validated, with an R²X (cum) value exceeding 0.5, signifying a well-fitted model, as detailed in Table 2.

The analysis of serum and fecal differential metabolites between groups was conducted using OPLS-DA VIP >1 and p < 0.05. In the M group, compared to the WKY group, 221 serum differential metabolites were identified at the cationic end and 225 at the anionic end. The top 30 metabolites, based on VIP values, were selected, including Geranic acid, PGB2, and 5,6-Dihydroxyprostaglandin F1a (Figures 8A, 9A). Similarly, fecal differential metabolites were also identified and ranked by the top 30 VIP values (Figures 8B, 9B). These findings suggest that SHRs exhibit more pronounced in vivo metabolite differences compared to normal rats. In the EA group, relative to the M group, 189 serum differential metabolites were identified at both the cationic and anionic termini. The top 30 metabolites, based on VIP values, were selected, including Taurodeoxycholic acid, Taurocholic Acid, and Deoxycholic Acid (Figures 8C, 9C). Similarly, fecal differential metabolites were identified and ranked by the top 30 VIP values (Figures 8D, 9D). In the MA group, compared to the M group, 226 serum differential metabolites were identified at the cationic end and 190 at the anionic end. The metabolites with the highest 30 VIP values, including Corchoroside B, Bz-Arg-OEt, and Mibefradil. (Figures 8E, 9E). Similarly, fecal differential metabolites were identified and screened based on the top 30 VIP values (Figures 8F, 9F). These findings suggest that both EA and MA influence the metabolite profiles in the serum and feces of SHRs.

By examining the intersection of serum metabolites, we identified those that were decreased in group M compared to WKY and increased in group EA compared to group M, as well as those that were increased in group M compared to WKY and decreased in group EA compared to group M. This analysis yielded a total of 182 callback serum metabolites from the EA group (Figures 10A,B; Supplementary Table S4). Similarly, 198 callback serum metabolites were identified for the MA group (Figures 10C,D; Supplementary Table S5), of which 136 were common to both group EA and group MA, including 19(S)-HETE, PGB2, and 1-Methylnicotinamide (Supplementary Table S6). These findings indicate that manual acupuncture results in a greater number of callback serum metabolites and that there is a substantial overlap in the metabolites affected by both EA and MA.

By intersecting the fecal metabolites that were decreased in the M group compared to the WKY group and increased in the EA group compared to the M group, as well as those that were increased in the M group compared to the WKY group and decreased in the EA group compared to the M group, we identified a total of 318 callback fecal metabolites associated with EA (Figures 10E,F; Supplementary Table S7). Similarly, 293 callback fecal metabolites were identified for MA (Figures 10G,H; Supplementary Table S8). Among these, 160 fecal metabolites, such as 6-Keto-prostaglandin F1a, Indole-3-Carboxaldehyde, 7-α,27-Dihydroxycholesterol, 5-Hydroxy-N-formylkynurenine, were common to both EA and MA (Supplementary Table S9). These findings indicate that EA results in a greater number of callback fecal metabolites and that both EA and MA share a subset of these metabolites.

The serum and fecal metabolites from the EA callback were independently subjected to KEGG pathway enrichment analysis. The metabolic pathways implicated in the callback of serum metabolites associated with hypertension included glycerophospholipid metabolism, linoleic acid metabolism, arachidonic acid metabolism, along with 11 additional pathways (Figure 11A; Table 3). In contrast, the primary metabolic pathways involved in the callback of fecal metabolites pertinent to hypertension encompassed steroid hormone biosynthesis, aldosterone-regulated sodium reabsorption, and dopaminergic synapses, among 7 other pathways (Figure 11B; Table 3).

Similarly, serum and fecal metabolites from the MA callback were analyzed separately for KEGG pathway enrichment. The metabolic pathways involved in the callback of serum metabolites related to hypertension included linoleic acid metabolism, arachidonic acid metabolism, arginine biosynthesis, and 10 additional pathways (Figure 11C; Table 3). The principal metabolic pathways involved in the callback of fecal metabolites associated with hypertension comprised tryptophan metabolism, primary bile acid biosynthesis, dopaminergic synapses, and 7 other pathways (Figure 11D; Table 3).