Integrated bioinformatic analysis of dysregulated microRNA-mRNA co-expression network in ovarian endometriosis

D‐L and Y‐L: study conception and design. Y‐L, LG‐J, C‐L, T‐L, Z‐Z, CM‐X, SH‐L, YD‐Z and J‐R: data acquisition. Y‐L, LY‐F and LG‐J: analysis and data interpretation. Y‐L and LY‐F: drafting of the manuscript. All authors contributed to the article and approved the submitted version.

This study was funded by Beijing Municipal Excellent Talents Foundation (2018000021469G282) and Beijing Obstetrics and Gynecology Hospital, Capital Medical University Program (FCYY202010).

The mRNA‐sequencing of 13 tissue samples generated 1 427 710 818 total clean reads and the FastQC quality test showed that 1 416 264 272 reads reached Q30 for further analyses. In all, 1 387 396 603 reads (97.96%) were uniquely aligned to hg38, and 16 619 mRNAs were identified. The detailed filtering and mapping data are listed in Table  S2 . Compared with normal endometria, 1467 DEGs were identified in EC tissues (Figure  2A ) and 224 DEGs in EU tissues (Figure  2B ). A total of 1524 DEGs were identified in EC tissues compared with EU tissues (Figure  2C ). Specifically, six mRNAs gradually increased from normal endometria to EU, compared with EC tissues: PIP4P2, MPV17L, KPNA5, ZNF420, EIF5A2 and ESD. In contrast, 24 mRNAs were gradually decreased (Figure  2D–F , Table  S3 ). Differentially expressed genes (DEGs) in normal endometria of three volunteers and paired ectopic/eutopic endometria of five patients with OEM. Volcano plot of DEGs in (A) EC ( n  = 5) compared with normal endometria ( n  = 3); (B) in EU ( n  = 5) compared to normal endometria ( n  = 3); (C) in EC ( n  = 5) compared to EU ( n  = 5). Venn diagrams of (D) down‐regulated (D) and (E) up‐regulatedDEGs in normal, eutopic and ectopic endometria. (F) Heat map of gradually increased and decreased DEGs from normal to EU, and then in EC tissues. The relative expression is depicted by color: red indicates up‐regulation; green indicates down‐regulation. EC, ectopic endometrium; EU, eutopic endometrium; Normal, normal endometrium. The small RNA‐seq generated 473 156 630 reads with an average of 36 396 663 reads per sample. The FastQC quality test showed that 467 372 834 reads had a Q‐score ≥30; these reads were considered in further analyses (Table  S4 ). The length of 119 103 727 reads were within 18–35 nucleotides and uniquely aligned to hg38. In all, 2756 miRNAs were identified after mapping to the miRBase database (v22.0) and miREvo. Compared with normal endometria tissues, 282 DEMs were identified in EC tissues (Figure  3A ) and 99 DEMs in EU tissues (Figure  3B ). A total of 220 DEMs were identified in EC tissues compared with EU tissues (Figure  3C ). Importantly, expression of seven miRNAs (miR‐141‐3p, miR‐144‐3p, miR‐363‐3p, miR‐484, miR‐577, miR‐767‐5p and miR‐96‐5p) were gradually decreased from normal endometria to EU, compared with EC tissues. In contrast, miR‐374a‐3p, miR‐382‐3p, miR‐493‐5p and miR‐592 were gradually increased (Figure  3D–F , Table  S5 ). Differentially expressed miRNAs (DEMs) in normal endometria of three volunteers and paired ectopic/eutopic endometria of five patients with OEM. Volcano plot of DEMs in (A) EC ( n  = 5) compared with normal endometria ( n  = 3); (B) in EU ( n  = 5) compared with normal endometria ( n  = 3); (C) in EC ( n  = 5) compared with EU ( n  = 5). Venn diagrams of down‐regulated (D) and up‐regulated (E) DEMs in normal, eutopic and ectopic endometria. (F) Heat map of gradually increased and decreased DEMs from normal to EU, and then in EC tissues is presented. The relative expression is depicted by color: red indicates up‐regulation; green indicates down‐regulation. EC, ectopic endometrium; EU, eutopic endometrium; Normal, normal endometrium. The miRNA and mRNA‐sequencing data of EC and EU tissues from five patients were normalized with sequencing datasets ( GSE105764 for mRNA sequencing and GSE105765 for miRNA sequencing) of EC and EU tissues from eight patients (Figure  S1 A,B). In these 13 paired EC and EU samples, 1237 DEGs were identified in EC tissues, of which 656 genes were up‐regulated and 581 down‐regulated (Figure  4A , Table  S6 ). Importantly, eight mRNAs (DYNLRB1, TRPM4, SEPT9, CAPNS1, RASEF, DSP, CREB3L1, MARVELD3) were gradually decreased and one mRNA (MPV17L) was gradually increased from normal endometria to EU compared with EC tissues (Figure  2F ) in these 1237 DEGs (Figure  4B ). A total of 348 DEMs were found (Figure  4C ), including 190 up‐regulated miRNAs and 158 down‐regulated miRNAs, between EC and paired EU. Among them, miR‐202‐5p was the most up‐regulated (Log 2 FC = 10.20) and miR‐375‐3p the most down‐regulated (Log 2 FC = −7.73) (Table  S7 ). Five miRNAs (miR‐141‐3p, miR‐363‐3p, miR‐577, miR‐767‐5p, miR‐96‐5p) were gradually decreased and two miRNAs (miR‐493‐5p, miR‐592) gradually increased from normal endometria to EU compared with EC tissues (Figure  3F ) in these 348 DEMs (Figure  4D ). Differentially expressed genes (DEGs) and miRNAs (DEMs) in EC compared with EU tissues of OEM patients in integrated sequencing data. The miRNA and mRNA‐sequencing data of EC and EU tissues from five patients were normalized with sequencing datasets ( GSE105764 for mRNA sequencing and GSE105765 for miRNA sequencing) of EC and EU tissues from eight patients. (A) Volcano plot of all DEGs in paired EC ( n  = 13) compared with EU tissues ( n  = 13). (B) Heat map of nine mRNAs that gradually increased/decreased in Figure  2F in paired EC ( n  = 13) compared with EU tissues ( n  = 13). (C) Volcano plot of all DEMs in paired EC ( n  = 13) compared with EU tissues ( n  = 13). (D) Heat map of 12 miRNAs that gradually increased/decreased in Figure  3F or the top five miRNAs in node degree analysis in paired EC ( n  = 13) compared with EU tissues ( n  = 13). The relative expression is depicted by color: Red indicates up‐regulation; green indicates down‐regulation. EC, ectopic endometrium; EU; eutopic endometrium. In all, 348 DEMs and 1237 DEGs identified in the previous steps were considered for the miRNA/gene binary interaction and regulatory network analysis. According to the standard of absolute PCC >0.6, correlation analysis revealed 1850 paired miRNA‐mRNA regulatory interactions between 87 miRNAs and 482 mRNAs. Although miRNAs usually suppress target mRNAs, positively paired miRNA‐mRNA correlations are occasionally observed under certain conditions. 21 , 35 For this reason, we identified and constructed the miRNA‐mRNA regulatory network of all differentially expressed miRNA‐mRNA pairs, including 934 negative correlating pairs and 916 positive correlating pairs. To identify key miRNAs in the integrated networks, we also calculated the node degree of each involved miRNAs via directly counting the links of one node. 36 As a result, the node degrees of the selected miRNAs ranged from 1 to 72; the top five key miRNAs in the miRNA‐mRNA network of OEM were miR‐17‐5p, miR‐20b‐5p, miR‐106a‐5p, miR‐20a‐5p and miR‐216a‐3p (Figure  4D , Table  S8 ). The key miRNA‐mRNA regulatory network is presented in Figure  5 . The integrated ovarian endometriosis‐related key miRNA‐mRNA co‐expression regulatory network: (A) miRNA‐mRNA co‐expression regulatory network of 12 important miRNAs and their target mRNAs in integrated sequencing data; (B) miRNA‐mRNA co‐expression regulatory network of nine important mRNAs and their corresponding miRNAs in integrated sequencing data. Red indicates up‐regulation; green indicates down‐regulation; orange indicates no significant change. Circles indicate mRNAs and triangles indicate miRNAs. Functional enrichment analysis of the target mRNAs in the key miRNA‐mRNA regulatory network (Figure  5 ) provided us with a general clue regarding their functional roles in OEM. GO analysis revealed that the top 10 GO biological process (GO‐BP) categories were mainly enriched in movement of cell or subcellular component, cell adhesion, regulation of cell communication, locomotion and extracellular matrix organization (Figure  6A ). The GO cellular component (GO‐CC) category was mainly involved in cell–cell junction, cell–cell adherens junction, cation channel complex and transmembrane transporter complex (Figure  6B ). The GO molecular function (GO‐MF) category was mainly involved in cAMP response element binding, chromatin binding, nucleoside binding, calcium channel activity, calmodulin binding (Figure  6C ). KEGG pathway analysis was significantly associated with Tight junction, PI3K‐Akt signaling pathway, AMPK signaling pathway, human papillomavirus infection and TNF signaling pathway (Figure  6D ). Gene Ontology terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis of miRNA‐targeted mRNAs that involved in the key regulatory network. (A) Biological processes. (B) Cellular components. (C) Molecular functions. (D) KEGG pathways.

