SHL samples were purchased from Harbin Pharmaceutical Group Sanjing Pharmaceutical Co., Ltd. (China, batch number: 20081311). LPS was purchased from Shanghai Yuanye Bio-Technology Co., Ltd. (China, batch number: M11GS141189). Aspirin was purchased from Cisen Pharmaceutical Co., Ltd. (China, batch number: 2005112512). Cholic acid-2,2,3,4,4,-d5 (CAS: 53007-09-3; purity ≥98%) and succinic acid-2,2,3,3,-d4 (CAS: 14993-42-6; purity ≥99%) were purchased from Sigma Company (United States). DL-tryptophan-2,3,3-d3 (CAS: 340257-61-6; purity ≥98%) and DL-methionine-3,3,4,4,-d4 (CAS: 93709-61-6; purity ≥99%) were purchased from CDN Company (China). L-phenylalanine (Ring-D5) (CAS: 63-91-2; purity ≥98%) and choline chloride (Trimethyl-D9) (purity ≥98%) were purchased from CIL Company (United States). Agarose (batch number: 75510-019) was purchased from Invitrogen (United States, batch number: DL15000). The marker was purchased from Takara (Japan). TAE (batch number: AM9870) was purchased from Invitrogen (United States). Quant-iT PicoGreen dsDNA assay kit (batch number: P7589) was purchased from Invitrogen (United States). Q5® High-Fidelity DNA polymerase (lot number: M0491L) was purchased from Beijing NBE (China). Agilent High Sensitivity DNA kit and MiSeq Reagent Kit V3 (600 cycles, Illumina, United States) were used. Methanol and acetonitrile (purity ≥99.0%) were purchased from Thermo (United States). Ammonium formate (purity ≥99.0%) was purchased from Sigma (United States). Formic acid was of LC-MS grade, and other reagents were of analytical purity. All rat-specific enzyme-linked immunoassay (ELISA) kits were purchased from Multi Science (Lianke) Biotech, Co., Ltd., Hangzhou, China.

According to the original composition and preparation method of Shuang-Huang-Lian oral liquid recorded in the Chinese Pharmacopeia (The State Commission of Chinese Pharmacopoeia, 2020), it was prepared by the following procedure. 375 g of Scutellaria baicalensis Georgi (Lamiaceae; Radix Scutellariae), 375 g of Lonicera japonica Thunb. (Caprifoliaceae; Flos Lonicerae japonicae), and 750 g of Forsythia suspensa (Thunb.) Vahl. (Oleaceae; Fructus Forsythiae) were decocted by boiling for 2 h once and then boiling for 1 h twice, concentrated, extracted with ethanol, adjusted PH value with HCl and NaOH, distilled to eliminate solvent, and the residue was dissolved and diluted with water to 1000 ml in volume. The extracted solution was filtered and made to a concentration of 1.5 g raw material per milliliter. In this study, the concentration of SHL preparation (batch number: 20081311) was 1.5 g mL⁻¹. LPS and aspirin were weighed accurately and dissolved in normal saline into a 50 ml volumetric flask separately to obtain 100 μg mL⁻¹ LPS stock solution and 10 mg mL⁻¹ aspirin solution.

Healthy adult male Wistar rats (weight: 220 ± 20 g, purchased from Jinan Pengyue Experimental Animal Breeding Co., Ltd., Shandong, China) were adaptively fed for 3 days in standard laboratory conditions at Shandong University of Traditional Chinese Medicine. After passing the quarantine, they were randomly divided into six groups, namely the normal control group (CG), the model group (MG), positive drug group (PG), SHL low (LG), medium-dose group (ZG), and high-dose group (HG), with eight rats each group. For three consecutive days prior to the LPS intervention, an SHL oral dose of 6.3, 12.6, and 25.2 ml kg⁻¹ was administrated once a day to the LG, ZG, and HG groups, respectively. After fasting for 24 h, the LG, ZG, and HG groups were orally administered with SHL again. The CG and PG were injected with an equal volume of normal saline and aspirin (100 mg kg⁻¹) after 0.5 h. Then, LPS solution (100 μg kg⁻¹) was intraperitoneally injected into each group except CG. The body temperature level of rats in each group was detected every 30 min. Blood samples were collected from the inner canthus of the rats under ether anesthesia at different time points. Then, 300 mg of intestinal excreta were collected in the sterile tube and stored in liquid nitrogen immediately. Serum samples were placed at room temperature for 30 min and then centrifuged at 3000 rpm for 10 min. All samples were stored at −80°C before analysis.

