Eighteen female Sprague-Dawley (SD) rats of 5-week-old (Charles River Laboratory Animal Technology, Beijing, China) were housed in Animal Centre of Zhejiang Chinese Medical University (Hangzhou, Zhejiang, China), accredited by Association for Assessment and Accreditation of Laboratory Animal Care (AAALAC). The experimental animals were maintained in a specific-pathogen-free (SPF) grade conditions with temperatures ranging from 22°C to 24°C and humidity ranging from 40°C to 65°C. The rats had unrestricted access to regular feed and drinking water in their enclosures at the Laboratory Animal Research Center. The Animal Ethics Committee of Zhejiang Chinese Medical University authorized all experimental procedures (Approval No. 20220815-11).

After adaptive feeding for 1 week, the rats were equally and randomly assigned three groups (by Ying Zhu): control, darkness, and melatonin group (n=6). Figure 1A illustrates the experimental design of the animal study. The control group was kept under a normal 12-hour light/dark (L/D) alternating exposure (illumination beginning at 7:30 AM, Zeitgeber Time 0, ZT0). Both the darkness and melatonin group were housed in constant darkness for 8 weeks. In the melatonin group, rats were intraperitoneally administered melatonin daily at ZT11, with a dosage of 10 mg/kg/day, which was widely preferred in rodent models (10, 28). The control and darkness group were intraperitoneally administered with saline of equal volume. Whole rat body weight was measured weekly throughout the 8-week period. Rats were sacrificed under the effects of 5% isoflurane anesthesia during the diestrus II stage for tissue sample collection.

For eight consecutive days during rat model development, vaginal smears were examined to determine the estrous cycle stages (by Ran Cheng and Ying Zhu). The procedure involved collecting vaginal secretions from each rat between 9:00 and 10:00 a.m. by rinsing the vagina with a saline solution. Collected secretions were then placed on a glass slide for microscopic analysis. The estrous cycle stages were classified into four categories: metestrus (diestrus I), characterized by a small number of cornified cells and numerous leukocytes; diestrus (diestrus II), primarily composed of leukocytes; proestrus, distinguished by a high concentration of nucleated epithelial cells and estrus, dominated by large, flat, irregularly shaped a nucleated cornified epithelial cells.