Because the etiology of OEM is unknown, the treatment options available are limited, with a high recurrence rate, which adds to its large socioeconomic impact. 37 OEM refers to the growth of cells derived from the endometrium in ovaries. DEGs 15 , 16 , 18 and DEMs 19 , 20 between EC and EU ovary tissues have been individually characterized in several microarray studies investigating the molecular mechanisms of OEM. Moreover, integrated bioinformatic analysis of these studies using a non‐biased approach recently confirmed that cell adhesion, inflammatory response and EMT were associated to the development of OEM. 17 , 38 , 39 However, the joint analysis of the miRNA‐mRNA regulatory network involved in OEM is rarely reported. In the present study, we detected miRNA and mRNA expression profiling of EC and EU tissues using high‐throughput sequencing techniques and performed integrated bioinformatic analysis of our data and sequencing data obtained by Zhao et al. 21 Several studies confirmed that a large number of miRNAs are crucial for the establishment and progression of OEM via directly regulating genes involved in angiogenesis, inflammation and immune regulation. 40 , 41 , 42 , 43 Thus, miRNA‐mRNA network construction is essential for investigating their functional roles in OEM. Our study revealed 348 DEMs and 1237 DEGs by integrated bioinformatic analysis. The large differential expression profiling was refined by PCC to 87 functionally DEMs and 482 target mRNAs. The constructed miRNA‐mRNA network included 1850 paired miRNA‐mRNA regulatory interactions and some miRNAs (eg miR‐17‐5p, miR‐20a‐5p, miR‐34c, miR‐449a, miR‐141‐3p) are similar to previous reports concerning the development of endometriosis. Here, we discuss these key regulators in the network, with the hope of providing valuable information for understanding the molecular mechanisms of OEM. MicroRNA‐17‐5p and miR‐20a‐5p are two members of the miR‐17‐92 cluster, a down‐regulated cluster in endometriotic lesions. 44 , 45 In the present study, the downregulated miR‐17‐5p/20a was directly connected to 72 and 70 corresponding targets, respectively. Several potential targets of miR‐17‐5p/20a, eg BCL2, CDKN1A/p21, HIF‐1α, VEGF‐A, TGF‐β and IL‐8, are associated with diverse biological events contributing to OEM development, eg proliferation, apoptosis, inflammation, EMT and angiogenesis. 46 Three members of the miR‐34/449 family (miR‐449a, miR‐449b and miR‐34c) were decreased in ectopic endometrial tissues, similar to previously published results. 47 Although the miR‐34/449 family was found to be associated with cell cycle, inflammation and actin network formation, 48 its function in endometriosis is still largely unclear. MicroRNA‐141‐3p has been reported to inhibit EMT in endometriosis via target inhibition of the TGF‐β1/SMAD2, Notch and Hippo signaling pathways. 49 , 50 , 51 , 52 In addition, miRNA‐141‐3p belongs to the miR‐200 family, which includes miR‐200a, miR‐200b, miR‐200c, miR‐141 and miR‐429, all key negative regulators of EMT. 53 In the present study, it is worth noting that all the miR‐200 family members were remarkably decreased in EC tissues compared with EU tissues (Table  S7 , Figure  4D ), and were also included in the miRNA‐mRNA network. We also observed a significantly gradual decrease of miR‐141‐3p from normal endometria to EU compared with EC tissues (Figures  3F and 4D ). As EMT is important for the development of endometriosis, 17 , 54 these results highlight the importance of miR‐200 family members in the pathogenesis of OEM. Endometriosis is a benign disease but it has invasive and metastasizing characteristics, and can lead to tumor formation in a small number of cases (0.5%–1%). 55 Aberrant activation of the RAS/MAPK and PI3K‐Akt signaling pathways is an inherent feature and is a potential driver contributing to endometriosis pathology outside of the context of cancer; 56 these pathways were enriched in KEGG analysis in the present study (Figure  6 ). The three members of miR‐106a‐363 cluster (miR‐363, miR‐20b and miR‐106a) were down‐regulated in ectopic lesions according to our results. More importantly, all of these miRNAs were reported to inhibit the activation of the MAPK signaling pathway in coronary arterial endothelial cells, 57 papillary thyroid carcinoma cells 58 and endometrial epithelial cells, 59 respectively. Although no study has reported the function of the miR‐106a‐363 cluster in endometriosis, miR‐363 could serve as tumor suppressor in different cancers 60 and was gradually decreased from normal endometrium to EU, and then in EC tissues (Figures  3F and 4D ). MicroR‐20b contained the highest node degree (connected to 72 targets) in the identified DEMs and showed a 6.9‐fold decrease in EC compared with EU. We therefore speculate that the miR‐106a‐363 cluster may be another important set of miRNAs in OEM. The expression of miR‐183/96/182 cluster members, including miR‐182‐5p, miR‐183‐5p and miR‐96‐5p, was previously reported to be down‐regulated in endometriosis. 19 , 61 , 62 In addition, functional analysis revealed that miR‐182‐5p inhibits the proliferation, migration and invasion of endometrial stromal cells in endometriosis via direct targeting of RELA of the NF‐κB signaling pathway. 63 MicroRNA‐96‐5p blocks the TGF‐β/SMAD signaling pathway via targeting of TGFBR1 and reverses the TGF‐β1‐induced EMT in endometriosis. 62 MicroRNA‐183‐5p suppresses the invasiveness of endometrial stromal cells via targeting of ITGB1, a critical factor for cell adhesion and invasiveness. 61 In the present study, these three miRNAs were decreased in EC compared with EU (Figure  4D , Table  S7 ), and miR‐96‐5p showed gradual decrease from normal endometrium to EU compared with EC tissues (Figure  3F ). These findings indicate that aberrant down‐regulation of the miR‐183/96/182 cluster is an important molecular event of pathogenesis and the development of endometriosis. Besides these decreased miRNA clusters, we also identified some novel increased miRNAs: namely, miR‐382‐3p, miR‐592, miR‐577 and miR‐493‐5p (Figures  3F and 4D ). Specifically, miR‐493‐5p showed a gradual increase from normal endometrium to EU compared with EC tissues in integrating data. Intriguingly, previous studies have confirmed that miR‐494‐5p suppresses proliferation and induces apoptosis of different cancer cells by direct targeting of multiple target genes (IGF2, HMGB1, AKT2, STK38L, HMGA2, ETS1, PHLPP2 and E2F5). 64 , 65 , 66 Thus, miR‐493‐5p may play a feedback role in OEM. MicroRNA‐382‐3p has dual functions, acting either as tumor suppressor or as oncogenic miRNA. On the one hand, miR‐382‐3p inhibits angiogenesis, proliferation, invasion and migration via targeting of VEGFA, EZH2, STAT1 and ZIC4 in various cancers. 67 , 68 But miR‐382‐3p is also a responsive miRNA of TGF‐β1, the well‐known inducer of EMT and chronic inflammation. 69 This unexpected finding requires further validation to identify the underlying mechanism of miR‐382‐3p and miR‐493‐5p in endometriosis. Although we constructed the miRNA‐mRNA network by integrating our high‐throughput sequencing data with datasets obtained by Zhao et al., 21 the limitations of this study should also be mentioned. First, OEM is a complex disease that is influenced by multiple factors, such as the menstrual cycle and CA125 levels, and this might impact the expression profiling. Secondly, the present results are preliminary and descriptive. The identified miRNA clusters and their direct targets need to be holistically validated in OEM by other experiments, such as qRT‐PCR, cross‐linking immunoprecipitation, in situ FISH or functional experimental validation. Thirdly, some direct targets of miRNAs identified in previous reports were not included in our miRNA‐mRNA network. These differences may be the result of the typical Pearson correlation coefficient utilized in this study, which might have led to differences compared with other computational methods (eg Bayesian network learning, regression models and causal inference) and the inevitable loss of information.

