Microbiome and metabolomics analyses of the effect of heat-sensitive moxibustion on allergic rhinitis in rats

Heat-sensitive moxibustion was purchased from the Affiliated Hospital of Jiangxi University of Traditional Chinese Medicine (Nanchang, China); ovalbumin (OVA; American Sigma, St. Louis, Missouri) and aluminum hydroxide (Tianjin Damao, Tianjin, China) were used. Plasma was used to determine IL-4 and IgE levels using an enzyme-linked immunosorbent assay (ELISA) kit (Shenzhen Xinbosheng); nasal mucosa was used to determine histamine levels using an ELISA kit (Jiangsu Meimian, Jiangsu, China).

Thirty-six male Sprague–Dawley (SD) rats (5 weeks old, 151–175 g) were purchased from Hunan Sheng Hua Laboratory Animal Co., Ltd. (Hunan Province, China). The rats were housed at the Experimental Animal Center of Jiangxi University of Chinese Medicine under controlled conditions: 12-hour light/dark cycles, 22 °C–26 °C room temperature, 50%–65% humidity, and adequate ventilation. All procedures were approved by Jiangxi Medical University's Institutional Animal Care and Use Committee (IACUC; Approval No. JZLLSC20240539) and conducted in accordance with institutional guidelines.

A total of 36 rats were randomly divided into two groups: blank group (n = 9) and OVA group (n = 27) (16). According to the previous literature, an AR model was established (17). During the sensitization phase (days 1–14), the control group received standard feed without interventions. The model group received intraperitoneal injections of sensitizing solution [0.3 mg ovalbumin (OVA) + 30 mg aluminum hydroxide in 1 mL saline] every other day for seven doses. On day 15, the model group was further sensitized with daily intranasal administration of 25 μg 10% OVA for seven consecutive days. Of 27 rats subjected to OVA sensitization, one died due to improper injection and one from nasal hemorrhage during modeling. Post-treatment symptom scores (sneezing, nasal scratching, and discharge) were assessed by blinded observers. A total score >5 indicated successful AR modeling. Rats were randomized into a model group (n = 8) and a moxibustion treatment group (n = 17).

To facilitate moxibustion, rats were gently restrained in a ventilated cage that allowed natural posture while minimizing movement. The room temperature was kept at 25.2 °C during the whole experiment. Lung acupuncture point (BL13), which is very important for AR, is the commonly used acupuncture point for AR clinical treatment and rat experiments (18, 19). Moxibustion was applied to depilated skin at a height of 3 cm using animal-specific moxa sticks (12-cm length, 0.6-cm diameter; produced by the Affiliated Hospital of Jiangxi University of Chinese Medicine). Both the heat-sensitive and standard moxibustion groups received daily 40-min sessions at the Feishu acupoint (BL13) for three consecutive weeks. During this period, the moxa stick was fixed onto the moxibustion stand.

The mid-tail temperature of rats was measured using a digital thermometer (Shanghai Medical Equipment Factory, Shanghai, China). Measurements were taken in a quiet environment with controlled room temperature (25 °C ± 2 °C). Before the experiment, rats were acclimated in their cages for 30 min. During the treatment phase, all three groups received ongoing sensitization with 5% OVA solution every 5 days. After heat-sensitive moxibustion treatment, one rat was excluded because its tail temperature response was indeterminate and deemed to be at the threshold for defining heat-sensitive responsiveness, ensuring the purity of the HM group phenotype. Based on tail temperature changes, the moxibustion group was subdivided into two subgroups—a non-elevated subgroup (≤1 °C, OM group, n = 8) and an elevated subgroup (>1 °C, HM group, n = 8 (16))—determined by average measurements over 7 days. The experimental protocol for OVA-induced AR and heat-sensitive moxibustion is summarized in Figure 1. After 40 min of treatment, there were some rats exhibiting an increase in tail temperature (TTI) (more than 2°C on average) in the moxibustion group (Table 1). The occurrence rate of TTI is eight out of 17 (47.05%). Furthermore, the temperature began to rise within approximately 10 to 15 min and then decreased within 30 to 40 min, lasting for approximately 20 to 25 min each time. There was no change in the tail temperature of the model group without treatment.