Peripheral blood mononuclear cells (PBMCs) were obtained from whole blood by standard Ficoll-Hypaque density centrifugation using rat peripheral blood lymphocyte isolation fluid (#KGA3103, KeyGEN Biotechnology) in accordance with the producer’s instructions (by Ran Cheng and Ying Zhu). The mononuclear cells underwent three washing cycles using PBS before being preserved at -80°C for subsequent extraction of total RNA.

Rat ovarian tissues were harvested and fixed with 4% paraformaldehyde prior to paraffin embedding (by Sheng-kai Wang). Sections of tissue, 5 μm in thickness, were prepared and subsequently subjected to dewaxing and rehydration using a graded series of ethanol solutions. Tissue samples were colored with hematoxylin for 5 min, and differentiated with hydrochloric acid for 30 s. Finally, the specimens were stained with eosin for 2 minutes and mounted on slides for microscopic examination.

The concentrations of LH and testosterone in rat sera were detected using the Rat LH ELISA Kit (#ELK2367, ELK Biotechnology) and Testosterone ELISA Kit (#ELK1332, ELK Biotechnology), respectively (by Chun-xiao Zong). The protocols were performed in accordance with the guidelines.

Total RNA from PBMCs and rat ovary tissues was isolated with TRIzol reagent (Invitrogen) and then reverse-transcribed into cDNA (#RR036A, Takara Bio, Dalian, China). RT-qPCR (#RR820A, Takara Bio, Dalian, China) was employed to measure the mRNA expression of target genes (by Ying Zhu). The ΔΔCt method was subsequently utilized for data analysis. The target mRNA level was calculated as the ratio of the target gene to GAPDH. Table 1 tabulates the primer sequences used for the target genes.

To obtain the spatial configuration and Simplified Molecular-Input Line-Entry Specification (SMILES) notation of melatonin, a search was conducted in the PubChem database (https://pubchem.ncbi.nlm.nih.gov/) (29). Subsequently, the 3D structure was input into the PharmMapper-server for potential target prediction, whereas SMILES was submitted to SwissTargetPrediction (http://www.swisstargetprediction.ch/) (30, 31). Furthermore, the DrugBank and Traditional Chinese Medicine Systems Pharmacology (TCMSP) (https://tcmspw.com/tcmsp.php) were utilized to recognized melatonin targets (32, 33). The retrieved target names were verified against authorized names in the UniProt database (https://www.uniprot.org/), and duplicate entries were eliminated. Cytoscape 3.8.2 was employed to visualize the results (by Ran Cheng and Chun-xiao Zong).

To identify targets associated with PCOS, researchers (by Ran Cheng and Sheng-kai Wang) have utilized three gene databases: GeneCards, DisGeNET, and NCBI (34, 35). The search terms “Polycystic ovary syndrome” and “PCOS” were employed. To ensure result reliability, GeneCards had a minimum relevance score of minimum 0.7, and a gene-disease score of 0.1 or higher, was applied. The STRING database (version 11.5) was used to predict the network of interactions between proteins (PPI) for the given gene sets with a high-confidence interaction score (0.900). Cytoscape software 3.8.2 was then used to visualize the network construction.

To identify potential medicinal targets, researchers (by Ran Cheng and Sheng-kai Wang) intersected melatonin and PCOS targets using the Venn Diagram tool. The PPI network of the gene lists was predicted using the STRING database with a moderate confidence interaction score (0.400) (36). The Cytoscape plugin “CytoHubba” was employed to identify hub genes using three topological analysis methods: Maximal Clique Centrality (MCC), Maximum Neighborhood Component Centrality (DMNC), and Neighborhood Component Centrality (MNC) (37).

DAVID (Database for Annotation, Visualization, and Integrated Discovery) equips researchers with an extensive array of functional characterization instruments, enabling them to uncover the biological relevance of substantial gene collections (38). We employed DAVID tools to conduct GO and KEGG pathway enrichment analyses to elucidate the functions of the potential medicinal targets as well as signaling pathways (by Ran Cheng and Sheng-kai Wang). The results of these analyses were visualized using a freely available online bioinformatics platform (https://www.bioinformatics.com.cn/).

The Protein Data Bank (39) supplied the 3D crystal structures of the target proteins in the PDB format. The PyMOL software (40)was used to eliminate water molecules and excess ligands. Subsequently, AutoDock tools (1.5.6) were utilized to incorporate protons and generate the PDBQT format of the proteins (41).

The 2D SDF structure of melatonin was obtained from PubChem (42). Energy minimization was performed using Chem 3D software (21.0.0) and the mol2 format was created (43). AutoDock tools (version 1.5.6) were then used to create the PDBQT format of melatonin.

Molecular docking with AutoDock Vina (40) was employed to assess melatonin binding and the hub pivot targets. Prior to molecular docking, proteins were represented as secondary structures. AutoDock Vina computed the binding affinity, while PyMol 2.5.4 (open-source version) was displayed to visually render the docking results (by Ran Cheng).

The configurations exhibiting the two leading scores derived through molecular docking were subsequently subjected to molecular dynamics (MD) simulations utilizing Gromacs 2022.3 software (44, 45). The preparation of small molecules involved pretreatment with AmberTools22 and hydrogenation using Gaussian 16 W to determine the response potential. The computed potential data were incorporated into the MD system’s topology file. Throughout the simulation, the conditions were kept constant at 300 K and 1 bar pressure. The simulation employed the amber99sb ildn force field, and Tip3p was chosen as the water model for solvation. To achieve a neutral total charge in the simulated system, an appropriate amount of Na+ ions was introduced. The energy of the simulation system was minimized through a two-step process: first, the steepest descent technique was employed, followed by equilibration procedures. These included 100000 steps each of NVT (constant volume and temperature) and NPT (constant pressure and temperature) equilibration, utilizing a coupling constant of 0.1 ps over a 100ps time period. After equilibration, a free molecular dynamics simulation was performed, comprising five million steps with a 2fs step size, amounting to 100ns. A post-simulation assessment was conducted using the built-in tools of the software. This analysis encompassed the calculations of RMSD, RMSF, and rotational radius of the protein for each amino acid trajectory, along with free energy (MMGBSA) computations and generation of free energy topographic maps, among other relevant data (by Ran Cheng). These structural predictions provide a theoretical basis for the potential mechanism of action, but they need to be comprehensively analyzed in conjunction with the expression changes at the transcriptional level.

In the Network Pharmacology section, topological analysis was performed using Cytoscape 3.8.2. The DAVID tools were applied for GO and KEGG enrichment analyses. AutoDock tools (1.5.6) were utilized to construct molecular docking simulations with PyMol 2.5.4, which were constructed using Gromacs 2022.3 software. Data evaluation was performed using GraphPad Prism 8 and SPSS version 25. As for the enrichment analysis, statistical significance was determined by a P -value < 0.05. In the section of experimental validation, results are expressed as mean ± SD. Each experiment was performed three times on an average. Data were first analyzed using the Kolmogorov–Smirnov test to verify the standard distribution. In cases where the data exhibited a normal distribution, statistical analyses were performed using either a paired Student’s t-test or one-way ANOVA, with subsequent Newman-Keuls post hoc testing. As for non-normally distributed data, the Kruskal–Wallis test followed by Dunn’s post-hoc test was applied. The threshold for statistical significance was established at P -value < 0.05(two-sided). This is denoted by *P < 0.05, **P < 0.01, and ***P < 0.001.

To examine how melatonin affects PCOS-like characteristics induced by circadian rhythm disruption, we administered melatonin injections to female Sprague–Dawley rats over an 8-week period while housing them in continuous darkness ( Figure 1A ). Rat-body-weights in the different groups were identical ( Figure 1B ). In the PCOS-like rat model subjected to constant darkness, we noted several characteristic features. These included irregular estrous cycles ( Figure 1C ), an increase in the quantity of abnormally large follicles, downregulation of ovarian corpora lutea ( Figure 2 ), and elevated serum testosterone levels ( Figure 1E ), which were consistent with previous studies (10, 26, 27). But the serum LH level remained unchanged ( Figure 1D ). The melatonin administration improved the disrupted estrus cycles and decreased the number of antral follicles in darkness rats. Besides, melatonin treatment also significantly alleviated the high serum testosterone level ( Figures 1E , 2 ).

A comprehensive set of 233 melatonin targets was compiled from PharmMapper (119 targets), SwissTargetPrediction (100 targets), DrugBank (17 targets), and TCMSP (7 targets). Among them, 10 duplicated targets were removed and subsequently illustrated using Cytoscape (version 3.8.2) ( Supplementary Figure S1 ).

Based on DisGeNET (262 targets), GeneCards (436 targets), and NCBI Gene (457 targets), 803 PCOS targets were identified. Utilizing Cytoscape’s analysis network feature in the string interaction network, we identified 103 significant targets for PCOS. These targets exhibited above-average values for a degree of centrality of 14.00000, betweenness centrality of 0.00348, and closeness centrality of 0.33856 ( Supplementary Figure S2 ).

By applying the Venn Diagram tool, researchers (Ran Cheng and Sheng-kai Wang) identified 37 potential targets related to melatonin and PCOS ( Figure 3A ). Subsequently, we utilized the STRING database to generate a 37-co-targets’ PPI-network ( Figure 3B ). The MCC, DMNC, and MNC scores were calculated by “CytoHubba” plugin, and the intersection of hub targets [AR, CYP19A1, PGR, IGF1R, MDM2, NR3C1, and CYP1A1] was obtained.

DAVID tools were employed to analyze the 37 common targets, with the top 10 Gene Ontology (GO) terms selected based on P-value and count frequency. Figure 3C displays the GO results, which were represented through a platform (free online: https://www.bioinformatics.com.cn/). The five most prominent biological processes (BP) were recognized in line with their count frequency and P-value, as shown in Figure 3D . These key BPs were found to be associated with melatonin treatment for PCOS: the intracellular steroid hormone receptor Signaling pathway, the positive regulation of the MAPK cascade, the steroid catabolic process, and the positive and negative regulation of transcription from the RNA polymerase II promoter. Importantly, four of the seven hub genes (AR, PGR, NR3C1, and MDM2) were enriched among the five leading BPs.

DAVID tools were employed to conduct a KEGG pathway enrichment analysis on the 37 therapeutic targets, identifying 18 pathways with significant associations, as depicted in Figure 4A . The targets enriched within individual pathways are shown in Figure 4B . These 18 pathways were organized into four functional modules, as shown in Figure 4C . Module 1 incorporated two pathways: circadian entrainment (hsa04713) and neuroactive ligand-receptor interaction (hsa04080). Module 2 encompasses six pathways: steroid hormone biosynthesis (hsa00140), ovarian steroidogenesis (hsa04913), prolactin signaling (hsa04917), estrogen signaling (hsa04915), GnRH secretion (hsa04929), and thyroid hormone signaling (hsa04919). Module 3 contained four pathways: diabetic cardiomyopathy (hsa05415), regulation of lipolysis in adipocytes (hsa04923), maturity-onset diabetes of the young (hsa04950), and endocrine resistance (hsa01522). Module 4 was composed of six pathways: the PI3K-Akt Signaling pathway (hsa04151), VEGF Signaling pathway (hsa04370), Rap1 Signaling pathway (hsa04015), FOXO Signaling pathway (hsa04068), oocyte meiosis (hsa04114), and progesterone-mediated oocyte maturation (hsa04914).

Seven hub genes sorted by MCC, DMNC, and MCN were identified for molecular docking analysis with melatonin ( Figure 5 ). Table 2 outlined the specifications of the molecular docking parameters. A stronger binding ability was indicated by a lower docking affinity. The binding strength between melatonin and its target molecules was substantial, with affinity values ranging from −7.5 to −5.4 kcal/mol. Melatonin may effectively treat PCOS by modulating the activities of these significant targets.

MD simulations were performed to verify the molecular docking simulation results for AR and CYP19A1. The RMSD curves for the AR/CYP19A1-melatonin complexes showed equilibrium after 60 ns, with average RMSD values of 0.2 nm and 0.23 nm ( Figure 6A ). Evaluation of the molecular stability using RMSD measurements in relation to the starting conformation ( Figure 6B ) showed that melatonin exhibited the smallest positional shifts when complexed with AR. The Rg values for the two complexes-maintained equilibrium with means of 1.85 nm and 2.25 nm ( Figure 6E ). The average SASA values remained consistent at 128 nm² and 210 nm² ( Figure 6F ). Throughout the 100 ns simulation, the count of hydrogen-bonding interactions for the AR/CYP19A1-melatonin complexes ranged from to 1–4 and 1-3, respectively ( Figures 6G, H ). The RMSF analysis revealed crucial flexible areas within the AR and CYP19A1 proteins, particularly in the vicinity of specific amino acid residues 200–40 for AR-melatonin; 790–800 for CYP19A1-melatonin ( Figures 6C, D, I, J ). The complex’s Gibbs free energy profiles exhibited clear minimum energy basins ( Figures 6A, B, E ), accompanied by hydrogen bond counts of 3 and 2, respectively ( Figures 6B–D, G–F ). Table 3 displayed the binding free energy of the complexes formed between AR/CYP19A1 and melatonin. The MD simulation generated movie-like files depicting the trajectories ( Supplementary Movie S1 . AR and Supplementary Movie S2 . CYP19A1).

As illustrated in Figure 7 , the minimum value of the free energy is indicated by the blue region, while the metastable conformation is represented by the cyan and green regions. A higher prevalence of blue coloration was observed for AR/CYP19A1-melatonin complexes, suggesting favorable binding energetics and conformational stability within the simulation timeframe. Molecular dynamics (MD) simulations thus indicated the potential for stable physical interaction between melatonin and the identified target proteins. While this computational finding supports the plausibility of melatonin binding to these targets, it is important to note that MD simulations primarily assess structural dynamics and binding stability, not functional consequences. The potential biological significance of these interactions, such as modulation of enzymatic activity or transcriptional regulation relevant to hyperandrogenism, remains hypothetical and requires empirical validation through functional assays.

As illustrated in Figure 8 , the mRNA expression of BMAL1 decreased significantly in the darkness contrasting to the control group (P < 0.05). Notably, Melatonin administration substantially reversed the mRNA expression of BMAL1 (P < 0.05). The mRNA expression of CLOCK and PER2 did not differ among the three groups (P < 0.05).

As illustrated in Figure 9 , the mRNA expression of MTNR1A, MTNR1B, and CYP19A1 decreased significantly in the darkness group compared to that in the control group (P < 0.05), whereas the mRNA expression of AR was markedly increased in the darkness group (P < 0.05). Compared to the darkness group, the mRNA expression of MTNR1A, MTNR1B, and CYP19A1 increased significantly in the melatonin group (P < 0.05), whereas the mRNA expression of AR was markedly reduced in the melatonin group (P < 0.05).