We provide the miRNA‐mRNA co‐expression network in OEM using integrated bioinformatic analysis with a non‐biased approach. The constructed miRNA‐mRNA network included 87 functionally DEMs, 482 target mRNAs and 1850 paired miRNA‐mRNA regulatory interactions. The present study confirms and extends the significance of prior studies, thereby defining the crucial roles of miRNA clusters, including miR‐17‐92 cluster, miR‐200 family, miR‐106a‐363 cluster and miR‐183/96/182 cluster, in the pathogenesis of OEM. The identified interactions could also have implications for targeted therapeutic strategies for OEM. Further in‐depth functional studies are encouraged to confirm our results.

Endometriosis is a frequently occurring gynecological disease estimated to affect 10%–15% of reproductive age females 1 and around to 20%–48% of infertile women. 2 The symptoms of endometriosis mainly include pelvic pain, dysmenorrhea and infertility, although asymptomatic cases do arise. 3 Ovarian endometriosis (OEM) is the most common type compared with the other two well‐recognized phenotypes superficial peritoneal lesions (SUP) and deep infiltrating endometriosis (DIE), accounting for 17%–44% of all endometrioses. 4 , 5 As a chronic estrogen‐dependent inflammatory disease, 6 the main therapeutic option for OEM is laparoscopic conservative surgery followed by hormonal suppression. 7 According to previous reports, the recurrence rate remains to be 21.5% at 2 years and 40.0%–50.0% at 5 years after primary surgery. 4 , 8 With the potential to inhibit ovulation, menstruation and downstream inflammation, postoperative hormonal suppression is effective in reducing recurrence, and improving patient‐reported pain and pregnancy rate. 9 , 10 , 11 Specifically, gonadotropin‐releasing hormone agonist (GnRHa) plus add‐back (norethindrone acetate/estrogens) treatment has been reported to be effective in 61% of cases. 12 However, some cases can be refractory to treatment with hormonal suppression such as combined hormonal contraceptives (CHC), progestins and GnRHa. 13 Thus, there is an urgent need to clarify the etiology and pathophysiology of OEM. Several theories have been suggested to be involved in OEM, including the most well‐accepted retrograde menstruation hypothesis, Müllerian metaplasia, ectopic implantation, proliferation of endometrial stem cells and bone marrow progenitors, genetic and epigenetic factors, immune and inflammatory factors. 3 To provide a molecular basis for these theories, several studies have been conducted to identify the gene 14 , 15 , 16 , 17 , 18 and microRNA (miRNA) 19 , 20 expression profiles of the ectopic endometrium (EC) and the eutopic endometrium (EU) from OEM patients. Aberrant miRNA and mRNA expression associated mainly with cell adhesion, inflammatory response, epithelial‐mesenchymal transition (EMT) and extracellular exosome, were confirmed to be involved in the development of OEM. Specifically, Zhao et al. 21 recently characterized the integrative miRNA‐mRNA co‐expression networks and found that some key regulators including the miR‐449 and miR‐106a‐363 cluster, miR‐200 family, miR‐182/183, FOX family, GATA family and E2F family, as well as CEBPA, SOX9 and HNF4A may play important roles in the pathogenesis of OEM. However, due to heterogeneity of independent experiments as a result of variations in tissue or specimens and/or different data processing methods (microarray or sequencing), the identification of these differentially expressed genes (DEGs)/miRNAs (DEMs) is inconsistent. In the present study, with the aim to understand better the molecular basis of the pathogenesis of OEM, high‐throughput sequencing techniques were used to detect miRNA and mRNA expression profiling in normal endometrium, paired EC and EU tissues. The miRNA‐mRNA regulatory network involved in the development of OEM was constructed by integrating our data with the sequencing data GSE105764 and GSE105765 21 with a non‐biased approach and deep bioinformatic analysis, including identifying DEGs and DEMs, gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment. The flow chart of the study design is presented in Figure  1 . Integrated bioinformatic analysis of the miRNA‐regulatory network can help us to understand the mechanism of OEM at the molecular level and enable the discovery of potential candidate diagnostic biomarkers or therapeutic targets in endometriosis. Flow chart of the construction of the miRNA‐mRNA regulatory network in ovarian endometriosis. DEG, differentially expressed gene; DEM, differentially expressed miRNA; EC, ectopic endometrium; EU; eutopic endometrium; OEM, ovarian endometriosis.

