Network pharmacology integrated molecular docking reveals the bioactive components and potential targets of Morinda officinalis–Lycium barbarum coupled-herbs against oligoasthenozoospermia

Oligoasthenozoospermia (OA) is one of the most common types of male infertility affecting sperm count and sperm motility. Unfortunately, it is difficult for existing drugs to fundamentally improve the sperm quality of OA patients, because the pathological mechanism of OA has not been fully elucidated yet. Morinda officinalis–Lycium barbarum coupled-herbs (MOLBCH), as traditional Chinese Medicines, has been widely used for treating OA over thousands of years, but its molecular mechanism is still unclear. For this purpose, we adopted a comprehensive approach integrated network pharmacology and molecular docking to reveal the bioactive components and potential targets of MOLBCH against OA. The results showed that MOLBCH alleviated apoptosis, promoted male reproductive function, and reduced oxidant stress in the treatment of OA. Ohioensin-A, quercetin, beta-sitosterol and sitosterol were the key bioactive components. Androgen receptor (AR), Estrogen receptor (ESR1), Mitogen-activated protein kinase 3 (MAPK3), RAC-alpha serine/threonine-protein kinase (AKT1), Glyceraldehyde-3-phosphate dehydrogenase (GAPDH) were the core potential targets. PI3K/Akt signaling pathway, prostate cancer, AGE-RAGE signaling pathway in diabetic complications were the most representative pathways. Moreover, molecular docking was performed to validate the strong binding interactions between the obtained core components and targets. These observations provide deeper insight into the pathogenesis of OA and can be used to design new drugs and develop new therapeutic instructions to treat OA.


