The experimental protocol and animal use were approved by the Animal Use and Management Ethics Committee at the Animal Experimentation Management Center of Inner Mongolia University of Nationalities (Lot No. NM-LL-2019-10-12-01). The experiment was designed and implemented in accordance with the National Institute of Health Guidelines for the Ethical Use of Animals. All rats were acclimatized for 1 week prior to the start of the study.

After acclimatization, the rats were assigned to a blank group, model group, or Sanwei sandalwood decoction group (i.e., treatment group) using the random number table method, with 10 rats per group. Animals in the blank group were injected intraperitoneally with saline once a week. To prepare the heart failure model, doxorubicin (saline preparation, 2 mg/mL) was administered intraperitoneally to the remaining rats at a dose of 1.5 mL/kg, once a week for 7 weeks. A combined dose of 21 mg/kg was used to establish the model, in accordance with the literature (Yuhong et al., 2018). After 7 weeks of administration, rats in the treatment group were gavaged with Sanwei sandalwood decoction at a dose of 2 g/kg for 30 consecutive days, while rats in the blank and model groups were given an equivalent volume of saline (Zhou, 2017).

Rat feces were collected 24 h after the last dose and transferred to lyophilization tubes, then snap frozen in liquid nitrogen and stored in a −80°C freezer (Shao et al., 2020).

The sample stored at −80°C refrigerator was thawed on ice. A 400 μL solution (methanol: water = 7:3, V/V) containing internal standard was added into 20 mg sample, and vortexed for 3 min. The sample was sonicated in an ice bath for 10 min and vortexed for 1 min, and then placed in −20°C for 30 min. The sample was then centrifuged at 12,000 rpm for 10 min (4°C). And the sediment was removed, then centrifuged the supernatant at 12,000 rpm for 3 min (4°C). A 200 μL aliquots of supernatant were transferred for LC-MS analysis.

The FASTP software was used for raw data quality control. Default software parameters were selected to preprocess the raw data obtained by Illumina HiSeq sequencing, and the resultant clean data was used for subsequent analysis (Karlsson et al., 2013). To remove host sequences from the samples, reads were filtered by comparison against the host database using Bowtie2 with the following parameter settings: --sensitive, −I 200, −x 400 (Karlsson et al., 2012; Scher et al., 2013).

After the above preprocessing steps, metagenome assembly was performed using the MEGAHIT software. The assembly parameters were as follows: --k-min 35, --k-max 95, --k-step 20, --min-contig-len 500. For each sample, the preprocessed clean data was compared against the assembled contigs to identify paired-end reads that were not utilized, using Bowtie2 with the following parameters: −I 200, −x 400 (Karlsson et al., 2013; Nielsen et al., 2014). Unutilized reads from each sample were combined for mixed assembly, with the same assembly parameters as for the single samples.

Based on the assembled contigs (≥500 bp) from initial and mixed assembly, MetaGeneMark was used with default parameters to predict ORFs (open reading frames), and after initial analysis, predicted genes shorter than 100 nt were filtered out (Mende et al., 2012; Li et al., 2014; Oh et al., 2014; Qin et al., 2014). The filtered ORF prediction results of all samples (including mixed assembly) were then combined, and CD-HIT (Li and Godzik, 2006; Fu et al., 2012) was used to remove redundancy and obtain a nonredundant initial gene catalog (with each gene defined by the longest nonredundant continuous nucleic acid sequence corresponding to that gene). The following parameter settings were used: −c 0.95, −g 0, −AS 0.9, −g 1, −d 0 (Li et al., 2014; Qin et al., 2014). By default, 95% identity and 90% coverage were used for clustering; and as mentioned above, the longest sequence for each gene was selected as the representative sequence. Next, Bowtie2 was used to compare the clean data for each sample against the initial gene catalogue (Li et al., 2014), and the number of reads per gene was calculated for each sample, using the following settings: --end-to-end, --sensitive, −I 200, −x 400. Genes with ≤2 reads in each sample were filtered out to obtain a final gene catalog (of Unigenes) that was used for subsequent analysis. Based on the number of reads and gene length, an abundance value was calculated for each gene in each sample, according to the following equation (Cotillard et al., 2013; Le Chatelier et al., 2013; Villar et al., 2015), where R is the number of reads compared to the gene, and L is the length of the gene.