Differences in gut microbiota among the CG, MG, and the optimal concentration of SHL administration group were assessed by collecting intestinal excreta samples from four randomly selected rats in each group for 16S rRNA gene analysis on the Illumina MiSeq platform. Total DNA was obtained using a DNA gel extraction kit. Furthermore, the purity and concentration of total DNA employed in this experiment were analyzed using a spectrophotometer, and DNA integrity was detected by 1.2% agarose gel electrophoresis. The V3–V4 region of the 16S rRNA gene was amplified by polymerase chain reaction (PCR) with the following primers: 338F (5′-ACT​CCT​ACG​GGA​GGC​AGC​A-3′) and 806R (5′-GGACTACHVGGGTWTCTAAT-3′). PCR assay was performed at 95°C for 3 min, followed by 27 cycles of 95°C for 30 s, 55°C for 30 s, and 72°C for 45 s, with the last step at 72°C for 10 min. The amplified PCR products were excised from 2% agarose gels and purified using the AxyPrep DNA Kit. The concentration of libraries was determined using a Qubit 2.0 fluorometer (Invitrogen) with the Quant-it PicoGreen dsDNA Assay Kit (Thermo Scientific).

According to the overlapping relationship between the paired-end reads, the sequence data achieved by HiSeq sequencing were merged into a sequence of tags, and the quality of reads was controlled and filtered. After discriminating samples, community bar charts were analyzed at phylum and genus levels. β-diversity was further studied by the Bray–Curtis distance algorithm and visualized by principal coordinate analysis (PCoA) (Bray and Curtis, 1957). Phylogenetic investigation of communities by reconstruction of unobserved states (PICRUSt) was used to predict the abundance of gene categories (COGs) and microbiota-derived pathways (Gavin et al., 2019). The description and functional information of each COG were parsed from the eggNOG database to obtain a functional abundance profile.

Internal standards solution includes succinic acid-2,2,3,3,-d4 (30 ppm final concentration), cholic acid-2,2,3,4,4-d5 (45 μM final concentration), L-phenylalanine-d5 (15 μM final concentration), DL-methionine-3,3,4,4-d4 (30 μM final concentration), DL-tryptophan-d3 (15 μM final concentration), and choline chloride-d9 (15 μM final concentration). Serum samples (100 μL) were added to methanol solution (400 μL) and the mixed solution of internal standards (100 μL). Then, the mixture was vortexed for 1 min with centrifugation for 10 min at 12,000 rpm, 4°C. The supernatant (500 μL) was dried under a vacuum concentrator. It was re-dissolved in 80% methanol solution (150 μL), thoroughly mixed and centrifuged (12000 rpm, 10 min, 4°C). Quality control (QC) samples were prepared to balance the chromatography–mass spectrometry for evaluating the stability of the system during the entire experiment. 20 µL of each extract supernatant was mixed to form a single QC sample for ultra-high-field X orbitrap mass spectrometry analysis (Sangster et al., 2006; Zelena et al., 2009; Dunn et al., 2011).

Serum samples were analyzed on the Waters Acquity UPLC System (Milford, MA, United States). The Acquity UPLC HSS T3 C18 column (2.1 mm × 150 mm, 1.8 μm; Waters, Milford, MA, United States) was selected. In positive mode, the mobile phases were as follows: A1: 0.1% aqueous formic acid solution and B1: 0.1% formic acid acetonitrile solution. The mobile phases in negative ion mode were A2: 5 mM ammonium formate water solution and B2: acetonitrile solution. The injection volume was 2 μL, and the flow rate was 250 μL/min. The column temperature was set to 40°C. The gradient elution procedure was as follows: (0–1 min), 2% B1/B2; (1–9 min), 2–50% B1/B2; (9–12 min), 50–98% B1/B2; (12–13.5 min), 98% B1/B2; (13.5–14 min), 98–2% B1/B2; (14–20 min), 2% B1-positive mode; (14–17 min), and 2% B2-negative ion mode. Ultra-high-field X Orbitrap mass spectrometer (Thermo, United States) and electrospray ionization (ESI) were used for mass spectrometry detection with positive and negative ion modes. The spray voltages of positive ions and negative ions were 3.50 and 2.50 kV, respectively. Other parameters were consistent, and the setting was as follows: sheath gas, 30 arb; auxiliary gas, 10 arb; capillary temperature, 325°C; resolution, 60000; and scan range, 81–1000 m/z. The secondary cracking was performed with HCD at a collision voltage of 30 eV, and unnecessary MS/MS information was removed by using the dynamic exclusion method.

