Adult female Sprague–Dawley rats, weighing 220-250 g, were obtained from the Laboratory Animal Center of Shanxi Cancer Institute (animal production certificate # SCXK (Jin) 2017-0001, Shanxi, China). The rats were housed in a controlled microbial environment at a constant temperature of 23 ± 2°C, with a 12-hour light/dark cycle, and had ad libitum access to sterile food and autoclaved water. After a week of acclimation, the rats were randomly assigned to the control group (CON, 7 rats), postmenopausal osteoporosis group (OVX, 7 rats), and Isaria felina treatment group (IF+OVX, 7 rats), with seven rats per group. The OVX and IF+OVX groups underwent bilateral ovariectomy to induce postmenopausal osteoporosis. Following the successful establishment of the model, the rats continued under the same conditions for 12 weeks. During this period, the IF+OVX group received 100 mg/kg of Isaria felina daily by oral gavage, while the control group received an equal volume of saline. The 100 mg/kg dose for Isaria felina was selected based on our previous study (27). After 12 weeks, the rats were euthanized by CO2 inhalation, and samples were collected for further analysis ( Figure 1 ). All animal experiments adhered to the National Institutes of Health (NIH) Guidelines for the Care and Use of Experimental Animals and were approved by the Ethical Committee of Experimental Animal Care of Shanxi Medical University (permit number: 2021014).

High-resolution micro-computed tomography (MicroCT) (vivaCT80, Scanco, Switzerland) scans of distal femurs were performed to analyze differences in trabecular bone volumes and structures between groups. One hundred contiguous cross-sectional slices above the femoral growth plate were manually selected for analysis. Histomorphometric parameters were computed using Scanco Medical software, determining bone mineral density (BMD), bone volume per tissue volume (BV/TV), trabecular number (Tb.N), trabecular separation (Tb.Sp), cortical bone thickness (Ct.Th), and cortical bone volume (Ct.V) for each sample.

Femoral samples were fixed in 4% paraformaldehyde (Solarbio, China) for 24 hours and decalcified in 20% ethylenediaminetetraacetic acid (EDTA) solution (Solarbio, China) at 37°C for 6–7 weeks until softened. The femurs were then dehydrated, embedded in paraffin, sectioned into 5-mm longitudinal slices, dried, and stained with hematoxylin and eosin (H&E) and tartrate-resistant acid phosphatase (TRAP) staining kits (Beyotime, China). Morphological characteristics were examined using light microscopy.

A three-point bending test was conducted to assess the mechanical properties of the femurs at the middiaphysis using an electronic universal testing machine (ElectroForce 3200 Series, TA Instruments, USA). The femoral samples were positioned and secured on a test stent with two fixed loading points spaced 20 mm apart. A constant displacement rate of 3 mm/min was applied until fractures occurred. Measurements of the internal and external major and minor axis lengths at the fracture points were taken. Peak load and maximum displacement values were recorded.

Blood samples were collected from the abdominal aorta of rats in each group. Three milliliters of peripheral blood were drawn into EDTA anticoagulant tubes (BD Vacutainer, USA) and mixed thoroughly. Subsequently, 30 µl of whole blood was transferred to a flow cytometry tube. For T-helper 17 (Th17) cell detection, CD4 (Mouse-Anti-Rat, BD, USA) and interleukin-17 (IL-17) (Thermo Scientific, USA) antibodies were sequentially added. For regulatory T (Treg) cell detection, CD4 (Mouse-Anti-Rat, BD, USA) and CD25 (Mouse-Anti-Rat, BD, USA) antibodies (4.5 µl each) were sequentially added, mixed thoroughly, and incubated for 15 minutes. Next, 200 µl of red blood cell (RBC) lysis buffer (BD Biosciences) was added, mixed thoroughly, and incubated for 10 minutes. The mixture was washed with saline, mixed thoroughly, and centrifuged at 1000 r/min for 5 minutes. The supernatant was discarded, and 400 µl of saline was added, mixed thoroughly, and then analyzed. Flow cytometry determined the proportions of IL-17+CD4+ T cells (Th17) and CD4+CD25+ T cells (Treg).

Blood samples were collected from the abdominal aorta of rats and allowed to stand at room temperature for approximately 30 minutes. Following this, the samples were centrifuged at 4°C and 3000 rpm for 20 minutes to separate the serum, which was then stored at -80°C. Interleukin-17 (IL-17), tumor necrosis factor-alpha (TNF-α), interleukin-4 (IL-4), and transforming growth factor-beta (TGF-β) levels were measured using enzyme-linked immunosorbent assay (ELISA) kits (IL-17, TNF-α, IL-4, and TGF-β kits from RayBiotech, USA).

Fecal samples were directly collected from the rats in both groups. At least two fecal pellets were obtained from each rat: one for microbial analysis and the other for metabolic analysis. The samples were placed in sterile centrifuge tubes, immediately frozen in liquid nitrogen, and stored at -80°C for further sequencing.

16S rDNA sequencing was performed at Lc-Bio Technologies Co., Ltd. Genomic DNA from fecal samples was extracted using the E.Z.N.A. Stool DNA Kit (Omega Bio-tek, Inc., USA) according to the manufacturer’s instructions. The V3-V4 region of the bacterial 16S rRNA gene was amplified using universal primers 338F (5’-ACTCCTACGGGAGGCAGCAG-3’) and 806R (5’-GGACTACHVGGGTWTCTAAT-3’). The amplified target bands were verified by 1% agarose gel electrophoresis at 170V for 30 minutes. PCR products were purified using the Agencourt AMPure XP nucleic acid purification kit (Beckman Coulter, Inc., USA). The library fragment size was verified using the 2100 Bioanalyzer (Agilent Technologies, Inc., USA). Finally, the library was sequenced on the Illumina Novaseq 6000 platform (Illumina, Inc., USA) with a PE250 sequencing strategy. The 16S rDNA sequencing methodology followed established protocols (17).

Sequencing data were split into different samples based on barcode sequences. The data were filtered and merged using Pear (v0.9.6) software. Sequences containing ambiguous bases (N) and those with a quality score below 20 were removed. During merging, the minimum overlap was set to 10 bp, with a p-value of 0.0001. Post-merging, sequences shorter than 230 bp were removed using Vsearch (v2.7.1) software. Chimeric sequences were identified and removed by comparison with the Gold Database using the uchime method. High-quality sequences were clustered into Operational Taxonomic Units (OTUs) using the uparse algorithm in Vsearch (v2.7.1) software, with a sequence similarity threshold of 97%. Alpha diversity and beta diversity analyses were performed using QIIME (v1.8.0) software based on the OTUs and their abundance. The OTU representative sequences were compared with the Silva138 database using the BLAST algorithm, with an e-value threshold of 1e-5, to obtain species classification information for each OTU. Based on the OTU and abundance data analysis, an in-depth exploration of the differences between comparison groups was conducted. To further optimize the data mining process and obtain more valuable information, various advanced analysis methods were employed. This analysis approach was informed by recent methodological advancements (17).

Initially, 25 mg of samples were mixed with 500 μL of extract solution (methanol: acetonitrile: water = 2:2:1) in an EP tube. The samples were homogenized at 35 Hz for 4 minutes and sonicated in an ice-water bath for 5 minutes, repeated three times. The mixture was centrifuged at 12,000 rpm for 15 minutes at 4°C, and the supernatant was transferred to a fresh glass vial for LC/MS analysis. A QC sample was prepared by combining equal aliquots of the supernatants from all samples. UHPLC-MS/MS analyses were conducted using a UHPLC (Vanquish, Thermo Fisher Scientific, USA) coupled to an Orbitrap Exploris 120 mass spectrometer (Thermo Fisher Scientific, USA) by Beijing Allwegene Technology Co., Ltd. Raw data were converted to mzXML format using ProteoWizard and processed with an in-house program based on XCMS, developed using R, for peak detection, extraction, alignment, and integration. The R package and AllwegeneDB were used for metabolite identification. The methodology used for UHPLC-MS/MS analysis aligns with previously validated techniques (28).

Metabolite features detected in more than 50% of experimental samples were excluded from data analysis. Missing values were imputed with half of the minimum value. Internal standard normalization was applied. Features with RSD greater than 30% were omitted from further analysis. The resulting three-dimensional data (peak number, sample name, and normalized peak area) were analyzed using MetaboAnalystR for orthogonal partial least squares discriminant analysis (OPLS-DA). OPLS-DA facilitated better group separation and identification of variables responsible for classification, with R2 and Q2 values indicating explained variation and prediction accuracy, respectively. To validate the robustness and predictive ability of the OPLS-DA model, 200 permutations were performed. The first principal component of variable importance in the projection (VIP) was obtained to refine the analysis. VIP values highlight each variable’s contribution to the model. Metabolites with VIP > 1, P < 0.05 (Student’s t-test), and fold change > 1.5 or < 0.67 were considered significantly altered. Metabolite pathways were identified using the KEGG (http://www.kegg.jp) and MetaboAnalyst (http://www.metaboanalyst.ca/) databases. This bioinformatics workflow follows best practices in metabolomics research (28).

All collected data were thoroughly processed using GraphPad Prism 9.0 statistical software. Student’s t-test was employed for statistical analysis between two groups. Differences in gut microbiota between two groups were identified using the Wilcoxon rank-sum test and the Kruskal-Wallis test. Differences in metabolites between groups were determined through univariate analysis, fold-change, and Student’s t-test. A P-value of less than 0.05 was considered statistically significant.

Micro-computed tomography (Micro-CT) analysis was conducted on the femurs of rats in each group. The results demonstrated that oral administration of Isaria felina mitigated bone loss in postmenopausal osteoporotic rats. Figure 2A presents the Micro-CT 3D reconstruction images of trabecular and cortical bone for each group. Various micro-CT parameters of trabecular and cortical bone, including bone mineral density (BMD), bone volume/total volume (BV/TV), trabecular number (Tb.N), trabecular separation (Tb.Sp), cortical thickness (Ct.Th), and cortical volume (Ct.V), are shown in Figures 2B–G . The findings indicated that oral Isaria felina significantly improved trabecular bone indices in ovariectomized (OVX) rats but did not significantly affect cortical bone.

Histopathological examinations of the femur, including hematoxylin and eosin (HE) staining and tartrate-resistant acid phosphatase (TRAP) staining, were performed to further understand the effects of oral Isaria felina on bone morphology in OVX rats. Compared to the control (CON) group, OVX rats exhibited trabecular bone fractures and thinning ( Figure 3A ), a significantly reduced trabecular area ratio, and decreased osteoblast numbers ( Figures 3B, C ). These conditions improved in the IF+OVX group. Additionally, the trabecular bone surface of OVX rats had numerous osteoclasts, indicating significant bone resorption, while the IF+OVX group showed a reduction in osteoclasts on the trabecular surface, resembling the CON group ( Figures 3D, E ). Oral Isaria felina effectively inhibited osteoclast activity and improved bone loss caused by postmenopausal osteoporosis.

Stress tests using a three-point bending test were conducted on the femurs of each group of rats. The maximum load-bearing capacity, fracture load, and stiffness of the femurs in the OVX group was significantly reduced compared to the CON group. In the IF+OVX group, the maximum load-bearing capacity, fracture load, and stiffness of the femurs improved, approaching that of the CON group. There was no significant difference in ultimate deflection between the groups ( Figures 4A–D ). The bone resorption and formation markers include CTX-I and PINP, respectively. The ELISA showed that the serum PINP was relatively higher in the OVX model because high bone turnover is associated with postmenopausal OP. The IF+OVX group showed significantly higher serum PINP levels than the OVX group, suggesting improved bone formation with Isaria felina treatment (Figure 4E). In contrast, serum CTX-I levels ( Figure 4F ) were significantly higher in the OVX group than in the CON group, indicating increased bone resorption. The IF+OVX group exhibited significantly lower serum CTX-I levels compared to the OVX group, indicating a reduction in bone resorption following Isaria felina treatment.

Throughout the oral administration of Isaria felina, the body weight of each rat group was consistently monitored. Compared to the CON group, the average body weight of OVX rats was notably higher, while the IF+OVX rats showed significantly lower body weight than the OVX rats ( Supplementary Figure S1A ). This suggests that Isaria felina effectively alleviates obesity resulting from postmenopausal osteoporosis, potentially contributing to mitigating bone loss. Following the treatment, uterine weight was assessed in each rat group. Both the OVX and IF+OVX groups exhibited markedly reduced uterine weight compared to the CON group ( Supplementary Figure S1B ), confirming the successful establishment of the OVX rat model.

A Venn diagram visually represented the shared and unique operational taxonomic units (OTUs) detected in each group through 16S rDNA technology, illustrating the OTU distribution across environmental samples and the intersections among samples or groups. The analysis revealed a total of 4423 OTUs across the three groups: 2729 OTUs were shared among all groups, 345 OTUs were exclusive to the CON group, 180 OTUs were unique to the OVX group, and 123 OTUs were distinctive to the IF+OVX group ( Figure 5A ).

Alpha diversity analysis, based on OTU clustering outcomes, was conducted. The rarefaction curve, generated by randomly selecting individuals from the samples and plotting them against species richness, served to evaluate diversity within samples at different sequencing depths and assess data collection adequacy. A leveling off of the curve indicated sufficient sequencing depth, where increasing data volume yields only a marginal increase in OTUs. Conversely, a continuing rise suggested potential for uncovering additional unknown OTUs with further sequencing, indicating the reasonableness of our sequencing depth ( Figure 5B ).

The species accumulation (Specaccum) curve, a graphical representation of species accumulation with increased sampling effort, was employed to evaluate species composition and estimate diversity within biological communities. This tool confirmed the adequacy of the sample size and estimated potential species richness in biodiversity and community structure investigations ( Figure 5C ).

In studies of microbial diversity, examining diversity within a single sample is pivotal for understanding the richness and diversity characteristics of microbial populations. Using Shannon and phylogenetic diversity (PD)_whole_tree indices, a significant reduction in gut microbial community diversity was observed in OVX rats compared to the CON group, while oral Isaria felina supplementation increased community diversity ( Figures 5D, E ).

Following OTU clustering, beta diversity analysis and species composition assessment were performed. Partial least squares discriminant analysis (PLS-DA) analysis, a supervised discriminant method, unveiled notable discrepancies among categories. The 2D and 3D PLS-DA analyses depicted distinct distributions of microbial communities across each group, highlighting significant differences in gut microbiota composition and abundance among the CON, OVX, and IF+OVX groups ( Figures 6A, B ). Subsequent species composition analysis spanned phylum, class, order, family, genus, and species levels. Genus-level species analysis findings, depicted through bar charts and heatmaps, revealed noteworthy alterations in gut microbiota abundance across various taxonomic levels among the groups ( Figures 6C, D ).

A comparative analysis of gut microbiota among the groups showcased significant differences at the genus level through bar charts and heatmaps ( Figures 7A–D ). Specifically, at the genus level, 20 distinct species differentiating the CON and OVX groups, along with 21 distinguishing species between the OVX and IF+OVX groups, were identified. Notably, oral Isaria felina supplementation reversed alterations in these species, including genera such as Bacteroides, Tuzzerella, Ruminococcus_torques_group, Barnesiella, Parabacteroides, Helicobacter, Alloprevotella, and Frisingicoccus.

Metabolites act as intermediaries through which gut microbiota influences the host. Untargeted metabolomics analysis was conducted on fecal samples from the three rat groups. Orthogonal partial least squares discriminant analysis (OPLS-DA) analysis revealed notable distinctions in fecal metabolite profiles between the CON and OVX groups, as well as between the OVX and IF+OVX groups ( Figures 8A, C ). In this investigation, 108 distinct metabolites were identified between the CON and OVX groups, with 69 downregulated and 39 upregulated in the OVX group ( Figure 8B ). Between the OVX and IF+OVX groups, 51 differential metabolites were identified, with 31 downregulated and 20 upregulated in the IF+OVX group ( Figure 8D ).

To explore the variations in fecal metabolites among the groups, KEGG pathway enrichment analysis was performed ( Figures 9A, B ). The results highlighted that the metabolic pathways influenced by oral Isaria felina primarily involved lipid and nucleotide metabolism. These pathways included alpha-linolenic acid metabolism, steroid biosynthesis, nucleotide metabolism, and purine metabolism.

Peripheral blood was collected from rats via the abdominal aorta for flow cytometry analysis following oral Isaria felina treatment. The findings revealed a significant increase in interleukin-17-positive (IL-17+)+CD4+ T cells (Th17) and a notable decrease in CD4+CD25+ T cells (Treg) in the OVX group compared to the CON group. Subsequent treatment with oral Isaria felina significantly reduced Th17 cells and elevated Treg cells in the IF+OVX group compared to the OVX group, although the increase in Treg cells was not statistically significant ( Figures 10A–D ). Furthermore, ELISA assessed cytokines associated with Th17 and Treg cell secretion, including IL-17 and TNF-α for Th17 cells, and IL-4 and TGF-β for Treg cells. The results indicated that, relative to the OVX group, levels of IL-17 and TNF-α significantly decreased, while IL-4 and TGF-β levels increased in the IF+OVX group post-treatment, although the increases were not statistically significant ( Figures 10E–H ). These findings suggest that oral Isaria felina has the potential to modulate and rectify the Th17/Treg immune imbalance in OVX rats.