UF-related expression profiles for miRNA (GSE159959) and mRNA (GSE593) and GSE18096) were collected from the gene expression omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). The associated information from the above three datasets was shown in Figure 1A. Due to the different process platforms of the included datasets, we normalized these data using R’s limma package.

The normalized expression datasets were analyzed using the GEO2R (http://www.ncbi.nlm.nih.gov/geo/geo2r/) for comparison between UF tissues and matched myometrium tissues. The screening criteria were set as p-value <0.05 and |log2-fold change (FC)| > 1, the top 100 regulated miRNAs were selected as DEMs. With the same threshold, the upregulated and downregulated overlapping DEGs were obtained between GSE593 and GSE18096. The DEMs and DEGs in the respective profiles were visualized in volcano plots. Furthermore, heatmaps of DEMs and DEGs were visualized by TB tools, which is a toolkit that integrates various biological data-handling tools.

The online STRING database (https://string-db.org) was used to visualize the positive functional interaction between DEGs and the threshold (combined scores ≥0.4). Cytoscape software vision 3.8.2 was utilized to screen the upregulated and downregulated DEGs and construct the PPI network. Additionally, molecular complex detection (MCODE) v1.5.1.16 (Bader and Hogue, 2003) was used to identify hub genes with the following cut-off criteria: degree cut-off = 2, max. depth = 100, node score cut-off = 0.2, and k-Core = 2.

To further clarify the functional annotation of DEGs and hub genes, Gene Ontology (GO) functional annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis were conducted in the Sangerbox database (http://sangerbox.com/Tool). GO functional analysis was mainly classified into cellular components (CC), biological processes (BP), and molecular functions (MF).

In order to explore the target drugs of hub genes, the drug-gene interaction database (DGIdb, https://dgidb.genome.wustl.edu/) was used, which is the online tool for drug-gene predication. Finally, Cytoscape was used to construct a drug-gene network.

Recently, artificial intelligence brought revolutionary changes to the field of protein structural biology through convolutional neural networks. As an AI prediction tool, Alphafold (https://alphafold.ebi.ac.uk/), as an AI prediction tool, predicts highly accurate protein structures competitive with experiments based on their amino acid sequence. Therefore, despite the crystal structures of human BIRC5 and TYMS having yet to be solved in the RCSB protein data bank (PDB), we adopted Alphafold to build the structure of STAT6 through in silico modelling. From the amino acid sequence of human BIRC5 and TYMS (accession numbers: O15392 and P04818, respectively) from the UniProt-KB database (http://www.uniprot.org/) and predicted protein structures starting from the sequence through AlphaFold v2.0. The stereochemical quality of the predicted model was evaluated by the local distance difference test (LDDT) score in the AlphaFold database.

Molecular docking simulation is a technique to predict the first-rank orientation of small-molecules to macromolecule targets to form stabilized complexes. Associated small-molecules were both downloaded from the Pubchem database (https://pubchem.ncbi.nlm.nih.gov). The AutoDockTools and PyMol were applied to preprocess the input file, including hydrogenation and deletion of crystallographic water and ligands. Molecular docking between small-molecules and BIRC5, TGFBR2, and TYMS binding pockets was performed by Autodock Vina with default parameters. The predicted binding interaction geometries of irinotecan (aka irinotecan hydrochloride) and BIRC5 or TYMS were visualized and the docking affinity between small-molecules and protein targets was scored. The optimal docking conformation and related results were analyzed by PyMol.

MiRWalk3.0 (http://zmf.umm.uni-heidelberg.de/apps/zmf/mirwalk3/) is an online tool to predict the downstream target genes of DEMs. Then, differentially expressed target genes (DETGs) were obtained by overlapping the predicted genes and DEGs for further analysis. In addition, the miRNA-target gene regulatory networks were visualized by Cytoscape. Among the DEMs in this network, we selected the DEMs with the highest number of regulatory genes ranking top two as the hub DEMs.

To further analyze the function of hub DEMs, GO and KEGG pathways were performed using miRPath v3.0 (http://snf-515788.vm.okeanos.grnet.gr/). An adjusted p < 0.05 was considered significantly enriched.

It has been confirmed that TF-mediated transcriptional regulation and miRNA-mediated post-transcriptional regulation are tightly coupled, which implies that they share a common regulatory mechanism (Lewis et al., 2005). Transcription factors and miRNAs are known to regulate each other as well as their target genes. Therefore, further research on the regulatory relationship between TFs and miRNAs can help elucidate the underlying mechanism of tumorigenesis and progress. Based on the hub DEMs, we predicted the potential TFs via the iRegulon plug-in (Lewis et al., 2005). Relevant parameters were as follows: enrichment score threshold = 5.0, ROC threshold for AUC calculation = 0.05, rank threshold = 5,000, minimum identity between orthologous genes = 0.05, and maximum FDR on motif similarity = 0.001. The Cytoscape software was subsequently utilized to visualize miRNA-TF-target gene regulatory networks.

Ten cases of UF tissue and matched normal myometrial tissue samples were collected from patients who underwent myomectomy in Xiangya Hospital from January 2021 to January 2022, aged 36–49 years old. All patients had regular menstruation without other complications, and they did not take hormone drugs within 3 months before surgery. All samples were taken with the consent of the patients, and informed consent was signed. Consent was granted by all women included in the study, and informed consent was signed. Tissues were stored at −80°C immediately after harvest until further use. The pathological diagnosis of the tissue was confirmed by at least two independent pathologists.

The tissues were fixed for 24 h and dehydrated with gradient alcohol. Paraffin-embedded tissues were prepared and subsequently cut into 5 µm slices, which were then stained using the H&E staining kit (Beijing Solarbio Science and Technology Co., Ltd., Beijing, China). The pathological changes in UF were observed and captured with five random images under a light microscope (Leica DMi6-M, Leica Microsystems Co., Ltd.).

The human UF tissue was washed three times with PBS solution and minced into small pieces (3–6 mm³) with sterilized eye scissors. To obtain single-cell suspensions, the tissue fragments were incubated with Dulbecco’s modified eagle media (DMEM, Hyclone, United States) medium supplemented with 4 mg/ml type I collagenase (Sigma-Aldrich, United States) and 1 mg/ml DNase I (Sigma-Aldrich, United States) at 37°C in a humidified incubator with 5% CO2 for 2 h. Subsequently, the suspensions were passed through a 40-μm cell strainer (Corning, United States) to eliminate undigested tissue. The primary cells were cultivated in DMEM medium containing 10% fetal bovine serum (FBS, Gibco, United States), 100 U/ml penicillin, and 100 μg/ml streptomycin at 37°C in a 5% CO2 incubator. Cultured human primary cells were passaged every 4–6 days and analyzed at the third passage for further study.

To explore the function of DEMs in leiomyoma cells, the primary cells were seeded in 6-well plates at an appropriate density. When cells were in log-phase growth, cell transfection was performed with negative control (NC) or hub DEM mimics or inhibitors (50 nM) using Lipofectamine 3000 (Invitrogen, United States). The sequences were synthesized by GenePharma (Shanghai, China).

RNA was extracted from tissues using the TRIZOL method (TransGen Biotech, Beijing, China) and the RNA was reverse transcribed into cDNA according to the instructions of the NovoScript® II reverse transcriptase (Novoprotein, Shanghai, China). MiRNA was extracted using the miRNA extraction kit (Tiangen, Beijing, China) and reversely transcribed into cDNA through the polyA tailing method. RT-PCR was carried out by using SYBR® Green MasterMix (Takara, Japan) for hub DEGs and miRNA universal SYBR qPCR Master Mix (Vazyme Biotech, Nanjing, China) for hub DEMs in an ABI PRISM® 7500 Sequence Detection System (Applied Biosystems Inc: Foster City, CA). RT-PCR reaction system: For hub DEGs, 45 s at 95°C, 40 cycles for 30 s at 95°C, 45 s at 60°C, 60 s at 72°C and a final cycle for 10 min at 72°C; For hub DEMs, 5 min at 95°C, 40 cycles for 10 s at 95°C, 30 s at 60°C. Fluorescence values were collected. Replicate wells were set up for each sample. GAPDH and U6 were used as the internal reference for hub DEGs and DEMs detection, respectively, and the experiment was repeated three times independently. The mean CT value was calculated. Relative quantification (ΔΔCT method) was applied for semi-quantitative analysis, and the target gene expression was represented by 2^−ΔΔCT values. The sequences of the primer were synthesized by Sangon Biotech (Shanghai, China).

Cells were transfected for 48 h and then inoculated into 96-well plates at a density of 4 × 10³/well. The cell proliferation was measured at 1–3 days, respectively. At each time point, 10 μl of CCK-8 solution (New Cell and Molecular Biotech, Jiangsu, China) was added to each well and incubated in the incubator for 2 h. Finally, the absorbance at 450 nm, as the optical density (OD) of the samples, was detected by spectrophotometry at different time points.

After 24 h of transfection of mimic or inhibitor, we trypsinized and resuspended to obtain a single cell suspension. The upper layer of the Transwell chamber was added with FBS-free suspension at a cell density of 3 × 10⁴/well. The lower layer was added with 20% FBS medium to drive cell migration. After 24 h of culture, we used cotton swabs to remove the nonmigrated cells in the upper chamber, and migrated cells in the lower chamber were fixed and stained with 0.1% crystal violet. Five randomly selected microscopic views were recorded through a light microscope and analysed by ImageJ software to count the number of cell migrations.

To determine the effect of hub DEMs on the apoptosis of leiomyoma cells, cells were harvested and centrifuged for annexin V-FITC/PI apoptosis detection (Boster Biotech Co., Ltd., Wuhan, China) after 24 h of transfection. After being washed twice, transfected cells were incubated with 5 µl of 488-annexin V and 5 µl of PI for 15 min at room temperature and protected from light. Subsequently, the flow cytometer (BD, United States) was applied to detect the apoptosis rates in each group. The apoptosis rate = sum of the early and late apoptotic rates.

Statistical analysis was carried out using SPSS 16.0 statistical software. Data was expressed as mean ± SD, and the differences between the two groups were compared using the independent sample t-test. The differences between the two groups were compared using the independent sample t-test, and the differences were considered statistically significant at p < 0.05.

Before screening DEGs and DEMs in UF, the three independent expression datasets processed on different platforms were normalized. Based on the mentioned threshold, differential expression analysis was applied to obtain DEGs and DEMs. 1185 and 2054 DEGs and 229 DEMs were identified in GSE593, GSE18096, and GSE159959, respectively. To shrink the scope of DEM screening, DEMs whose |log2 FC| ranked in the top 100 were chosen as candidate DEMs. As shown in Figures 1B–D, 2A–C, volcano plots and heatmaps exhibited the overall distribution of DEGs and DEMs for UF versus normal myometrium. Furthermore, a total of 397 candidate DEGs were identified by Venn diagram (Chen and Boutros, 2011) analysis through the intersection of these two mRNA datasets (Figure 2D). Excluding DEGs with inconsistent expression trends in these two datasets, 373 genes containing 192 upregulated and 181 downregulated genes were obtained as candidate DEGs.

To identify potential interactions among 373 candidate DEGs, the PPI network was constructed using the STRING database and visualized by Cytoscape (Figure 3A). The network included 297 nodes (genes) and 1,111 edges (interactions). Furthermore, hub DEGs was validated by MCODE, and then we obtained 8 cluster modules. As indicated in Figure 3B, the highest-rated module included 24 nodes (genes) and 99 edges (interactions). The 24 hub DEGs included 19 upregulated and 5 downregulated genes.

To further reveal the vital function and biological pathways of DEGs and hub DEGs. We performed GO and KEGG analysis of 373 candidate DEGs and 24 hub DEGs. As illustrated in Figures 4A–D, for the BP-correlated category, DEGs were mainly enriched in biological adhesion, locomotion and cell migration; for CC-correlated category, DEGs were primarily involved in collagen containing extracellular matrix (ECM), an intrinsic component of the plasma membrane and cell surface; and for the MF-correlated category, DEGs were principally concerned with extracellular matrix structural consitituent, growth factor binding and signaling receptor binding. Furthermore, KEGG pathway results demonstrated that DEGs primarily participated in ECM-receptor interaction, focal adhesion, and the PI3K-AKT pathway.

As for hub DEGs, BP-enriched terms included external encapsulating structure organization and collagen fibril organization; CC-enriched terms were involved in supramolecule complexes and collagen trimers; MF-enriched terms referred to extracellular matrix structural consitituent and structural molecule activity. KEGG results contained ECM-receptor interactions and pathways in cancer (Figures 5A–D).

To seek effective target-drugs for hub DEGs, drug-gene interactions, including 75 potential drugs for UF, were obtained through the DGIdb database. For ease of reading, drug-gene networks were visualized by Cytoscape (Figure 6). However, the underlying mechanism between most potential drugs and hub DEGs remains unrevealed.

Molecular docking was applied to further screen potential drugs and elucidate possible molecular mechanisms. The crystal structure of human BIRC5 and TYMS was predicated by Alphafold v2.0 based on AI technology. As shown in Figures 7A,B, the stereochemical quality of the predicted structure was evaluated by the LDDT score. Autocdock Vina was used for drug-protein molecular docking to screen the optimal potential target drugs. The Affinity score can be used as a standard to judge the merits of docking. The high absolute value of the score indicated stronger binding between small molecules and proteins. The docking scores of potential drugs were presented in Figures 7C–F and suggested that irinotecan had the strongest binding affinity towards BIRC5 (7.0–8.6 kcal/mol, |interval score|) and TYMS (7.7–10.0 kcal/mol, |interval score|). The perfect conformation exhibited that irinotecan interacted with residues of BIRC5 (THR21 and HIS45) and TYMS (ARG50, ILE108, ASN226 and TYR258) through hydrogen bonds. After screening, irinotecan was found to be the optimal drug for intervention of both BIRC5 and TYMS and was used to treat UF patients in low doses.

Based on UF and matched normal myometrial tissues in patients at our hospital, we validated the expression level of hub DEGs through RT-PCR (Figures 8A, B). Compared with normal myometrium, the expression levels of ASPM, SDC1, SNAI2, GINS1, CCNB2, KIAA0101, BIRC5, BUB1, TYMS, FOXM1, DTL, KIF20A, COL3A1 were increased, which was consistent with microarray analysis results. Moreover, the expression levels of TGFBR2, ITGA6, CYR61, ITGB4, COL11A2, COL4A5, COL4A6, and COL4A4 were decreased in comparison to normal myometrium. There were no significant differences in COL13A1 and COL5A2 expression.

Recent research indicated that tumorigenesis and development were regulated by miRNAs through binding with pairing bases of 3′UTR, resulting in target mRNAs degradation. Therefore, predication of target mRNAs can contribute to revealing the function and internal mechanism of DEMs. MiRWalk 3.0 was used to predict the target genes of DEMs. A total of 1,507 genes were identified as the potential target genes of the DEMs. After intersecting with 373 candidate DEGs, 28 overlapping genes were identified as DETGs for further analysis (Figure 9A). The miRNA-target gene regulatory network was visualized using Cytoscape. As indicated in Figure 9B, triangle represents DEMs, orange ellipses represent upregulated DETGs, and purple ellipses represent downregulated DETGs. The regulatory network included 21 DEMs, 16 upregulated DETGs, and 12 downregulated DETGs. The DEMs were ranked by the number of regulated DETGs, and the top two were identified as hub DEMs. In vitro and vivo experiments validated the upregulated hsa-miR-181b-5p and hsa-miR-381-3p as well as five screened DETGs.

GO and KEGG pathway analysis were performed to elucidate the function of hub DEMs through miRPath v3.0. As illustrated in Figures 9C–F, these hub DEMs were enriched in ion binding, cellular nitrogen compound metabolic process, organelle, biosynthetic process and transcription, DNA-templated in the BP terms; Organelle and cellular component in CC terms; Ion binding, nucleic acid binding transcription factor activity and molecular function in MF terms. Furthermore, KEGG results indicated that hsa-miR-181b-5p and hsa-miR-381-3p involved in amphetamine addiction and TGF-β signaling pathway and proteoglycans in cancer respectively.

The transcriptional regulatory relationships among these miRNAs, TFs, and target genes were identified through the iRegulon plug-in. Then, miRNA-TF-target gene networks were established through Cytoscape software. In the transcriptional regulatory network, hsa-miR-181b-5p consisted of 2 TFs (STAT3 and ONECUT1) and hsa-miR-381-3p consisted of 2 TFs (FOXA1 and KAT2A) (Figures 10A, B).

In order to assess the expression of miR-181b-5p and miR-381b-3p in UF, RT-PCR was performed in UF tissues compared to normal myometrial tissues (n = 5). Consistent with our analysis of hub DEMs, we found that miR-181b-5p and miR-381b-3p were markedly upregulated in UF tissues compared to normal myometrium tissues (Figure 10C).

In the functional verification section, the hub DEMs, mimics, and inhibitors were transfected into UF cells. The CCK-8 assay was used to examine the effects of hsa-miR-181b-5p and hsa-miR-381-3p on cell proliferation. As demonstrated in Figures 11A,B, overexpression of both hub DEMs enhanced the proliferation of UF cells compared with NC mimics, while hsa-miR-181b-5p and hsa-miR-381-3p inhibitors could reverse the above results.

Although UF is a benign tumor with poor metastatic ability, it may locally migrate. Therefore, the evaluation of cell migration still reflects hub DEMs’ function. Many studies have also shown that UF cells have the metastatic ability. After transfection with hsa-miR-181b-5p and hsa-miR-381-3p mimics, the number of migrational cells was higher than NC mimic group. On the contrary, the metastatic ability of transfected cells with hsa-miR-181b-5p and hsa-miR-381-3p inhibitors was remarkably reduced, with fewer cells crossing through the membranes. Above, results suggested that there were positive correlations between the metastatic ability of UF cells and hsa-miR-181b-5p and hsa-miR-381-3p expression (Figure 11C)

Although surgical excision is the undisputed treatment for benign tumors, apoptosis of tumors induced by non-invasive procedures can effectively avoid operative trauma and tends to improve the outcome. MiRNA therapeutics provides a novel horizon for tumor treatment. To verify whether hub DEMs’ silence can induce UF cells’ apoptosis, we transfected NC, hsa-miR-181b-5p, and hsa-miR-381-3p inhibitors to UF cells. As indicated in Figure 11D. Downregulation of hsa-miR-181b-5p and hsa-miR-381-3 could significantly contribute to cell apoptosis through comparison of cell apoptosis rates in each group by flow cytometry.