Reagents of the LC-MS purity grade were acquired from Fisher Scientific (Fair Lawn, NJ, United Staes), including formic acid (purity >98%), acetonitrile, and methanol. Estradiol benzoate and progesterone injection were purchased from Zhejiang Xianju Pharmaceutical Co., Ltd. (Taizhou, China). ELISA kits of estradiol (E₂), testosterone (T), and progesterone (P) were obtained from Boster Biological Technology Co., Ltd. (Shanghai, China). McCoy’s 5A medium was provided by Gibco^® (Grand Island, NY, United States). Follicle-stimulating hormone (FSH) and luteinizing hormone (LH) were acquired from Shanghai Linc-Bio Science Co. Ltd. (Shanghai, China). Dexamethasone (DEX) and formononetin (FT) were supplied by Chengdu Must Biotechnology Co., Ltd. (Chengdu, China). The purity of FT was above 98%. Cell Counting Kit-8 (CCK-8) was purchased from Shanghai Beyotime Biotechnology Co., Ltd. (Shanghai, China). HHXYT (batch number: 190106), the same batch named as “PZ-5” in our previous study (Mi et al., 2020), was supplied by Puzheng Pharmaceutical Co., Ltd. (Jiangxi, China). The contents of fourteen quality marker compounds of HHXYT were quantified by HPLC-QQQ-MS (Table 1). The chromatograms and chemical structure were shown in Supplementary Figure S1.

All sexually mature and non-pregnant female Sprague-Dawley rats were obtained from Shanghai Xipuer-Bikai Experimental Animal Co. Ltd., Shanghai, China (License No.: SCXK-Shanghai-2018-0006). They were bred in a stabilized environment for 1 week with room temperature 25 ± 2°C, relative humidity 55 ± 5%, 12 h each of light and dark, and free access to drinking water and feed. All procedures were approved by the Animal Care and Use Committee of Second Military Medical University and abided by the Guide for the Care and Use of Laboratory Animals.

The raw data files were imported into the Progenesis QI 2.0 software (Waters, Milford, MA, United States) for data preprocessing, including alignment, peak picking, and compound identification. StatTarget was utilized to remove system errors as well as signal drift (Luan et al., 2018). The data was then imported into SICMA 14.1 (Umetric, Umeå, Sweden) for multivariable statistical analysis. Variates with p-values < 0.05, fold change (FC) ≥1.2 or FC ≤ 0.8, and VIP >1 (variable importance in the projection), were considered as endogenous differential metabolites, and further identified with online databases of HMDB (Wishart et al., 2018), LIPID MAPS (Liebisch et al., 2020) and KEGG (Kanehisa et al., 2017).

We applied the diffusion kernel algorithm to detect metabolites that locate at the proximity of the endogenous differential metabolites in the genome-scale metabolic network of rats (Luan et al., 2018). The diffusion parameter was set as 0.1. All endogenous differential metabolites of the MGH, HHXYT-L, and HHXYT-H groups were input as seed nodes to perform the diffusion kernel algorithm. The nodes with high diffusion scores were taken as compounds close to the endogenous differential metabolites in the genome-scale metabolic network. Then, the metabolic network affected by mammary glands hyperplasia and HHXYT treatment was constructed by extracting these nodes and links between them from the background network.

For prototypes and their metabolites in the plasma of rats after intragastric administration of HHXYT, we collected or predicted their putative targets from 3 databases and web servers: ETCM (Encyclopaedia of Traditional Chinese Medicine) (Xu et al., 2019), HERB (High-throughput Experiment- and Reference-guided dataBase of traditional Chinese medicine) (Fang et al., 2021), and TargetNet (Yao et al., 2016). Target information in HERB was collected by manual reference mining and from multiple TCM databases, including SymMap (Wu, Zhang, Yang et al., 2019), TCMID (Xue et al., 2013), TCMSP (Ru et al., 2014), and TCM-ID (Chen et al., 2006). The threshold of confidence score was set as 0.9 to filter the compound-protein interactions predicted by TargetNet.

Disease genes of MGH were collected from the GeneCards database (Stelzer et al., 2016). Inputting the keyword “mammary gland hyperplasia” generated 4,251 results with relevance scores in the area [0.53, 248.35]. We set the threshold of relevance score as 20 to filter the results and obtained 420 top-ranked genes with high relevance to MGH, taking about 10% of the results. These 420 genes were considered as MGH associated genes for further analysis.

We used the STRING platform to conduct functional enrichment analysis for gene groups under study (Szklarczyk et al., 2019).