Symptom scores after the final treatment were assessed by observation of sneezing, nose scratching, and runny nose in rats and were unknown to the observer. A total score of more than 5 indicates that the AR model is successfully established, and the rats fluctuating around 1 °C were not included in the observation objects. After treatment in the moxibustion group, one rat was excluded because it could not be determined whether its tail temperatures were heat-sensitized.

After overnight fasting, 12-hour urine samples were collected from individually housed rats in metabolic cages (8:00 PM to 8:00 AM). Urine was centrifuged at 4 °C and 800 × g for 10 min to remove debris, and the supernatant was retained. After centrifugation at 12,000 ×g and 4 °C for 10 min, the clarified middle layer was aspirated, aliquoted into 1.5-mL tubes (1 mL per aliquot), and stored at −80 °C. Prior to dissection, rats were restrained by securing the neck, spine, and tail. Fecal samples were collected via stress-induced defecation onto sterile paper (replaced for each rat) and transferred to sterile microcentrifuge tubes using disinfected toothpicks. Fecal pellets (one to two per rat) were placed in sterile tubes, flash-frozen in liquid nitrogen for 5 min, and stored at −80 °C. Rats were euthanized via intraperitoneal injection of xylazine (200 μg/rat) and ketamine (2 mg/rat). Blood was collected from the brachial artery, centrifuged to isolate serum, and stored at −80 °C for ELISA analysis. The nasal septa of the rats were collected and placed in the fixative solution.

Blood was collected in a tube containing Ethylenediaminetetraacetic acid (EDTA) and centrifuged (3,500 × g, 4 °C, 10 min) to separate the blood. IL-4 and IgE levels were measured in plasma using an ELISA kit.

Mouse nasal septa were fixed in 4% paraformaldehyde, followed by decalcification in EDTA solution for 9 days, and embedded in paraffin. Tissue sections (5-μm thickness) were prepared and stained with hematoxylin and eosin (H&E).

Total microbial genomic DNA was extracted from rat dropping samples using the FastPure Stool DNA Isolation Kit (MJYH, Shanghai, China) according to the manufacturer's instructions. The quality and concentration of DNA were determined by 1.0% agarose gel electrophoresis and a NanoDrop2000 spectrophotometer (Thermo Fisher Scientific, Waltham, Massachusetts, USA) and kept at −80 °C prior to further use. The hypervariable regions V3–V4 of the bacterial 16S rRNA gene were amplified with primer pairs 338F (ACTCCTACGGGAGGCAGCAG) and 806R (GGACTACHVGGGTWTCTAAT) using the T100 Thermal Cycler PCR thermocycler (BIO-RAD, Hercules, California, USA). The PCR mixture included 4 μL 5 × Fast Pfu buffer, 2 μL 2.5 mM dNTPs, 0.8 μL each primer (5 μM), 0.4 μL Fast Pfu polymerase, 10 ng of template DNA, and ddH₂O to a final volume of 20 µL. PCR amplification cycling conditions were as follows: initial denaturation at 95 °C for 3 min, followed by 27 cycles of denaturing at 95 °C for 30 s, annealing at 55 °C for 30 s, extension at 72 °C for 45 s, and a single extension at 72 °C for 10 min, and ended at 4 °C. The PCR product was extracted from 2% agarose gel and purified using the PCR Clean-Up Kit (YuHua, Shanghai, China) according to the manufacturer's instructions and quantified using Qubit 4.0 (Thermo Fisher Scientific, USA).

Urine sample (200 μL) was transferred to a 1.5-mL centrifuge tube, and 400 μL of extract (acetonitrile:methanol = 1:1) was added. After vortex mixing for 30 s, cold ultrasonic extraction was performed for 30 min (5 °C, 40 kHz), and the sample was allowed to stand at −20 °C for 30 min and then at 4 °C. After centrifugation at 13,000 × g for 15 min, the supernatant was removed and dried using nitrogen, and 100 µL of compound solution (acetonitrile:water = 1:1) was reconstituted. Low-temperature ultrasonic extraction for 5 min (5 °C, 40 kHz) and centrifugation at 13,000 × g for 5 min were performed to remove the supernatant, which was transferred to a vial. In addition, 20 µL of the supernatant from each sample was used as a quality control sample to investigate the repeatability of the whole analysis process.

