Twenty-one-day-old female Sprague–Dawley (SD) rats (specific pathogen-free (SPF) grade) were raised in the Animal Experiment Center of Guang’anmen Hospital, Chinese Academy of Chinese Medical Sciences, and their experimental animal license number is SCXK (Beijing) 2019-0008. The feeding conditions were as follows: room temperature, 22°C ± 2°C; relative humidity, 50% ± 10%; and 12-h light/12-h darkness, alternating cycles. All rats drank and ate freely. There were two kinds of feed, namely, ordinary feed (energy supply ratio: protein 22.47%, fat 12.11%, and carbohydrate 65.42%) and high-fat feed (energy supply ratio: protein 19.5%, fat 23.8%, and carbohydrate 56.7%); the main ingredients are casein, starch, maltodextrin, sucrose, fiber vegetarian, soybean oil, etc. The rats and feed were provided by Beijing Huafukang Biotechnology Co., Ltd. This experiment was reviewed by the Experimental Animal Ethics Committee of the China Academy of Chinese Medicine (Approval Number: IACUC-GAMH-2017-002).

Modified Banxia Xiexin Decoction (ingredients: Rhizoma Pinelliae (Qing Banxia) 9 g, Scutellaria (Huanqin) 20 g, Rhizoma Coptidis (Huanlian) 10 g, Rhizoma Zingiberis (Ganjiang) 9 g, Codonopsis pilosula (Dangshen) 12 g, Licorice (Gancao) 12 g, Fructus Ziziphi Jujubae (Dazao) 9 g, Epimedium (Xianlingpi) 15 g, and Fruit of Chinese wolfberry (Gouqizi) 30 g. All Chinese herbs were supplied by the pharmacy of Guang’anmen Hospital, Chinese Academy of Chinese Medical Sciences, and then made into Banxia Xiexin Tang extract (3.9 g of crude drug per ml of extract); the extract was stored in a refrigerator at 4°C for later use.

Seventeen rats were randomly allocated to 3 groups: control group (group Z, n = 6), model group (group M, n = 6), and treatment group (group A, n = 5). To induce the PCOS-IR model, rats in the model group and treatment group were gavaged with 2 ml LET-CMC solution (1 mg/kg/day of letrozole dissolved in 1% carboxymethyl cellulose), while the control group was gavaged with 2 ml of normal saline, lasting 3 weeks. Then group A was given modified Banxia Xiexin Decoction (26.46 g/kg/day) for 2 weeks; at the same time, two other groups were given the same volume of normal saline. The dose was calculated according to the “equivalent dose conversion of body surface area,” (Xu et al., 2002) and it is administered at 2 times the clinically equivalent dose.

The morning after 2 weeks of treatment, blood samples were collected from the abdominal aorta of anesthetized rats using a procoagulation tube and centrifuged at 3,000 rpm for 20 min using the 4°C centrifuge (ICT15RE, Hitachi, Tokyo, Japan). Then the upper supernatants were carefully transferred to the EP tube using a pipette gun and immediately stored at a −80°C refrigerator for further detection. The levels of estradiol (E2, JEB-13691, Jinyibo, Wuxi, China), testosterone (T, ab10866, Abcam, Cambridge, MA, USA), LH (JEB-13706, Jinyibo, China), follicle-stimulating hormone (FSH; JEB-13680, Jinyibo, China), fasting insulin (FINS, ab100578, Abcam, USA), LPS (JEB-13984, Jinyibo, China), and tumor necrosis factor-α (TNF-α, ERC102a.96, Xinbosheng, China) were analyzed by ELISA kits, which were operated in strict accordance with each instruction manual. Blood glucose (Glu) was measured using Roche glucometer before anesthesia. The homeostasis model assessment of IR index (HOMA-IR) and HOMA of beta-cell function index (HOMA-beta) was calculated using the following formulas:

Rats were euthanized by cervical dislocation after anesthesia. The left ovaries of rats were collected, then fixed with 4% paraformaldehyde solution, embedded in conventional paraffin (EG140H, Leica, Wetzlar, Germany), sectioned at a thickness of 5 μm (RM2135, Leica), and processed according to the standard procedure of H&E staining, and ovary morphological changes were observed using a light microscope (Olympus, Tokyo, Japan).

Fecal samples were collected by rectal pressing method before the rats were anesthetized; 3–5 fecal samples were collected from each rat and immediately stored at −80°C.

Microbial community genomic DNA was extracted from rat fecal samples using the E.Z.N.A.^® soil DNA Kit (Omega Bio-Tek, Norcross, GA, USA) according to the manufacturer’s instructions. The DNA extract was checked on 1% agarose gel, and DNA concentration and purity were determined with NanoDrop 2000 UV-vis spectrophotometer (Thermo Scientific, Wilmington, DE, USA). The hypervariable region V3–V4 of the bacterial 16S rRNA gene was amplified with primer pairs 338F (5′-ACTCCTACGGGAGGCAGCAG-3′) and 806R (5′-GGACTACHVGGGTWTCTAAT-3′) by an ABI GeneAmp^® 9700 PCR thermocycler (ABI, Los Angeles, CA, USA). The PCR amplification of 16S rRNA gene was performed, 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, single extension at 72°C for 10 min, and end at 4°C.

Purified amplicons were pooled in equimolar and paired-end sequenced on an Illumina MiSeq PE300 platform (Illumina, San Diego, CA, USA) according to the standard protocols by Majorbio Bio-Pharm Technology Co. Ltd. (Shanghai, China). The raw 16S rRNA gene sequencing reads were demultiplexed, quality-filtered by fastp version 0.20.0 (Chen et al., 2018), and merged by FLASH version 1.2.7 (Magoč and Salzberg, 2011). Operational taxonomic units (OTUs) were clustered by UPARSE version 7.1 with ≥97% similarity (Edgar, 2013), and chimeric sequences were identified and removed. The taxonomy of each OTU representative sequence was analyzed by RDP Classifier version 2.2 (Wang et al., 2007) against the 16S rRNA database (silva138/16s_bacteria) using a confidence threshold of 0.7. In order to observe the α-diversity of the intestinal microbial community, including richness and diversity, the number of Ace and Shannon indexes was calculated by the built-in commands of Mothur (https://www.mothur.org/wiki/Download_mothur). Principal coordinate analysis (PCoA), based on weighted UniFrac distance, and partial least squares discriminant analysis (PLS-DA) were used to reflect the β-diversity of intestinal microbiota, which mainly express the significance of intestinal microbiota difference among groups. To obtain significantly different bacterial species among groups, linear discriminant analysis (LDA) effect size (LEfSe) was applied. Here, LDA values >4.0 were set as the thresholds for biomarker identification. Meanwhile, in order to predict the functional profile of the intestinal microbiota community, the Phylogenetic Investigation of Communities by Reconstruction of Unobserved States (PICRUSt, http://picrust.github.io/picrust/) was applied to perform microbiota functional prediction. After OTU sequences were clustered according to 97% similarity, the OTU abundance was standardized by PICRUSt, which meant that the influence of copy number of 16S marker gene in species genome was removed, and the information of Clusters of Orthologous Groups of proteins (COG) family and Kyoto Encyclopedia of Genes and Genomes (KEGG) corresponding to each OTU was obtained through the Greengene ID (http://greengenes.secondgenome.com/) and then according to the different databases, i.e., EggNOG (evolutionary genealogy of genes: Non-supervised Orthologous Groups, http://eggnog.embl.de/) and KEGG (http://www.genome.jp/kegg/); the COG function and each hierarchy level of KEGG pathways could be obtained.

SPSS statistical package 22.0 and GraphPad Prism version 8.3.0 were used for the statistical analysis of clinical data, and all data were expressed as means ± SD or as medians with interquartile ranges. The Kolmogorov–Smirnov test of normality was applied to all data. One-way ANOVA was used to evaluate the statistical significance of differences for normally distributed data. Then homogeneity of variances was tested, and if the variances were homogeneous, Scheffe’s test was used; if not, Dunnett’s method was used. For non-normally distributed data, the Mann–Whitney U-test (two-group comparison) and the Kruskal–Wallis test (more than two groups) can be used.

For bioinformatics analysis, Wilcoxon rank-sum test was used to analyze the differences in α-diversity index (including Ace and Shannon indexes) and the abundance of intestinal microbiota community between two groups. To analyze the significant difference in intestinal microbiota community (i.e., β-diversity) among groups, analysis of similarity (ANOSIM), based on weighted UniFrac distance, was performed with 999 iterations. p < 0.05 was considered statistically significant.

The general characteristics of all rats, including body weight, female hormones, and related indicators of glucose metabolism, are summarized in Table 1 and Figure 1 . Compared with those in the control group, the body weight, serum T, and Glu significantly increased in the model group (p = 0.000, 0.004, and 0.04); FSH, LH, HOMA-IR, and LPS had an increasing trend (p > 0.05); HOMA-beta had a downward trend (p > 0.05), showing the change of endocrine and metabolic disorders. Compared with those in the model group, the body weight in the treatment group significantly decreased (p = 0.043), and the FSH, HOMA-IR, and TNF-α had a downward trend (p > 0.05). There was no significant difference in E2, LH/FSH, and FINS among these groups (p > 0.05).

Under the light microscope, follicles in different developmental stages, including a mature follicle, and many corpora lutea were observed in the control group. Meanwhile, in the mature follicle, oocytes could be seen, and granulosa cells were arranged neatly and had multiple layers ( Figures 2A, B ). In the model group, many cystic dilated follicles with irregular structures can be seen, and the granulosa cell layer was reduced and even disappeared. Mature follicle and corpora lutea can hardly be seen ( Figures 2C, D ). The morphology of ovarian tissue in the treatment group was similar to that in the control group. A mature follicle and corpora lutea could be observed in the treatment group, and the layers of granulosa cells in follicles were more than those in the model group ( Figures 2E, F ).

α-Diversity index analysis can reflect the richness and diversity of intestinal microbiota communities in a specific regional environment; the greater the value, the higher the community diversity. The results are shown in the box diagram ( Figures 3A, B ). Compared with those in the control group, the Ace and Shannon indexes in the model group had an increasing trend (p > 0.05); meanwhile, the Ace index and Shannon index in the treatment group decreased significantly (p < 0.01, p < 0.05), which indicated that modified Banxia Xiexin Decoction can regulate the α-diversity of PCOS-IR rats.

The β-diversity of the intestinal microbiota was expressed by PCoA and PLS-DA. The PCoA results showed that the intestinal microbiota was significantly different among the three groups (R = 0.4193, p = 0.001, Figure 3C ). From PLS-DA, we can clearly see the distance of samples between groups. The greater the distance, the bigger the difference. Distances between groups meant that the samples in the model group can be clearly distinguished from other groups, while the samples in the treatment group can also be completely discriminated from the control group ( Figure 3D ).

On the phylum level, 5 phyla were identified, arranged in descending order of relative abundance, namely, Firmicutes, Verrucomicrobiota, Bacteroidota, Actinobacteriota, and Proteobacteria ( Figure 4A ). Compared with that in the control group, Verrucomicrobia in the model group showed a significant decrease (p = 0.005, Figure 4C ). Compared with those in the model group, Firmicutes and Actinobacteriota showed a significant decrease (p = 0.014 and 0.022), and the abundance of Verrucomicrobiota and Proteobacteria was significantly increased in the treatment group (p = 0.008 and 0.013, Figure 4D ).

On the genus level, identified bacterial taxa are represented in Figure 4B . The abundance of genera Akkermansia, Allobaculum, and norank_f:norank_o:Clostridia_UCG-014 in the model group was significantly less than that in the control group (p = 0.005, 0.013, and 0.008, Figure 4E ). Compared with those in the model group, genera Akkermansia and Escherichia-Shigella were significantly increased in the treatment group (p = 0.008 and 0.017, Figure 4F ), while Clostridium_sensu_stricto_1 showed an obvious decrease (p = 0.008, Figure 4F ).

Next, we used LEfSe to evaluate significant differences in the relative abundance of bacterial orders, classes, and genera among the three groups. In the control group, order Clostridia_UCG-014, order Staphylococcales, order Oscillospirales, and genus Allobaculum presented significant enrichment. The PCOS-IR model group was enriched in genus Clostridium_sensu_stricto_1 (order Clostridiales) and order Coriobacteriales (class Coriobacteriia). In the treatment group, the abundance of genera Akkermansia (phylum Verrucomicrobiota), Escherichia-Shigella (phylum Proteobacteria), and Blautia was more plentiful (LDA score >4.0, p < 0.05, Figure 4G ).

In the fecal samples of three groups, the COG function of gut microbiota was predicted including translation, ribosomal structure, and biogenesis; amino acid transport and metabolism; carbohydrate transport and metabolism; transcription; replication, recombination, and repair; cell wall/membrane/envelope biogenesis; inorganic ion transport and metabolism; and energy production and conversion ( Figure 5A ). At KEGG level 1, metabolism accounted for a large proportion of pathway prediction, followed by environmental information processing, genetic information processing, cellular processes, etc. ( Figure 5B ). At KEGG level 3, the predicted results indicated that the microbiota may be mainly involved in functions such as metabolic pathways, biosynthesis of secondary metabolites, microbial metabolism in diverse environments, biosynthesis of amino acids, and carbon metabolism to exert a therapeutic effect ( Figure 5C ).