To generate species annotations, the DIAMOND software (Buchfink et al., 2015) was employed to compare Unigene sequences against those of bacteria, fungi, archaea, and viruses from the NCBI NR database, using the lowest common ancestor (LCA) algorithm with the BLASTP search protocol and E-value ≤10⁻⁵. Note that because each sequence may give multiple comparison results, many different species classifications can be obtained. To ensure biological significance, the MEGAN software (Huson et al., 2011) was used again with the LCA algorithm to obtain final species annotation information for the sequences. Based on the results of LCA-based annotation and the gene abundance values, species abundance information was derived for each sample at each taxonomic level (i.e., phylum and genus), and a gene number table for each sample at each level (phylum compendium genus) was also obtained. The abundance of a species in an individual sample was taken as the sum of the annotated gene abundance of that species. The gene number of a species in a sample was defined as the number of genes annotated from that species, for all genes with nonzero abundance. In addition, the following procedures were performed based on the abundance tables at each classification level (genus/species of phylum): Krona analysis (Ondov et al., 2011); relative abundance overview; abundance clustering heat map analysis; principal component analysis (PCA); nonmetric multidimensional scaling (NMDS) dimensionality reduction; analysis of similarities (ANOSIM) inter-/intra-group difference analysis (Rao, 1961); and linear discriminant analysis effect size (LEfSe) multivariate statistical analysis of intergroup differential species.

All samples were acquired by the LC-MS system followed machine orders. The analytical conditions were as follows, UPLC: column, Waters ACQUITY UPLC HSS T3 C18 (1.8 μm, 2.1 mm × 100 mm); column temperature, 40°C; flow rate, 0.4 mL/min; injection volume, 2 μL; solvent system, water (0.1% formic acid): acetonitrile (0.1% formic acid); The column was eluted with 5% mobile phase B (0.1% formic acid in acetonitrile) at 0 min followed by a linear gradient to 90% mobile phase B (0.1% formic acid in acetonitrile) over 11 min, held for 1 min, and then come back to 5% mobile phase B within 0.1 min, held for 1.9 min.