ProteoWizard software was used to convert the high-resolution mass spectrometry data into mzXML format and R software (v3.3.2, XCMS kernel) was used for retention time correction, peak identification, peak extraction, and peak alignment. The metabolites were identified based on accurate mass (mass measurement accuracy within 10 ppm). After that, metabolites were obtained according to the fragmentation pattern of MS/MS and confirmed by Metlin (http://metlin.scripps.edu), MoNA (https://mona.fiehnlab.ucdavis.edu/), and a self-built standard database. The dataset was manually corrected to remove system contaminants and uninformative data by comparing with the blank sample and then normalized to sample weight before various statistical analyses were conducted. Differential metabolite analysis consisted of UPLC-MS data preprocessing and multivariate analysis such as principal component analysis (PCA) and orthogonal partial least squares discrimination analysis (OPLS-DA). PCA was used to get an overview of sample distribution and possible outliers. OPLS-DA was carried out to identify significant metabolites. The quality of the OPLS-DA model was assessed by a 7-fold cross-validation method. R²Y and Q² values were presented to show the ability of the model to explain y matrix information and the predictability of the model, which was calculated by SIMCA-P⁺ software (v14.1, Umetrics AB, Umea, Sweden). In order to further verify the fitting degree of the model, permutation tests were performed many times (n = 200) randomly to get the corresponding random Q² values. According to the results of OPLS-DA, variable importance projection (VIP) > 1 and p-value < 0.05 were considered potential chemical markers. Characteristic metabolites were further identified by the combination of secondary fragment ions and parent ions.

The serum concentrations of pyrogenic cytokines such as IL-6, interleukin-1β (IL-1β), and TNF-α were determined using rat-specific ELISA. In this study, all samples were tested in parallel with the secondary hole. For non-parametric distribution, the Kruskal–Wallis test complemented with Dunn’s tests was performed to determine statistical significance (*indicates p-value < 0.05).

To investigate the antifebrile effects of SHL, we recorded the rectal temperatures of rats continuously for 8 h after LPS infusion (Figure 1A, Supplementary Table S1). The results showed that the rectal temperature significantly increased in a three-phase thermal curve (peaks at ∼1, 2, and 5.5 h postinjection) after injecting LPS. In addition, ZG and HG could markedly decrease the temperature at almost all the checking points, indicating the antifebrile effect of SHL. However, HG exhibited lower body temperature than CG in the prefebrile phase. To study the optimal time point for blood sampling, concentrations of pro-inflammatory cytokines in CG, MG, and ZG at different blood collection time points were recorded. As shown in Figure 1B, TNF-α was detected within 30 min after LPS invasion and reached the maximum value earlier than the other two cytokines, which was consistent with other studies (Wang et al., 2019). All three pro-inflammatory cytokines reached high response values at 2nd h, and the contents of the three cytokines decreased significantly after SHL treatment (Supplementary Figure S1). That was why the optimal time point for 2 h-point blood sampling was selected for subsequent metabolomic analysis. As shown in Figure 1C, the non-parametric Mann–Whitney U test of the comparison between different groups found that compared with the normal control group, the levels of IL-1β ( ^## p < 0.01), IL-6 ( ^## p < 0.01), and TNF-α ( ^## p < 0.01) pro-inflammatory cytokines at 2 h time point were remarkably increased in rats with inflammation. And the levels of three pro-inflammatory cytokines decreased significantly after SHL treatment with medium doses (*p < 0.05). Thus, the optimal dose group (ZG) was selected for further serum metabonomics and 16S rRNA microbial community analysis.

Serum metabolomics can be used to find all the different metabolites among samples and the characteristic markers of biological systems by analyzing the composition and content changes of small-molecule chemicals (metabolites) in samples (Tuttolomondo et al., 2008; Zhang et al., 2012). In the relative quantitative analysis of metabolomics, the chemical components of CG, MG, and ZG were detected by ultra-high-field X Orbitrap mass spectrometer. Typical base peak intensity chromatograms of serum samples in positive and negative ion modes are illustrated in Supplementary Figure S2. Batch normalization was used to screen real biological signals, resulting in a data matrix containing 24 samples and 20,514 variables in positive ion mode. The data matrix containing 24 samples and 12,377 variables was obtained in negative ion mode. After further analysis of the filtered data, a total of 58 metabolites in negative ion and 127 metabolites in positive ion were identified, including amino acids, lipids, carbohydrates, nucleotides, xenobiotics, cofactors, and vitamin compounds.

In order to identify the metabolic characteristics of serum samples from each group, multivariate statistical methods were used to analyze the normalized metabolic data. Principal component analysis (PCA) and orthogonal projections to latent structures discriminant analysis (OPLS-DA) were used to explore the differential variables. PCA score plot showed that the serum datasets of different groups could be separated in positive ion mode (Supplementary Figure S3A) but not significantly distinguished in negative ion mode (Supplementary Figure S3B). No abnormal samples were found according to the DModX plot, after which OPLS-DA analysis was performed. The scores plots of PCA and OPLS-DA in positive and negative ion modes between CG and MG are shown in Figure 2. According to the PCA score plots (Figures 2A and B), there was an incomplete separation between CG and MG, and one outlier was found. Then, the OPLS-DA model was often used to analyze standardized data and identify marker compounds because of its ability to minimize the differences between groups. According to the OPLS-DA score plots (Figures 2C and D), the samples could be divided into two categories clearly, and the model established had no risk of overfitting in the 200 permutation tests (Supplementary Figures S4A and B). Furthermore, OPLS-DA combined with VIP value (threshold >1) and p-value (threshold <0.05) was often used to determine the potential chemical markers. Comparing the inflammatory MG with the CG, 31 metabolites showed significantly different levels in the positive ion mode and 8 in the negative ion mode (Table 1). Differentially expressed biomarkers mainly included lipid, carbohydrate, amino acid, nucleotide, pyrimidine, neurotransmitter, and vitamin compound. Changes in identified serum biomarkers associated with LPS-induced inflammatory hyperthermia were described in heatmaps (Figures 2G and H). Topological metabolic pathway analysis was performed for all 39 potential biomarkers using KEGG and MetPA databases. The metabolic network associated with the inflammatory hyperthermia model is shown in Table 1. For serum samples in rats with LPS-induced inflammation, the metabolic pathways were determined as D-glutamine and D-glutamate metabolism, riboflavin metabolism, cholesterol and arachidonic acid metabolism, arginine and proline metabolism, taurine and hypotaurine metabolism, butanoate metabolism, pyrimidine metabolism, alanine, aspartate, and glutamate metabolism, and tryptophan metabolism.

In order to explore the effect of SHL on the serum metabolomics of LPS-induced inflammation rats, OPLS-DA analysis was performed on data from all groups (CG, MG, and ZG) to obtain comprehensive metabolic characteristics. According to the OPLS-DA score charts of the three groups (Figures 3A–D), the ZG and the CG were significantly different from the MG on both sides of the Y-axis, indicating significant differences between the groups. In addition, in order to further verify the fitting degree of the model, permutation tests were performed many times (n = 200) randomly to get the corresponding random Q² values. The results showed that the established OPLS-DA model was reliable and had good predictability (Supplementary Figures S4C and D). After treatment with SHL, the eight endogenous metabolites disordered in rats with inflammatory hyperthermia returned to almost normal levels (Figures 3E and F). As shown in Figure 4 and Table 2, the non-parametric Mann–Whitney U test of the comparison between the groups found that the eight different metabolites had a significant callback trend after the administration, such as uric acid, riboflavin, pyrrolidonecarboxylic acid, N-acetylserotonin, cytidine, 4-hydroxycinnamic acid, 1-methylxanthine, and palmitoleic acid. These potential biomarkers were related to riboflavin metabolism, D-glutamine and D-glutamate metabolism, caffeine metabolism, ABC transporters, and vitamin rare earths and absorption. They were respectively attributed to coenzyme factor and vitamin metabolism, amino acid metabolism, exogenous biodegradation and metabolism, nucleotide metabolism, and other upstream metabolic pathways, so it is speculated as potential target pathways for SHL to exert antipyretic and anti-inflammatory effects.

In order to explore whether the antipyretic and anti-inflammatory effects of SHL are affected by the gut microbiota, the 16S rRNA gene sequence of fecal samples of rats in different groups was analyzed. At the phylum level, the global community structure of the sample showed that the proportion of Bacteroides and Firmicutes accounted for about 90%. Compared with the CG, the ratio of Firmicutes to Bacteroides in the MG increased, and the ratio of Proteobacteria to Actinobacteria decreased. After the treatment of SHL, the ratios all returned to the level of the CG, maintaining the steady state of the microbiota (Figure 5A). In addition, the alpha diversity of the samples was assessed at the OTU level. The index refers to the diversity within habitat or intra-community, and it is in direct proportion to species richness. The results showed that there had a high degree of microbial community richness, and the Kruskal–Wallis test showed no significant difference in community richness among groups (Supplementary Figure S5). Due to the low OTU abundance in the C1 sample and the large difference in the blank group, C1 was eliminated in subsequent analysis.

Principal coordinate analysis (PCOA) based on the unweighted UniFrac distance exhibited that there was a significant difference between CG and MG. After SHL treatment, the gut microflora structure of the MG was significantly changed. Principal components 1 (PC1) and PC2 accounted for 16.9 and 14.1% of the total variance respectively, explaining a total of 30% of the variation (Figure 5B).

Biomarker discovery associated with inflammation was conducted by linear discriminant analysis (LDA) effect size (LEfSe), which was used to determine the significant differences among different groups. As shown in Figure 5C and Supplementary Figure S6, compared with the CG, the diversity of community composition in the MG was mainly as follows: at the genus level, the abundance of Bifidobacterium, Olsenella, Lacticigenium, and Isobaculum in the MG was significantly higher than that in the CG (*p < 0.05) and the abundance of Bacteroides, Subdoligranulum, and Faecalibacterium decreased significantly (*p < 0.05). Bifidobacterium, Olsenella, Bacteroides, the unknown P_75_a5 species, Clostridium, Enterococcus, and Aerococcus levels were restored by treating with SHL. Meanwhile, the abundance of probiotics Roseburia was significantly increased (*p < 0.05), the pathogenic abundance of Turicibacter, Allobaculum, Jeotgalicoccus, Lactococcus, the unknown P_75_a5 species, Clostridium and Ralstonia decreased significantly (*p < 0.05) (Figure 5D, Supplementary Figure S7). In general, SHL can regulate the gut microbiota dysbiosis of LPS-induced inflammatory injury in rats.

These microbiota were ranked based on their variable importance in genus-level as follows: Olsenella > Clostridium > Enterococcus > Bifidobacterium > Aerococcus. MetagenomeSeq method was used to analyze the differential metabolic pathways among each group. According to the criteria of *p < 0.05, the differential metabolic pathways in the KEGG and MetaCyc pathways between the MG and the ZG were mainly reflected as follows: the Geraniol degradation pathway was significantly upregulated, and the proteasome and naphthalene degradation pathways were significantly downregulated in the ZG (Figure 5E).

Taking the content of cytokines IL-6, TNF-α, and IL-1β as the X variable, and the peak area of the different metabolites in each group as the Y variable, the data matrix was performed to z-score standardization for further analysis. Spearman correlation analysis was used to calculate the XY data matrix correlation and p-value, and the specific results are shown in Supplementary Table S2. N-acetylserotonin was significantly negatively correlated with IL-6 and IL-1β (*p < 0.05). 4-hydroxycinnamic acid was significantly negatively correlated with IL-6, TNF-α, and IL-1β (**p < 0.01, ***p < 0.001, and *p < 0.05). Uric acid was significantly positively correlated with IL-6 and TNF-α (*p < 0.05). N-acetylserotonin and 4-hydroxycinnamic acid showed a significant negative correlation with all pro-inflammatory cytokines. These two endogenous metabolites are significantly decreased in the MG (compared to the CG), resulting in increased pro-inflammatory cytokines and heat transfer in MG. In addition, 1-methylxanthine and uric acid had a significant positive correlation with IL-6 and TNF-α (*p < 0.05, **p < 0.01). However, the levels of TNF-α were significantly positively correlated with cytidine (**p < 0.01). Therefore, it is speculated that 4-hydroxycinnamic acid, uric acid, N-acetylserotonin, 1-methylxanthine, and cytidine are biomarkers exhibiting antipyretic and anti-inflammatory therapeutic effects of SHL. In general, the current analysis found strong evidence of a correlation between serum biomarkers and pro-inflammatory cytokines.

In order to comprehensively analyze the relationship between eight biomarkers and 12 kinds of gut microbiota, a correlation matrix was established by calculating the Spearman correlation coefficient, as shown in Supplementary Table S3. The results showed that pyrrolidone carboxylic acid was significantly positively correlated with Actinobacteria, Bifidobacteriumles, Bifidobacteriumceae, Enterococcaceae, and Bifidobacterium (*p < 0.05), and significantly negatively correlated with Bacteroidaceae and Bacteroides (*p < 0.05). N-acetylserotonin was significantly negatively correlated with Actinobacteria, Bifidobacteriumles, Bifidobacteriumceae, Enterococcaceae, Bifidobacterium, Olsenella, Enterococcus (*p < 0.05). Cytidine and Olsenella presented significant negative correlations (*p < 0.05). 1-methylxanthine was significantly positively correlated with Bifidobacteriumles, Bifidobacteriumceae, and Bifidobacterium (*p < 0.05) and significantly negatively correlated with Bacteroidaceae, Olsenella, Bacteroides, and Clostridium (*p < 0.05). Observations such as these underscored the close association between the metabolites and gut microbiota.