The instrument platform used for this Liquid Chromatography-Mass Spectrometry (LC–MS) analysis was Thermo Scientific ultra-high performance liquid chromatography–tandem Fourier transform mass spectrometry UHPLC-Q Exactive HF-X system. Chromatographic conditions were as follows: 2-μL samples were separated using an HSS T3 column (100 mm × 2.1 mm i.d., 1.8 µm) and then detected using mass spectrometry. Mobile phase A was 95% water + 5% acetonitrile (containing 0.1% formic acid), and mobile phase B was 47.5% acetonitrile + 47.5% isopropyl alcohol + 5% water (containing 0.1% formic acid). Separation gradient was performed as follows: at 0–3.5 min, mobile phase B increased from 0% to 24.5%, and the flow rate was 0.40 mL/min. At 3.5–5 min, mobile phase B increased from 24.5% to 65%, and the flow rate was 0.4 mL/min. At 5–5.5 min, mobile phase B increased from 65% to 100%, and the flow rate was 0.4 mL/min. At 5.5–7.4 min, mobile phase B was maintained at 100%, and the flow rate increased from 0.4 to 0.6 mL/min. At 7.4–7.6 min, mobile phase B decreased from 100% to 51.5%, and the flow rate was 0.6 mL/min. At 7.6–7.8 min, the mobile phase B decreased from 51.5% to 0%, and the flow rate decreased from 0.6 to 0.5 mL/min. At 7.8–9 min, the mobile phase B was maintained at 0%, and the flow rate decreased from 0.5 to 0.4 mL/min. At 9–10 min, 0% mobile phase B was maintained, and the flow rate was 0.4 mL/min. The column temperature was 40 °C.

Mass spectrum conditions were as follows: the sample quality spectrum signal was collected in positive and negative ion scanning modes, and the mass scanning range was 70–1,050 m/z. The flow rate of sheath gas was 50 psi, the flow rate of auxiliary gas was 13 psi, the heating temperature of auxiliary gas was 425 °C, the positive mode ion spray voltage was set to 3,500 V, the negative mode ion spray voltage was set to −3,500 V, the ion transport tube temperature was 325 °C, and the normalized collision energy was 20–40–60 V cyclic collision energy. The primary mass spectral resolution was 60,000, the secondary mass spectral resolution was 7,500, and the data were collected using Data-Dependent Acquisition (DDA) mode.