The original data file acquisited by LC-MS was converted into mzML format by ProteoWizard software. Unsupervised PCA (principal component analysis) was performed by statistics function prcomp within R. ² The HCA (hierarchical cluster analysis) results of samples and metabolites were presented as heat maps with dendrograms. For two-group analysis, differential metabolites were determined by VIP (VIP >1) and p-value (p-value <0.05, Student’s t test). Identified metabolites and annotated metabolites were annotated and mapped using KEGG Compound and Pathway database (http://www.kegg.jp/kegg/compound/ and http://www.kegg.jp/kegg/pathway.html).

The data was processed using SPSS 22.0 and values are expressed as (x ̅ ± s). The t-test was used for pairwise comparisons between groups if the data met both normality and chi-square; if not, the rank sum test was used. One-way ANOVA was used to compare multiple groups. Differences where p < 0.05 were regarded as statistically significant.

Basic statistics on the gene catalog are presented in Table 1. Using the gene abundance information, core genome and pan-genome analyses were performed, along with intersample correlation analysis and Venn diagram construction. As depicted in Figure 1A, the core and pan-genome analyses revealed the number of genes after sample dilution. The sparse curve tends to flatten out and approach with increasing data, indicating that the collected samples were suitable for subsequent bioinformatics analysis. The number of genes differed between groups, and there was also variability between groups (Figure 1B). The Venn diagram (Figure 1C) shows that there were 778,808 common genes between the blank and model groups, while 49,315 genes were lost and 521,008 were gained in the model group relative to the blank group, suggesting that heart failure induces changes in the number of microbial genes found in feces. Compared with the model group, this value was also significantly affected by treatment with Sanwei sandalwood decoction. Herein, in the rat heart failure model, the number of genes specific to the intestinal flora increased after the therapeutic intervention. Moreover, we determined the correlation between samples to assess the degree of similarity between them, i.e., in terms of gene expression levels. As shown in Figure 1D, the correlation coefficients for all pairwise comparisons were <1, indicating a low degree of similarity and implying a high number of differential genes between experimental groups.

Misuse of antibiotics has induced irreversible changes in human and environmental microbial communities, posing new risks to human health and ecological systems. The study of antibiotic resistance genes (ARGs) has therefore received considerable attention from researchers (Guo et al., 2018). In recent years, the CARD resistance gene database has been developed to assist these efforts. It has the advantages of comprehensive information coverage, user-friendliness, and timely updates and maintenance (Jia et al., 2017). Here, it was used to analyze the abundance of AROs and produce a bar chart of abundance, a heat map of abundance clusters, and a map of abundance circles, as well as identifying ARO differences between groups. It also provided information on the species affiliation of resistance genes (i.e., Unigenes annotated to AROs). The details are shown in Supplementary Figure S1.

To visualize the overall distribution of ARO abundance in each group, the top 10 most abundant AROs were selected and plotted in an overview circle (Supplementary Figures S2, S3A). Additionally, the top 30 AROs were used to generate an abundance clustering heat map; and to test whether between-group differences were significantly greater than within-group differences, ANOSIM analysis of the abundance of resistance genes was performed (Supplementary Figure S3B). To further examine differences between sample groups, box plots were made of differences in the number of resistance genes and the number of AROs between groups (Supplementary Figure S3C).

The relative abundance of the genera Bacteroides and Prevotella, and the phylum Bacteroidetes, was significantly reduced in the intestinal flora in the rat model of heart failure, while the relative abundance of the genera Bacteroides, Clostridium, and Lactobacillus, and the order Actinomycetales, was markedly increased. After treatment with Sanwei sandalwood decoction, the relative abundance of the phylum Bacteroidota, and the genera Aspergillus and Prevotella, increased significantly in the intestinal flora, while that of Lactobacillus spp. decreased. The intervention also led to a significant reduction in the relative abundance of Aspergillus spp. and the phylum Bacteroidota, while balancing the relative abundance of Lactobacillus spp. and other genera that had been altered by heart failure. The relative abundance of Phyllobacterium spp. in the intestinal flora of rats after heart failure, and the least clustered abundance of these species, were significantly rescued after the treatment intervention; and the gut microbiota became more enriched in Phyllobacterium spp., Lactobacillus spp., and Bifidobacterium spp. These are all important probiotic species in the intestinal flora that play key roles in metabolic processes in vivo.

The metabolic functional classes had the highest relative abundance in the intestinal flora, among KEGG class I pathways. Moreover, carbohydrate transport and metabolism, and amino acid transport and metabolism, were more abundant in the gut microbiota of all three experimental groups. Of the six protein functional modules in the CAZy database, the classes GHs, GTs, and CBMs had higher relative abundance; while KO clustering analysis revealed that the Sanwei sandalwood decoction group had a greater degree of energy metabolism, carbohydrate metabolism, amino acid metabolism, and membrane transport compared with the model group. Based on eggNOG functional annotation, the processes of biosynthesis; transport and catabolism of secondary metabolites; energy production and conversion; transcription; amino acid transport and metabolism; signal transduction mechanisms; and nucleotide transport and metabolism were more abundant in the treatment group. The intestinal flora of rats in the dispersion group was significantly higher than that in the model group. In summary, the most abundant functional annotations in the intestinal flora of rats with heart failure, and treated with Sanwei sandalwood decoction, were mainly distributed among CBMs/carbohydrate metabolism, amino acid metabolism, nucleotide transport, energy metabolism, and GTs.

Metabolomics enables qualitative and quantitative assessment of metabolites in living organisms (Guo, 2021). Metabolites are downstream of gene transcription and protein translation and are closely related to phenotypes, especially for multifactorial diseases. In recent years, metabolomics has been increasingly used to identify heart failure biomarkers (Laíns et al., 2019). Compared with conventional biomarkers, metabolomics is more accurate in assessing heart failure-related metabolic disorders, and metabolite changes are more conducive to assessing the prognosis of heart failure (Cheng et al., 2015; Lanfear et al., 2017). Plasma-targeted surrogomics studies in heart failure patients with reduced ejection fraction identified 13 metabolites independently associated with survival, further demonstrating the promise of metabolic impairment and metabolomics in assessing heart failure phenotype and risk stratification and identifying therapeutic targets. In this study, PCA and OPLS-DA S-PLOT found that the three groups were clearly distinguished, while the position of the Sandalwood soup group tended to be the normal group, which could improve the metabolism of heart failure, and the classification of KEGG differential metabolites in enrichment analysis showed that Sandalwood soup may regulate heart failure by regulating 15 metabolites such as Loxoprofen, Conifery-l-acetate, and Trichilin A. It can be concluded that arachidonic acid metabolism pathway occupies an important proportion in the pathway classification chart, and it is known that arachidonic acid (AA) is an important long-chain hyperunsaturated fatty acid in organisms, and its metabolism plays an important role in the occurrence and development of myocardial fibrosis (Zhang et al., 2018). Experiments further show that the issue of the possible effect of Sanwei sandalwood decoction on the diversity of intestinal flora in the body by regulating differential metabolites, such as increasing the expression of AA, is also worth further exploration.

It is already established that ARGs change in response to variations in the microbiome (Zhao et al., 2021), and that reductions in ARGs generally correlate with reductions in the bacterial community. In particular, research on β-lactam antibiotic resistance genes (ARGs) is an important area of study that cannot be ignored. In this work, resistance genes against other antibiotic classes—such as tetracyclines, macrolides, and quinolones—were identified less frequently. However, in a finding that has not been previously reported, ARGs for some commonly used disinfectants and antiseptics (triC and qacH) were detected in this study.