The genomic scale protein-protein interaction (PPI) network of human beings was downloaded from the HINT (High-quality INTeractomes) database (Das and Yu, 2012), whose PPI data was integrated from 8 databases (BioGRID, MINT, iRefWeb, DIP, IntAct, HPRD, MIPS, and the PDB) and filtered both systematically and manually to keep high-confidence and -quality interactions. It was used as a background network for the construction of HHXYT’s target network. We used Cytoscape version 3.6.0 to construct, analysis and visualize networks (Shannon et al., 2003). The Cytohubba app of Cytoscape was applied to compute topological measures of nodes (Chin et al., 2014).

Metabolite-gene network was constructed by the Network Analysis module of MetaboAnalyst platform version 5.0 (Chong et al., 2019).

The SDF file for a ligand’s structure was downloaded from the PubChem database (Kim et al., 2021). The PDB file for a protein receptor’s crystal structure was obtained from RCSB Protein Data Bank (PDB) database (Burley et al., 2021). For proteins without known crystal structures, the predicted crystal structures were downloaded from AlphaFold Protein Structure Database (Jumper et al., 2021). The preprocessing and preparing of ligands and receptors were conducted by AutoDock Tools (Tanchuk et al., 2016) following the tutorial and manual (http://vina.scripps.edu/manual.html). Each ligand was docked with corresponding receptors by AutoDock Vina (Trott and Olson, 2010), and a compound-target pair with the docking score of less than −5 kcal/mol was considered as a binding pair. PLIP platform was used to analyze the binding sites of the binding compound-target pair (Adasme et al., 2021). PyMol was applied for visualizing the results of AutoDock Vina and PLIP (Alexander et al., 2011).

Thecal cells were extracted and purified from the ovaries of female rats at 3–4 weeks old (Li and Hearn, 2000). 5 × 10³ thecal cells per well were seeded into 96-well plates and set as the control group and experimental groups. The control group was McCoy’s 5A medium containing thecal cell. The experimental groups were based on the control group supplemented with different concentrations of luteinizing hormone (LH), follicle-stimulating hormone (FSH), HHXYT, and its ingredient formononetin (FT). After 36 h of adherent culture at 37°C and 5% CO₂, these drugs were supplemented according to different groups and cultured for another 24 h. Then 10 µL of CCK-8 reagent was added to each well in the dark, and the absorbance was measured at a wavelength of 450 nm 2 h later to calculate the cell survival rate (Chen et al., 2021).

4 × 10⁴ Cells were added to each well of 48-well plates containing McCoy’s 5A medium, which consisted of 10% FBS and 1% penicillin-streptomycin. After culturing at 37°C and 5% CO₂ for 36 h, different concentrations of LH, dexamethasone (DEX), HHXYT, and FT were added into different groups of wells and cultured for another 24 h (Ong et al., 2019). ELISA kits were used to detect the levels of estradiol (E₂) or testosterone (T) in cell media.

To evaluate the therapeutic effect of HHXYT on MGH, we detected variations of nipple diameter and serum sex hormone level. After establishing the MGH model, the nipple diameter grew significantly (p < 0.01, MGH vs. Sham) from 0.858 ± 0.068 mm to 1.230 ± 0.086 mm, compared to the Sham group. Different doses of HHXYT, on the other hand, could significantly (p < 0.01, HHXYT vs. MGH) reduce nipple enlargement (Figure 2A). Meanwhile, in the MGH group, the estradiol level in serum increased prominently (p < 0.01, MGH vs. Sham) while the progesterone level decreased notably (p < 0.01, MGH vs. Sham). After treatment with different doses of HHXYT, the serum level of estradiol and progesterone was reversed, respectively (Figures 2B,C).

According to Hematoxylin-Eosin (HE) staining, mammary tissue in the Sham group had no aberrant proliferation, a tiny number of the acinus, nearly no dilatation, and no secretion of mammary ducts and lumens (Figure 2D). The mammary tissue in the MGH group showed distinct proliferation; the volume of mammary lobules was increased; the number of glandular bubbles was also increased and became close to fusion; the acinus and lumens were revealed dilatation and excessive secretions (Figure 2E). However, after HHXYT treatment, mammary gland hyperplasia was prominently alleviated due to a considerable reduction in the account of the volume of mammary lobules, the number of glandular vesicles, and the secretion of lumens (Figure 2F).

It is widely accepted that prototypes and their metabolites in blood after oral TCM intake have a therapeutic effect (Ye et al., 2019). In order to determine the active components of HHXYT, we comprehensively analyzed the prototypes and their metabolites of HHXYT in rat plasma. A total of 54 components were identified, comprising 22 prototypes (Figures 3A,B; Supplementary Table S1) and 32 metabolites (Supplementary Table S2). After HHXYT treatment, three prototype components, named desbenzoylpaeoniflorin, 1-O-β-d–glucopyranosylpaeonisuffrone, and bohenoside A, as well as 32 metabolites, were identified in rat plasma for the first time. These 32 metabolites include 2 monoterpenoids, 8 triterpenes, 11 organic acids, 10 flavonoids, and 1 other metabolite (Figure 3C; Supplementary Table S4-S7).

Based on 54 active ingredients identified in the rat plasma, we predicted their potential target through network pharmacology (hydroxysafflower yellow A was also added to the analysis because it was reported to be an active ingredient of Carthami Flos (Honghua) (Li et al., 2021)). Using 3 databases and web servers (ETCM, HERB, and TargetNet), 373 distinct protein targets were predicted as HHXYT’s putative targets regulated by the 55 active ingredients (Supplementary Table S3).

The STRING platform was utilized to analyze the 373 targets for functional enrichment. Setting FDR<10⁻⁴ as the threshold, we identified 93 statistically enriched KEGG pathways. The top enriched pathways were presented in Table 2, and these pathways were mainly involved in the biological processes of the endocrine, immune, and nervous system.

Then, based on the KEGG pathways in Table 2, we performed enrichment analysis for targets of HHXYT’s herbs and prototypes. When conducting enrichment analysis, we combined targets of metabolites derived from its prototype with targets of the prototype. The threshold for statistical significance was set as FDR<10⁻⁴. As shown in Figure 4, targets of 6 herbs and 9 prototypes from HHXYT were enriched in at least 3 of the 20 KEGG pathways in Table 2, respectively. These results suggested that the 6 herbs and 9 prototypes could play key roles in HHXYT’s treatment of MGH.

On the basis of research on prototypes and their metabolites of HHXYT, the metabolomics analysis was used to further identify the serum endogenous differential metabolites of rats intervened by HHXYT. The PLS-DA model was applied to characterize the metabolic disturbances. As shown in Figures 5A,B, the serum samples of Sham, MGH, HHXYT-L, and HHXYT-H groups were distinctly separated in both positive and negative modes. Based on the criteria of p-values < 0.05, FC ≥ 1.2 or FC ≤ 0.8, and VIP >1, 58 endogenous differential metabolites in serum were screened and identified, including 29 metabolites in the Sham/MGH group, 19 metabolites in the MGH/HHXYT-L group, and 18 metabolites in the MGH/HHXYT-H group (Figure 5C). There were 5 same endogenous differential metabolites between HHXYT-H and MGH (Figure 5D, red compound names), as well as 2 overlaps between HHXYT-L and MGH (blue compound names). Furthermore, the regulatory trends of HHXYT-H on the 5 same differential metabolites were all opposite to that of MGH. Meanwhile, only one differential metabolite had been oppositely regulated between HHXYT-L and MGH groups. These results indicated that the HHXYT-H group had more pronounced alterations in metabolites endogenous on treating MGH. The majority of these endogenous differential metabolites are lipid metabolites, accounting for 68%, including 31% steroids lipid (12% bile acids, 8% pregnane steroids, 7% hydroxysteroids, etc.), 8% fatty acids, 7% fatty acid esters, 8% prenol lipids, 8% glycerophospholipids and so on (Figure 5E). The above results suggested that HHXYT could prominently regulate lipid metabolism, especially hormone disorder in rats with hyperplasia of mammary glands.

To reveal the underlying significance of metabolites variation in vivo, we proceeded with metabolic module establishment and pathway analysis of endogenous differential metabolites. Inputting all the endogenous differential metabolites as seed nodes, we ran the diffusion kernel algorithm to score all nodes in the genome-scale metabolic network of rats. The metabolic network influenced by MGH and HHXYT therapy was then constructed by extracting the nodes with diffusion scores at least 0.05 and the links between them (Figure 6). This network consists of 19 isolated modules, in which only the first three large modules include endogenous differential metabolites affected by both MGH and HHXYT treatment, suggesting that they could be key metabolic modules by which HHXYT exerts its therapeutic effect on MGH.

We searched the KEGG database (Kanehisa et al., 2017) to annotate the first 3 modules in Figure 6. It was found that metabolites in the first module were mainly involved in pathways of steroid hormone biosynthesis, cortisol synthesis and secretion, and ovarian steroidogenesis; the compounds in the second module participated in primary bile acid biosynthesis and bile secretion; the third module was associated with glycerophospholipid metabolism and biosynthesis of unsaturated fatty acids. Notably, all the 3 pathways associated with the first metabolic module and 1 pathway of the second metabolic module, were also found in the pathways listed in Table 2, which were enriched by 373 HHXYT’s predicted targets. It demonstrated that the pathways involved by HHXYT’s putative targets discovered by network pharmacology were consistent with the analyzed regulation of metabolic disturbances in rat serum by metabolomics.

The first step was to conduct a network analysis of HHXYT’s potential targets. We calculated the intersections of HHXYT’s 373 target genes with the 420 diseases genes for MGH (Supplementary Table S4; Figure 7A). It could be seen that one-fourth of HHXYT’s targets were disease genes. Then we constructed the protein-protein interaction (PPI) network between the targets based on the human genomic-scale PPI network from the HINT database. As shown in Figure 7B, 114 of the 373 targets interacted with each other to form HHXYT’s target network, with 53 of them being disease genes associated with MGH.

The centrality measures for nodes in the target network were computed using the Cytohubba App of the Cytoscape software. The centrality of nodes was then ranked applying 3 global parameters based on shortest paths (Closeness, EcCentricity, and Radiality). As shown in Figure 7C, the 4 nodes with the highest global centrality were AR, ESR1, AKT1, and JUN, suggesting that they were central nodes whose impacts were easy to reach out into the whole target network.

The Network Analysis module of the MetaboAnalyst platform version 5.0 was then utilized to construct metabolite-gene networks corresponding to the metabolic modules affected by HHXYT. We input the endogenous differential metabolites in the first and second modules of Figure 6, together with HHXYT’s targets enriched in the corresponding pathways of the modules, respectively. Then we got the metabolite-gene interaction modules for the first and second metabolite module affected by HHXYT (Figures 7D,E). As shown in Figure 7D, the first module included 3 hub genes, i.e., CYP11A1, HSD3B2, and CYP17A1, which linked to all the metabolites in the network, implying that they performed a key role in HHXYT’s treatment to MGH.

The key targets of HHXYT for treating MGH were then obtained by combining the 4 central genes of the target network and the 3 hub genes in the first metabolite-gene module. We integrated the 7 key targets with the key herbs, prototypes (Figure 4), and their metabolites on HHXYT’s treatment to MGH, to construct the herb-ingredient-target network as shown in Figure 8A. There were 7 targets, 6 herbs, and 17 ingredients in this network. 8 of 17 ingredients were HHXYT’s prototypes absorbed in plasma, and 9 were metabolites derived from them. These herbs, ingredients, and targets constitute key substances for HHXYT to exert therapeutic effects on MGH.

At last, to verify the binding between 5 of HHXYT’s key targets and corresponding prototypes, we utilized AutoDock Vina software to conduct molecular docking. The crystal structures of 4 targets, i.e., AR (PDB ID:2q7i), CYP11A1 (PDB ID:3n9y), CYP17A1 (PDB ID:6ciz), ESR1 (PDB ID: 2qgw), which have higher resolution, were downloaded from RCSB Protein Data Bank. The predicted structure of HSD3B2 (ID: AF-P26439-F1) was downloaded from AlphaFold Protein Structure Database. The SDF files for chemical structures of the compounds (ferulic acid, formononetin, glycyrrhetinic acid, glycyrrhizic acid, isoliquiritigenin, liquiritigenin, and protocatechuic acid) were downloaded from the PubChem database. The docking scores between 14 pairs of compound-target are less than −5 kcal/mol (Table 3), suggesting that the corresponding compounds and targets have a higher binding affinity. Binding poses and binding sites for 5 pairs of compound-target interactions were shown in Figures 8B–F.

Firstly, the cell-based assay was performed to demonstrate that the extracted and purified cells in vitro were thecal cells. After LH (50 ng/ml) was added, the cell proliferation rate (maximum ratio 147%) and testosterone (T) level (from 1965.15 ± 96.65 to 2,329.43 ± 87.93 pg/ml) both increased prominently (Figures 9A,B). However, after the addition of FSH, no significant changes in cell quantity were observed, and estradiol (E₂) was not detected. It indicated that the purified cells were primarily thecal cells, with no obvious interference from granulosa cells, which were also presented in follicles of the ovary. The cell survival rate results of HHXYT and its ingredient formononetin (FT, Figure 9C) were shown in Supplementary Figure S9.

To evaluate the effect of HHXYT and its ingredient formononetin (FT) on the production of testosterone in thecal cells, thecal cells were stimulated with different concentrations of dexamethasone (DEX) and then treated with different concentrations of HHXYT and FT. The results showed that testosterone levels were markedly elevated at the DEX concentration of 10, 20, 25 µM (Figure 9D). Whereas, HHXYT (50 µg/ml) and FT (3.2, 6.3, and 12.5 µM) could distinctly inhibit the increase of testosterone level caused by DEX (10 µM) in thecal cells, suggesting that HHXYT and FT have an inhibitory effect on testosterone production in thecal cells.