All data analyses were performed on the Megge BioCloud Platform (https://cloud.majorbio.com↗). The details were as follows: the mothur (20) software (http://www.mothur.org/wiki/Calculators↗) was used to calculate alpha diversity, which was measured using Chao 1 and Shannon index. Alpha diversity was analyzed using the Wilcoxon rank sum test. Principal coordinates analysis (PCoA) (primary coordinates) was used in combination with the base-based Bray–Curtis distance algorithm, and it was combined with the Permutational Multivariate ANalysis Of Variance (PERMANOVA) non-parametric test to determine whether the difference in microbial community structure between sample groups was significant. Linear discriminant analysis (LDA) effect size (LEfSe) analysis (21) (http://huttenhower.sph.harvard.edu/LEfSe↗ (LDA > 2, p < 0.05) was used to identify bacterial taxa that significantly varied in phylum to genus level abundance between different groups. Network plot analysis (22) was performed to pick species based on Spearman's correlation |r| > 0.6, p < 0.05.

After computations were completed, LC–MS raw data were imported into the metabolomics processing software Progenesis QI (Waters Corporation, Milford, MA, USA) for baseline filtering, peak identification, integration, retention time correction, and peak alignment. Finally, a data matrix of retention time, mass–charge ratio, and peak intensity was made. The MS and Tandem Mass Spectrometry (MSMS) mass spectrum data were also matched with those in the metabolomics public databases HMDB (http://www.hmdb.ca/↗) and Metlin (https://metlin.scripps.edu/↗) and the Meiji self-built library to obtain the metabolite information. After searching, the matrix data were uploaded to Meiji Biological Cloud Platform (https://cloud.majorbio.com↗) for data analysis. First, for data preprocessing, a data matrix with 80% rules was used to remove the missing value, namely, to keep at least a set of non-zero values of variables, and then the vacancy value (the smallest value in the original matrix) was filled. To reduce the error of the sample preparation and instrument instability, the sample mass spectrum peak response intensity was normalized using the sum normalization method to obtain the normalized data matrix. Variables with relative standard deviation (RSD) >30% were also removed, and log10 logarithm was used to obtain the final data matrix for subsequent analyses.

Differential analysis was performed on the preprocessed matrix files. The R software package ropls (Version 1.6.2) was used to perform principal component analysis (PCA) and orthogonal partial least squares discriminant analysis (OPLS-DA), and a seven-cycle interaction validation was used to evaluate the stability of the model. Furthermore, a Student's t-test was performed.

Differences in metabolites were distinguished using the KEGG database (https://www.kegg.jp/kegg/pathway.html↗) on the metabolic pathways of annotation, as well as the differences in metabolites involved in pathways. The Python software package scipy.stats was used for pathway enrichment analysis, and the biological pathways most relevant to experimental treatment were obtained using Fisher's exact test.

The behavioral scores of the rats were recorded within 30 min after modeling. The behavioral scores were divided into three symptoms: nose scratching, sneezing, and runny nose. The scores were assigned according to the behavioral score table.

The rat model of AR induced by OVA exhibited clinical symptoms similar to those of AR, such as sneezing, runny nose, and nasal itching. OVA stimulation significantly increased the incidence of sneezing and friction in rats (p < 0.01), and the modeling was significantly successful (Figure 2A). In the OVA group, epithelial cells were disordered, and submucosal eosinophils infiltrated. In the moxibustion group, inflammatory cell infiltration was observed in the mucous lamina propria, epithelial cell arrangement was disordered to a certain extent, and submucosal eosinophils infiltrated. The structure of the heat-sensitive moxibustion group was clear, the mucosal surface cilia were dense, the lamina propria had slight capillary dilatation, and a small amount of eosinophils accumulated, which was significantly reduced compared with that in the model group. The above indicates that heat-sensitive moxibustion is more effective than ordinary moxibustion in reducing nasal symptoms (Figures 2D–G).

We measured serum IgE and IL-4 levels in different groups of mice. Serum IgE and IL-4 in the AR group increased significantly compared with those in the Con group. After moxibustion treatment, both of them decreased significantly (Figures 2B, C).

To determine the taxonomic composition of the gut microbiota and differences in microbial diversity, alpha- and beta-diversity analyses were performed. Alpha diversity was used to analyze the richness (number of taxa) and evenness (abundance distribution of groups) within the same microbial community. Ace, Chao, Shannon, and Simpson indices were assessed for each sample. The Ace index, analyzed using the Wilcoxon test, showed a significant difference between the HM and AR groups and no difference between the OM and AR groups. The Chao index showed a significant difference between the NM and AR groups and no difference between the OM and AR groups. There was no difference in the Shannon index and Simpson index between the HM group and AR group (Figures 3A–D). Next, beta-diversity analysis was used to generate a weighted UniFrac PCoA and showed the clusters between the HM, OM, and AR groups. As shown in Figure 3E, the HM group and AR group showed different trends, compared with the OM group, which showed a similar trend.

Considering the different composition of the microbiota between the HM and OM groups, to further identify the unique potential mechanism of heat-sensitive moxibustion for AR, we performed untargeted metabolomics determination of the three groups using UHPLC-Q Exactive HF-X. The RSD distribution plot demonstrates that the analysis process has excellent reproducibility and that the data quality is reliable (Figure 5A).

PLS-DA is a supervised discriminant analysis method. It was found that the PLS-DA plot showed a clear separation of the HM sample from the AR sample, with a slight trend of dissociation from the AR group (Figure 5B). To further screen out the significantly different metabolites between the two comparison groups, two-group comparisons between the HM and OM groups were also performed. In the positive ion mode, the PLS-DA plot showed a clear separation of the HM samples from the OM samples, a trend also reflected in the negative ion mode. The goodness of fit was satisfactory (Figures 5C, D). In the positive ion mode, a total of 32 differential metabolites were identified in the HM and OM groups, with 23 metabolites upregulated and nine metabolites downregulated (Figure 5F). A total of 718 differential metabolites were identified in the HM and OM groups, with 358 metabolites upregulated and 360 downregulated (Figure 5H); 34 were identified between the HM and OM groups, including 24 upregulated and 10 downregulated (Figure 5G). A total of 449 differential metabolites were identified between the HM and AR groups, with 253 upregulated and 196 downregulated (Figure 5I). To visually demonstrate common or unique metabolites among individual metabolic sets, the results are shown using Venn diagram analysis: there were 19 potential biomarkers between the HM group and AR group, and between the HM group and OM group, which play potentially important roles in the heat-sensitive moxibustion group, common moxibustion group, and model group (Figure 5E). These potential biomarkers were then used to generate the predictive receiver operating characteristic (ROC) curves (Figure 5J). All area under the curve (AUC) values of the biomarkers were greater than 0.70, indicating that these 19 metabolites could be considered as potential markers.

The results of hierarchical clustering analysis are shown in Figures 5K, L; by clustering the differentially expressed metabolites, we can intuitively see the change trend of different metabolites in different groups. The metabolites clustered in the same cluster have similar expression patterns, may have similar functions, or participate in the same metabolic processes or cellular pathways. The color in the figure indicates the relative expression of metabolites in this group. Compared with the AR group, the HM group had significantly changed metabolites, including citric acid and 16β-hydroxyestradiol; compared with the OM group, the HM group also had changed metabolites, including fructosyl-lysine, l-hypoglycin A, and cAMP.

We performed KEGG pathway enrichment analysis for differential metabolites. The results showed that the metabolites of MH and OM were mainly involved in the MAPK signaling pathway, glycine, serine, and threonine metabolism, histidine metabolism, etc. (Figure 5N). The metabolites of MH and AR were mainly involved in nucleotide metabolism, tryptophan metabolism, and linoleic acid metabolism. After screening, we found that there were seven intersection pathways of HM vs. AR and HM vs. OM: glutamatergic synapse, dopaminergic synapse, Kaposi sarcoma-associated herpesvirus infection, cocaine addiction, melanogenesis, alcoholism, and histidine metabolism (Figure 5M).

This finding suggests that the histidine metabolism pathway is a core mechanism through which moxibustion exerts its anti-allergic effects and that heat-sensitive moxibustion (HM) and conventional moxibustion (OM) may differentially regulate this pathway. To clarify the activity changes in this pathway, we measured the levels of histamine—a key downstream effector molecule—using ELISA. As shown in Figure 6, histamine concentration in the nasal mucosa was significantly elevated in the AR model group (AR) compared to the blank control group (Con) (p < 0.0001). Following intervention with either HM or OM, histamine levels decreased markedly compared to those in the AR group (p < 0.0001). Compared with OM, the histamine concentration in HM was significantly decreased (p < 0.05).

In order to study the correlation between changes in metabolite levels and microbiome composition, Spearman's correlations were performed on the differentiated species sets and differentiated metabolic sets screened from the microbiome (Figure 7). Multiple correlation analysis showed that Erysipelotrichaceae, Saccharimonadaceae, and Lachnospiraceae were positively correlated with norvaline, quisqualic acid, Asp Leu Ser Glu, and other metabolites. They were negatively correlated with daidzein, 9alpha-(3-methyl-2E-pentenoyloxy)-4S-hydroxy-10(14)-oplopen-3-one and "1,2,3-trihydroxybenzene". In addition, Lactobacillaceae were positively correlated with ethyl glucuronide, zileuton O-glucuronide, trichloroethanol glucuronide, and other metabolites. They were negatively correlated with 9alpha-(3-methyl-2E-pentenoyloxy)-4S-hydroxy-10(14)-oplopen-3-one; Prevotellaceae were negatively correlated with daidzein. These results show that there is a close correlation between the nine microorganisms selected and the 19 metabolites selected. Heat-sensitive moxibustion can regulate the disturbance of microorganisms and metabolites in AR treatment.