None.

In all, 13 tissue samples, including paired EC and EU tissues from five patients with OEM and normal endometrial tissues of three volunteers, were obtained from the Beijing Obstetrics and Gynecology Hospital. All of the patients were confirmed to have a moderate to severe (III–IV) stage of OEM by histological examination according to the revised American Fertility Society (r‐AFS) classification during laparoscopic surgery. 22 None of the women was taking any types of hormonal suppression and they were confirmed to have regular menstrual cycles using the method of Noyes et al. 23 Clinical characteristics of the participants are listed in Table  S1 . All tissue samples were divided into two parts: one half was fixed and prepared for pathological examination and the other half was snap‐frozen in liquid nitrogen and stored at −80°C until further use. Tissues were ground under liquid nitrogen in a mortar and pestle, and the resulting powder solubilized in 1 ml of TRIzol™ Reagent (Invitrogen, Carlsbad, CA, USA) in a FastPrep microfuge tube containing Lysing Matrix D ceramic beads. The RNA was subsequently isolated from the resulting TRIzol solution and further purified with an additional RNeasy step (QIAGEN, Redwood City, CA, USA). The quality and purity of RNA were examined by a NanoDrop™ One/OneC spectrophotometer (Thermo Scientific, Waltham, MA, USA) and Life Invitrogen Qubit® 3.0 Fluorometer via using Qubit™ RNA HS Assay Kit. RNA integrity was analyzed using Agilent 4200 TapeStation system (Agilent, Santa Clara, CA, USA). Samples with an absorbance wavelength ratio (A260/A280) ≥1.9 and an RNA integrity number ≥8 were included. For mRNA sequencing, 3 μg of RNA per sample was prepared for library preparation. The Ribo‐Zero rRNA Removal Reagent (H/M/R) (Illumina, San Diego, CA, USA), magmatic RNA Beads (Vazyme, Nanjing, China) and NEBNext UItra Directional RNA Library Prep Kit for Illumina (NEB, Ipswich, MA, USA) were used for rRNA removal and library construction following the manufacturers’ protocols. For high‐throughput sequencing, paired‐end 150‐bp sequencing of the cDNAs was performed on an Illumina PE150 system conducted by HaploX Genomics Center (Shangrao, China). The raw data were processed with fastp (v0.20.0) to ensure the quality of the data used in subsequent analyses. 24 HISAT2 (Hierarchical Indexing for Spliced Alignment of Transcripts, v2.1.0) was used to build the genome index, and mapped the clean data to the human genome (GRCh37) ( http://grch37.ensembl.org/index.html ). 25 HTSeq (v0.10.0) was used to quantify the read counts of each gene, 26 and the fragments per kilobase of exon per million fragments mapped (FPKM) were calculated to estimate the expression level of genes in each sample. DEGSEQ2 (v1.18.1) 27 was used to analyze differentially expressed genes (DEGs) with parameters of FDR  1. For small RNA sequencing, 3 μg of total RNA was used for library preparation. The purified libraries were concentrated by KAPA Library Quant kit (Illumina) and sequenced on an Illumina NovaSeq 6000 platform with matched NovaSeq S4 reagent kit by HaploX Genomics Center. The clean data were generated from raw data and mapped to the Ensemble database (GRCh37) ( http://grch37.ensembl.org/index.html ) using FASTP 24 and the Burrows–Wheeler Alignment tool (BWA). 28 The data were then compared with MIRBASE (Release 22) 29 to identify mature miRNAs and processed by MIREVO 30 integrated software platform to identify novel miRNAs. For each sample, the count and TPM (transcripts per million reads) 31 value of the miRNAs were collected. DEGSEQ (v1.18.0) 27 was used to analyze differentially expressed miRNAs (DEMs) with parameters of FDR  1. We used “endometriosis” as a keyword on the GEO database: one mRNA‐sequencing dataset ( GSE105764 21 ) and one miRNA‐sequencing dataset ( GSE105765 21 ) were collected. Eight paired EC and EU samples were included in GSE105764 and GSE105765 . The platform and series matrix files were downloaded. The sra files of the GSE105764 and GSE105765 original sequencing data were transformed into FASTQ files by SRATOOLKIT.2.11.0‐CENTOS (version 2.11). RSTUDIO software (version 3.6) was used to process and standardize the sequencing data of local and GEO datasets. SVA package (version 3.40.0) and bladderbatch package (version 1.30.0) were used to remove batch effect. DEGs and DEMs were identified by using the limma package (version 3.40.6) with a p‐ value  1. The miRNA‐mRNA interactions were predicted using the TargetScan database (v7.1) with a z ‐score ≥0.7. 32 The Pearson correlation coefficient (PCC) was used to measure the relationships with a PCC >0.6 and a p ‐value <0.05. The intersecting mRNAs were selected and the integrative miRNA‐mRNA regulatory network was visualized by CYTOSCAPE (v3.2). 33 To explore the functional roles of DEGs in the miRNA‐mRNA networks, we used DAVID, which integrates the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases to analyze biological functions. 34 The enrichment values of GO terms and KEGG pathways were implemented by the hypergeometric test, and the q ‐value (adjusted as p ‐value) <0.05 was considered to be significantly enriched. The study protocol was approved by the Local Ethical Committee of the Beijing Obstetrics and Gynecology Hospital, Capital Medical University (Approval number: 2021‐KY‐022‐01) on March 15, 2021. All participants provided written informed consent.

Figure S1 Click here for additional data file. Table S1 Click here for additional data file. Table S2 Click here for additional data file. Table S3 Click here for additional data file. Table S4 Click here for additional data file. Table S5 Click here for additional data file. Table S6 Click here for additional data file. Table S7 Click here for additional data file. Table S8 Click here for additional data file.