Results
MOLBCH component-target network. 354 components of MOLBCH were obtained from TCMSP and TCMID. Among them, 174 components were from MO, 202 components were from LB, and 22 common components were from MO and LB. Subsequently, OB ≥ 30% and DL ≥ 0. 18 were used as the screening criteria. Finally, we got 66 bioactive components of MOLBCH, including 20 from MO, 48 from LB, and 2 (beta-sitosterol, sitosterol) from MO and LB (Fig. 1a,b). Then, the structural information of 66 bioactive components was collected from PubChem and ALOGPS2.1 (Supplementary Table S2). Four public webservers, Swiss Target Prediction, SEA, TCMSP and Drugbank, were used to predict the potential targets of the bioactive components according to the similarity-based method. 671 potential targets were predicted from 65 obtained bioactive components after removing duplicates, except one non-effective component (cyanin). The MOLBCH component-target network was constructed by Cytoscape software, including 736 nodes and 3034 edges (Fig. 1c). Among the obtained nodes, 372 targets were from 20 components of MO, 504 targets were from 47 components of LB, and 205 common targets were from the common components of MO and LB. According to these results, we suggest that MO and LB act on an integrated effect on OAthrough the common targets to a certain extent.
MOLBCH-OA common-target network. The pathogenesis of OA is resulted from the co-occurring condition of oligospermia and asthenozoospermia. Therefore, taking the genes of these two diseases into consideration is essential to reveal the common targets of MOLBCH on OA. The number of targets got from oligospermia and asthenozoospermia is 993 and 683, while 473 overlapping targets were collected (Fig. 2a). 31 targets from the search term "oligoasthenozoospermia" were added to get total 495 OA-related targets (Fig. 2b). In order to ensure the comprehensiveness of target collection, we used five human genomic databases, namely DisGeNET, CTD, OMIM, GeneCards and NCBI. The number of targets from these databases was 51, 415, 36, 78 and 30, respectively.
Additionally, 136 common targets were obtained from the intersection of OA-related targets and MOLBCH targets (Fig. 2c). 19 bioactive components from MO and 43 bioactive components from LB were associated with 136 MOLBCH-OA common targets (Fig. 2d). Furthermore, the degree value of the MOLBCH-OA commontarget network was calculated by Network Analyzer, a plugin of Cytoscape, which was used as a screening condition to get the key bioactive components (Supplementary Table S6). We found four key bioactive components of the MOLBCH-OA common-target network, namely Ohioensin-A from MO, quercetin from LB, beta-sitosterol and sitosterol from both MO and LB, respectively. Meanwhile, the degree value of Androgen receptor (AR),    Table S9). The top 20 significant terms are shown in Fig. 4a-d. Specifically, BP is related to "response to oxidative stress (OS)", "reactive oxygen species metabolic process", "response to oxygen levels", "cellular response to oxidative stress", "response to hypoxia", "response to decreased oxygen levels", "regulation of reactive oxygen species metabolic process", "response to reactive oxygen species" and "positive regulation of reactive oxygen species metabolic process", indicating that MOLBCH has an anti-oxidant effect to regulate male reproduction on OA (Fig. 4a). CC is associated with "nuclear envelope", "membrane raft", "membrane microdomain", "membrane region", "organelle outer membrane", "outer membrane", "mitochondrial outer membrane", "transcription factor complex", "nuclear chromatin", "RNA polymerase II transcription factor complex", "protein kinase complex", "serine/threonine protein kinase complex" and "cyclin-dependent protein kinase holoenzyme complex", demonstrating that MOLBCH can regulate male reproductive function by acting on membrane and protein kinase complex in cellular (Fig. 4b). MF is relevant to "protein serine/threonine kinase activity", "protein tyrosine kinase activity", "transmembrane receptor protein tyrosine kinase activity" and "protein serine/threonine/tyrosine kinase activity", revealing that MOLBCH could affect protein kinase activity in the pathogenesis of OA (Fig. 4c).
The KEGG database was conducted to investigate the pathways related to the possible functions of the PPI network (Supplementary Table S9). The top 20 significant pathways were shown in Fig. 4d. The results indicate that MOLBCH regulates apoptosis process through "PI3K-Akt signaling pathway", "MAPK signaling pathway", "Apoptosis", "IL-17 signaling pathway" "TNF signaling pathway". In addition, "AGE-RAGE signaling pathway in diabetic complications" and "HIF-1 signaling pathway" are related to oxidant stress in the course of disease. Besides, MOLBCH could regulate the male reproductive function by affecting "Prostate cancer", "Endocrine resistance", "Relaxin signaling pathway", "EGFR tyrosine kinase inhibitor resistance" and "Prolactin signaling pathway". Therefore, through the GO and KEGG pathway enrichment analyses of MOLBCH-OA PPI network, we believe that MOLBCH might treat OA via promoting male productive function, reducing OS, and inhibiting apoptosis process.
GO and KEGG pathway enrichment analyses of the cluster. Network cluster is defined as a set of highly interconnected nodes, which is helpful to discover and reveal the hidden biological information in the network 19 . In order to identify the potential mechanism of the 136 common targets, the MOLBCH-OA PPI network was divided into 6 clusters (Fig. 5). According the result of Network Analyzer, the degree value of GAPDH, AKT1, CASP3, Interleukin-6 (IL6), VEGFA, Mitogen-activated protein kinase 3 (MAPK3), Myc proto-oncogene protein (MYC), ESR1, and Epidermal growth factor receptor (EGFR) is much higher than other proteins in cluster 1.
Since the score of cluster 1 was much higher than other clusters, so we performed GO and KEGG pathway enrichment analyses to further investigate the proteins of cluster 1 (Fig. 6). "Response to OS", "cellular response to OS", "reproductive structure development", "reproductive system development" in BP suggest that the proteins in cluster 1 are related to anti-oxidant effect and male reproductive regulation (Fig. 6a). "Nuclear envelope", "nuclear membrane", "membrane raft", "membrane microdomain", "membrane region", "organelle outer membrane", "outer membrane", "mitochondrial outer membrane", "nuclear inner membrane", "protein kinase complex", "RNA polymerase II transcription factor complex" and "cyclin-dependent protein kinase holoenzyme complex" in CC demonstrate that the proteins in cluster 1 are relevant to membrane and protein kinase complex in cellular (Fig. 6b). "DNA-binding transcription factor binding", "RNA polymerase II-specific DNA-binding transcription factor binding", "activating transcription factor binding" and "core promoter sequence-specific DNA binding" in MF reveals that the proteins in cluster 1 could affect DNA binding (Fig. 6c).
The KEGG database was conducted to investigate the pathways related to the possible functions of the PPI network (Supplementary Table S10). The top 20 significant KEGG pathways were shown in Fig. 6d. The KEGG results indicate that the proteins in cluster 1 are related to inhibit apoptosis process through "PI3K-Akt signaling pathway" and "IL-17 signaling pathway". The proteins in cluster 1 also have connections with promoting male reproductive function by affecting "Prostate cancer", "Endocrine resistance", "Relaxin signaling pathway" and "Prolactin signaling pathway". The GO and KEGG pathway enrichment analyses indicate that the proteins in cluster 1 are mainly related to oxidant stress, apoptosis, and male reproductive function, which consistents with the results of MOLBCH-OA PPI network analysis ( Fig. 6d- www.nature.com/scientificreports/  www.nature.com/scientificreports/  www.nature.com/scientificreports/ quercetin, respectively. The common ingredients beta-sitosterol and sitosterol from MO and LB is also of great importance during the course of drug treatment. As shown in Table 1, GO and KEGG analyses of MOLBCH-OA PPI network and clusters indicate that MOLBCH could treat OA by regulating male reproductive function, reducing apoptosis and OS. The main KEGG signaling pathways related to the above mechanism are PI3K-Akt signaling pathway, MAPK signaling pathway, Apoptosis, IL-17 signaling pathway, TNF signaling pathway, AGE-RAGE signaling pathway in diabetic complications, HIF-1 signaling pathway, Prostate cancer, Endocrine resistance, Relaxin signaling pathway, EGFR tyrosine kinase inhibitor resistance. The targets relevant to the above pathways were shown in Fig. 6, AR, ESR2, ESR1, Cytochrome P450 17A1 (CYP17A1), MAPK3, M-phase inducer phosphatase 2 (CDC25B), Nitric oxide synthase (NOS2), AKT1, Dual specificity mitogen-activated protein kinase kinase 1 (MAP2K1), and EGFR are the top 10 vital targets in core-component-target-pathway network. In addition, GAPDH, AKT1, CASP3, VEGFA, MYC, EGFR, IL6, MAPK3, and ESR1 also play an essential role in MOLBCH-OA PPI network and clusters. We found that the duplicated targets between them are ESR1, MAPK3, AKT1, and EGFR. Since AR is the most significant target in MOLBCH-OA common-target network and core component-target-pathway network, GAPDH is the most significant in the MOLBCH-OA PPI network and clusters, we suggest that AR, ESR1, MAPK3, AKT1, and GAPDH are the core potential targets of MOLBCH against OA. Based on the above results, we hold that PI3K-Akt signaling pathway, Prostate cancer, and AGE-RAGE signaling pathway in diabetic complications are the most important signaling pathways during the process of MOLBCH in treating OA. The targets related to these three signaling pathways are shown in Figs. 8, 9 and 10. Furthermore, we chose Ohioensin-A, quercetin, betasitosterol, sitosterol as the molecular docking ligands, AR, ESR1, MAPK3, AKT1, and GAPDH as the targets to reveal the interaction between the key bioactive components and core potential targets of MOLBCH against OA. Molecular docking. Molecular docking is a procedure using molecular modeling techniques to predict how a protein interacts with small molecules (ligands) 24 . In this study, we chose Ohioensin-A from MO, quercetin from LB, and the common ingredients beta-sitosterol and sitosterol from MO and LB as small molecules (ligands), AR, ESR1, MAPK3, AKT1, and GAPDH as proteins to perform the molecular docking. These proteins not only play an important role in the KEGG signaling pathways, but also serve as the key nodes of the PPI network and clusters. The obtained docking results indicate that the receptor-ligand interaction between drugs and proteins includes hydrophobic interactions and polar interactions. According to Tables 2 and 3, Ohioensin-A, quercetin, beta-sitosterol and sitosterol have strong binding interactions with AR, ESR1, MAPK3, AKT1, and GAPDH.

Discussion
With the changes in people's living habits, environmental pollution, and psychological factors, the incidence rate of infertility continues to rise, where male infertility accounts for 50% of the cases 25 . As one of the most common types of male infertility, OA has a complicated mechanism of action, thus the current treatment methods and drugs have some limitations 26 . So far, the treatment of OA has mainly focused on hormones, anti-infection, surgery, and ART 27 . However, the existing therapeutic drugs and surgical methods could not fundamentally improve the sperm quality of patients with OA. In addition, surgical treatment has brought economic and psychological pressure to infertile couples. Under the current influence of COVID-19 sweeping the globe, researchers have found that the COVID-19 virus can damage the male reproductive system [5][6][7][8] . Therefore, it is worth to be further explored to improve sperm quality of OA patients, and develop new drugs against OA.
Here, we adopted a comprehensive method integrated network pharmacology and molecular docking to reveal the bioactive components and potential targets of MOLBCH against OA. The experimental flow of this study was shown in Fig. 15. In our study, for the first time, Ohioensin-A, quercetin, beta-sitosterol and sitosterol were found to be the main bioactive ingredients of MOLBCH against OA. Specifically, the cytotoxic activity of Ohioensin-A has good effects against various cancer cell lines, including murine leukemia cell line and breast cancer cell line 20 . Ohioensin could reduce the TNF-α-induced production of intracellular reactive oxygen species (ROS) and phosphorylation of AKT in vascular smooth muscle cells (VSMCs) 28 . In previous study, quercetin was confirmed to indirectly affect the stimulation of the sex organs, both at the cellular and organ levels 21 , and showed outstanding beneficial effects on the serum total testosterone 22 . The supplement of quercetin could decrease the expressions of AKT, AR, cell proliferative and anti-apoptotic proteins on prostate cancer in the in vivo model 29 . The stimulation of cell proliferation by quercetin is proved to be mediated by ESR1 30 . In a previous study, researchers found quercetin could elicit apoptosis through an ESR1-dependent mechanism www.nature.com/scientificreports/ in cancer cell lines 31 . Quercetin has a protective effect against chronic prostatitis in rat model through NF-κB and MAPK signaling pathways 32 , and could attenuate cell migration and invasion by suppressing the protein levels of p-AKT1, MMP-2, and MMP-9 in HCCLM3 cells 33 . Beta-sitosterol and sitosterol are natural occurring phytosterols with steroidal moiety, which could inhibit tumor growth, modulates immune response, and has antioxidant capacity. Beta-sitosterol is regarded as a potential chemo preventive agent for treating a variety of cancer, including prostatic carcinoma and breast cancer 23 . It has been reported that beta-sitosterol could inhibit the growth and migration of prostate cancer cell and slow the growth of prostate tumors in mice. Its mechanism of action could be involved in AR 34 . Incorporation of beta-sitosterol into the cell membrane could increase the resistance to OS and lipid peroxidation via ESR1-mediated PI3K/GSK3β signaling 35 . Beta-sitosterol could increase the tyrosine phosphorylation of IRS-1, serine phosphorylation of AKT, threonine phosphorylation of AKT and threonine phosphorylation of Akt substrate of 160 KD in the adipose tissue of type-2 diabetic rat 36 .
In addition, we found that the core potential targets of MOLBCH on OA were AR, ESR1, MAPK3, AKT1 and GAPDH. AR is essential for the development and maintenance of the male phenotype and spermatogenesis 37,38 . Previous studies demonstrated AR might work through testicular Sertoli and peri-tubular myoid cells, maintaining spermatogonia numbers, blood-testis barrier integrity, completion of meiosis, adhesion of spermatids and spermiation 39,40 . ESR1 could regulate expression of genes during the process of spermiogenesis, and has been implicated in male infertility [41][42][43] . Besides, ESR1 was found to be associated with testicular germ cell cancer, which usually occurs in young men 44 . The MAPKs has been linked to the disturbances in spermatogenesis and dysfunction of germ cells and Sertoli cells, resulting in reduced semen quality and male reproductive dysfunction 45 . In human, MAPK3 may play a crucial role in cell cycle progression and apoptosis 46 . AKT1 is considered as the moderator of cellular growth, survival, metabolism and proliferation 47 . AKT1 could also suppress radiation-induced germ cell apoptosis in vivo 48 and enhance the effects of thyroid hormone on postnatal testis development 49 . In the testis, GAPDH is of particular importance for spermatogenesis, and could reduce sperm motility induced by male infertility 50 . Besides, EGFR, ESR2, MYC, CASP3, VEGFA, etc. are also important during the process of MOLBCH against OA. EGFR is located in the head and middle of the sperm, and participated in the acrosome reaction and the polymerization reaction of actin in the process of sperm capacitation 51,52 . ESR2 regulates the expression of genes related to germ cell apoptosis and spermatization 41 . c-MYC is an immediate growth early response gene, and might play a key role in cell proliferation and tumorigenesis 53,54 . The increasing expression of CASP3 in the testis may cause germ cell apoptosis in the seminiferous tubules 55 . VEGFA is related to male reproductive function and maintenance of spermatogonial stem cells 56 . www.nature.com/scientificreports/  www.nature.com/scientificreports/ Furthermore, 12 significant pathways were related to apoptosis, male reproductive functions, and OS in the molecular mechanism of MOLBCH against OA. In particular, the most represented signaling pathways are PI3K-Akt signaling pathway, Prostate cancer, and AGE-RAGE signaling pathway in diabetic complications. The PI3K/Akt signaling pathway plays an essential role in inhibiting cell apoptosis and promoting the survival of male infertility 57 . It also perturbs the intracellular redox equilibration by raising ROS 58 . The aberrant www.nature.com/scientificreports/ activation of PI3K-Akt signaling pathway may contribute to increase cell invasiveness and facilitate prostate cancer progression 59 . Prostate cancer is the second leading cause of cancer death in men, and key determinants of its cellular phenotype include carcinogen defense (GSTP1), growth factor signaling pathways (NKX3.1, PTEN and p27) and AR 60 . Activation of AGE/RAGE signaling pathway promotes ROS production, which can induce lipid peroxidation and eventually sperm DNA damage [61][62][63][64] . www.nature.com/scientificreports/ Besides, MAPK signaling pathway, IL-17 signaling pathway, TNF signaling pathway, Endocrine resistance, Relaxin signaling pathway, EGFR tyrosine kinase inhibitor resistance, Prolactin signaling pathway, HIF-1 signaling pathway, also play an important role in the process of MOLBCH on OA. The MAPK signaling pathway participates in many stages of germ cell development, including spermatogenesis, germ cell cycle progression, germ cell apoptosis, acquisition of motility in the epididymis, sperm capacitation and acrosome reaction before www.nature.com/scientificreports/ the fertilization of oocytes 65,66 . The aberrant IL-17 signaling pathway is of great importance to maintain the testicular immune, including cell immunity, mucosal immunity and cytokines [67][68][69] . TNF family is regarded to stimulate NF-κB, and further acts an effect on varicocele-mediated pathogenesis 70 . Estrogens may be involved in the pathophysiology of varicocele-associated male infertility 71 . In male mice, disruption of the relaxin gene results in the delayed reproductive tract development and arrested sperm maturation 72 . EGFR is a tyrosine kinase www.nature.com/scientificreports/ related to the regulation of cellular homeostasis. Mutations in the EGFR gene and protein overexpression could activate downstream pathways associated with cancer 73 . Prolactin receptor deficient models built in a previous study showed neuroendocrine and reproductive abnormalities for male rodents 74 . Hypoxia-Inducible Factor (HIF)-1 plays an integral role in responding to low oxygen concentrations or hypoxia in human 75 . Therefore, we believe that MOLBCH could improve male reproductive functions, decrease apoptosis and OS in the treatment of OA, which might further clarify the pathological mechanism of OA. In summary, we conducted a comprehensive approach integrated network pharmacology and molecular docking to demonstrate the key bioactive components and core potential targets of MOLBCH against OA. We found that MOLCH could alleviate apoptosis, promote male reproductive function, and reduce OS in the treatment of OA. The key bioactive components are Ohioensin-A, quercetin, beta-sitosterol and sitosterol. The core potential targets are AR, ESR1, MAPK3, AKT1 and GAPDH. The most representative pathways are PI3K/Akt signaling pathway, prostate cancer, and AGE-RAGE signaling pathway in diabetic complications. In order to further verify the results of network pharmacology, molecular docking was employed to confirm the strong binding interaction www.nature.com/scientificreports/ between the key bioactive components and core potential targets. This study provides deeper insights into the pathogenesis of OA and can be helpful to design new drugs and develop new therapeutic instructions to treat OA.

Methods
Data collection and processing. Components 77 were used to obtain the components of MOLBCH (Supplementary Table S1).
Bioactive components of MOLBCH. The adsorption, distribution, metabolism, and excretion (ADME)-related models, integrating oral bioavailability (OB) and drug-likeness (DL), were used to filter the bioactive components of MOLBCH. Oral bioavailability could represent the relative amount of orally administered drug that reaches the blood circulation 78 . Drug-likeliness was used to describe and optimize pharmacokinetic and pharmaceutical properties 79 . According to the published literature 16,19,80,81 , the compounds that meet the requirements of both OB ≥ 30% and DL ≥ 0.18, were used to identify the bioactive components of MOLBCH (Supplementary Table S2).
MOLBCH targets.  90 . Taking advantage of the different characteristics of each database, we selected "oligoasthenozoospermia", "oligospermia", "oligozoospermia", "asthenospermia" and "asthenozoospermia" as the keywords and criteria to search related targets. We integrated duplicated targets of "oligospermia" and "oligozoospermia", and defined the obtained targets as the oligospermia-related targets. Similarly, the same targets of "asthenospermia" and "asthenozoospermia" were incorporated, and defined the obtained targets as the asthenozoospermia-related targets. Afterwards, we took the intersection of the oligospermia-related targets and the asthenozoospermia-related targets, and the targets searched by "oligoasthenozoospermia" were added to get the OA-related targets. Finally, we standardized the names of OA-related targets using UniProtKB (Supplementary  Table S4).
Network construction. MOLBCH-OA common-target network. The MOLBCH-OA common target network was constructed by the common targets of MOLBCH and OA Cytoscape (http://www.cytos cape.org, version 3.7.2) 91 (Supplementary Table S5). Network Analyzer 92 , a plugin of Cytoscape, was used for calculating the degree value of the network, aiming to find the underlying components and targets of MOLBCH on OA (Supplementary Table S6).
PPI network and evaluation. The STRING database v11.0 (http://strin g-db.org) 93 was utilized to get the protein-protein interaction (PPI) information. The confidence score was set to 0.4 or higher. The PPI information were exported in TSV format, and then visualized by Cytoscape. CytoNCA, a plugin of Cytoscape, was used to evaluate the PPI network, including Betweenness (BC), Closeness (CC), Eigenvector (EC), Local Average Connectivity-based method (LAC), Network (NC), Subgraph (SC), Information (IC) 94 . First, the data was screened by using the screening criteria of 'DC ≥ 2 × median DC' . Then, the potential targets were obtained by the screening criteria of 'DC, BC, EC, CC, LAC, and NC greater than or equal to their median' 95 (Supplementary Table S7).
Cluster analysis. MCODE 96 , a cluster analysis algorithm in Cytoscape, was performed to analysis the sub-region of PPI network. The same or similar targets were treated as clusters to further explore the underly information the PPI network 97 (Supplementary Table S8). The conditions were set as node score cutoff = 0.2, K-core = 2, and degree of cutoff = 2.

Data availability
The data generated or analyzed during this study are included in the Supplementary Source files.