This study included patients who were diagnosed with endometriosis through laparoscopy at Peking University First Hospital. Exclusion criteria included individuals experiencing abnormal menstrual cycles, those who had taken hormonal medications or used intrauterine devices in the three months preceding the surgery, and individuals with medical complications or pelvic inflammatory disease. A group of individuals ranging in age from 27 to 40 years who were experiencing menstrual periods in the proliferative phase were chosen for inclusion in the disease group. The control group consisted of ten patients (aged 28–40 years) who had laparoscopy for ovarian teratoma, simple ovarian cyst, fallopian tube cyst, and uterine fibroids during the proliferative phase of menstruation. Supplementary Table 1 displays the characteristics of the patients. The research was carried out following the guidelines of the Declaration of Helsinki, and the Ethics Committee of Peking University First Hospital approved the protocol (approval number No. 2020KEYAN343). Following the acquisition of informed consent from every participant, 10 mL of plasma was extracted from both the patient and control groups.

As per the established procedure, plasma samples were spun at 2000 $\times$g for 10 minutes at 4 °C in order to eliminate deceased cells and cellular waste. The resulting liquid was moved into a fresh tube for another round of centrifugation at a force of 10,000 times the acceleration due to gravity for a duration of 30 minutes to eliminate larger vesicles and waste materials. Afterwards, the specimen underwent ultimate ultracentrifugation at a speed of 110,000 $\times$g for a duration of 75 minutes at a temperature of 4 °C. After being washed with phosphate buffered saline (PBS, catalog# P1020, Solarbio, Beijing, China), the extracellular vesicle (EV) pellets were passed through a 0.22 µm filter. Another round of ultracentrifugation (UC) was conducted at a speed of 110,000 $\times$g for an additional duration of 75 minutes at a temperature of 4 °C. The fraction enriched with electric vehicles was suspended in a volume ranging from 100 to 200 µL and then frozen at –80 °C for storage.

Initially, the exosome microstructure was examined using transmission electron microscopy (TEM). A copper grid was used to place a 10 µL sample of extracellular vesicles, which was then precipitated for 1 min. Afterwards, the excess liquid was eliminated using filter paper. Afterwards, phosphotungstic acid was introduced onto the copper grid and left to precipitate for a duration of 1 minute, following which the excess liquid was once more eliminated using filter paper. The specimen was permitted to air dry for a few minutes at ambient temperature. The specimen was left to dry naturally at ambient temperature for a few minutes. After the sample had dried, imaging results were obtained using a Hitachi HT7700 electron microscope (Hitachi, Ltd., Tokyo, Japan) with an acceleration voltage set between 40–120 kV (in 100 V increments). Ultimately, the instrument yielded imaging results through transmission electron microscopy.

Furthermore, the exosomes’ size distribution and quantification were examined using nanoparticle tracking analysis (NTA). The exosome samples, which were isolated, were diluted in 1$\times$ PBS buffer and then loaded into the cell. At 11 different positions, the device autonomously recorded and examined every sample. ZetaView Nanoparticle Tracking Analyzer (Particle Metrix, Meerbusch, North Rhine-Westphalia, Germany) was utilized to determine the concentration and diameter dimensions of the samples.

In the end, an analysis called western blotting was performed to investigate the presence of exosome-positive indicators (CD9, CD63, and TSG101) as well as an exosome-negative indicator (calnexin). Exosome lysates were acquired by combining radio immunoprecipitation assay (RIPA, Catalog# R0010, Solarbio, Beijing, China) solution with protease and phosphatase inhibitors, then measured using a bicinchoninic acid assay (BCA, Catalog# PC0020Solarbio, Beijing, China) Protein Assay Kit. Equivalent quantities of exosomal proteins from every group were analyzed using sodium dodecyl sulfate polyacrylamide gel electrophoresis (SDS‒PAGE, Catalog# P1200, Solarbio, Beijing, China), followed by the transfer of proteins on the gel onto a nitrocellulose membrane. The filters were obstructed using 5% low-fat milk and subsequently placed in an incubator with the primary antibodies for an extended period at a temperature of 4 °C. Following the wash with phosphate buffered solution (PBST, Catalog# P1033, Solarbio, Beijing, China), the membranes were subjected to incubation with a secondary antibody that had been modified with horseradish peroxidase.

Afterwards, the exosome size distribution and quantification were examined using nanoparticle tracking analysis (NTA). After being moderately diluted with 1$\times$ PBS buffer, the isolated exosome samples were introduced into cells. Samples from 11 distinct sites were automatically recorded and analyzed by the device. Using a ZetaView nanoparticle tracking analyzer (Particle Metrix, Meerbusch, North Rhine-Westphalia, Germany), the samples’ concentration and diameter size were determined.

In the end, western blotting detected exosome-positive indicators (CD9, CD63, and TSG101) as well as exosome-negative indicator (calnexin). To acquire exosome lysates, RIPA-containing protease and phosphatase inhibitors were added and then quantified using a BCA protein analysis kit. The same quantity of exosome proteins from both groups underwent SDS‒PAGE, followed by the transfer of the proteins onto a nitrocellulose membrane. The membranes were obstructed using 5% skim milk and subsequently subjected to incubation with primary antibody for an extended period at 4 °C. Following the washing process using PBST solution, the membranes were subjected to incubation with a secondary antibody conditioned with horseradish peroxidase.

The main antibodies utilized were rabbit anti-CD63 [1:1000; System Biosciences (SBI), Catalog# EXOAB-CD63A-1; Palo Alto, CA, USA], rabbit anti-CD9 (1:1000; SBI, Catalog# EXOAB-CD9-1; Palo Alto, CA, USA), rabbit anti-TSG101 (1:1000; SBI, Catalog# EXOAB-TSG101-1; Palo Alto, CA, USA), and rabbit anti-Calnexin (1:1000; Immunoway Bio-technology , Catalog# YM0089, Plano, TX, USA). Exosome-validated goat anti-rabbit IgG (CD63, CD9; 1:5000; SBI) and BeriBlot for IP Detection Reagent (Calnexin, 1:2000; Catalog# ab10286, Abcam, Cambridge, UK) were the secondary antibodies used.

The miRNeasy Micro Kit (QIAGEN, Catalog #217084; Hilden, North Rhine-Westphalia, Germany) was utilized for the isolation of total RNA. Next, the RNA samples were utilized to produce a biotinylated cRNA that targeted the gene chip® miRNA 4.0 array. In short, after that, marked examples were added to the mixture for hybridization and then incubated as per the instructions provided by the manufacturer before being injected into the miRNA array. Following a 16-hours hybridization at a temperature of 48 °C, the array underwent washing, staining, and scanning with a GeneChip scanner (3000 7G GeneChip Scanner, Affymetrix, Santa Clara, CA, USA). Wayen Biotechnologies Inc. (Shanghai, China) carried out the chip experiment in accordance with Thermo Inc. protocol.

First, the quality was assessed by scanning the microarrays and plotting box line plots. After annotating and counting the microarray data, the annotation information and expression values of miRNAs were derived. To compare and screen differentially expressed miRNAs between EM and control plasma exosomes, the R package ‘limma’ (http://www.bioconductor.org/packages/release/bioc/html/limma.html) [15] was utilized. The occurrence of false positive results was corrected by applying the Benjamini and Hochberg false discovery rate (FDR) method and adjusting the p value. To identify differential miRNAs (DE-miRNAs), the thresholds were set as an adjusted p value 0.05 and $|$log2-change (FC)$|$ $>$0.5. The DE-miRNAs identified by ‘limma’ were visualized using volcano plots and heatmaps.

For this research, we developed a weighted gene co-expression network analysis (WGCNA) by utilizing miRNA expression data to create a network with weighted coexpression. To begin with, a total of 20 samples were subjected to hierarchical clustering using the Euclidean distance of expression. The purpose was to identify any outliers within the samples and subsequently remove them [16]. Next, the similarity matrix was constructed by calculating the correlation degree of each miRNA pair. To establish the correlation threshold between miRNAs, a suitable power within the range of 1 to 20 was chosen as the soft threshold. The relationships between the scale-free network evaluation coefficient R2 and the soft threshold $\beta$, as well as between the average connectivity and the power of the soft threshold, were established. The weighted parameter of the neighborhood function is the soft threshold $\beta$, and the optimal value is obtained using the pick soft threshold function in the R package ‘WGCNA’. Next, the adjacency matrix underwent a conversion to become a topological overlap matrix (TOM). TOM has the ability to indicate the network connection of a gene, which refers to the total of its connections to all other genes generated by the network [17]. Hence, the identification of modules is achieved through hierarchical clustering of genes using the dissimilarity measure based on TOM. In the end, by utilizing the dissimilarity matrix diss TOM = 1 ‒ TOM, we formed a hierarchical clustering tree to group the miRNAs. The mRNAs exhibiting comparable expression patterns were categorized into the identical module. The higher the absolute value of module significance (MS) and gene significance (GS), the more relevant the miRNA module was to EM. Those modules that are highly correlated with clinical features can be identified as candidate modules [18].

The intersection of these DE-miRNAs and module miRNAs yielded the key DE-miRNAs. We examined the significant DE-miRNAs to determine their enrichment in signaling pathways using the miRNA Enrichment Analysis and Annotation tool (miEAA) [19]. We searched two databases: Gene Ontology (GO) Pathway and Kyoto Encyclopedia of Genes and Genomes (KEGG) Pathway. Retained were ontology terms and KEGG pathways that had a significant adjusted p value (0.05).

Evaluating the quality or performance of diagnostic tests was effectively done using the receiver operating characteristic (ROC) curve, which plotted the y coordinate as test sensitivity and the x coordinate as 1-specificity or false positive rate (FPR). To obtain the ROC curve and the area under the curve (AUC) value, the ‘pROC’ package was utilized [20]. The closer the AUC was to 1, the more accurate the diagnosis. The diagnostic effect was low when the AUC was between 0.5 and 0.7; there were some reference values in the diagnostic effect when the AUC was between 0.7 and 0.9; the diagnostic accuracy was high when the AUC was greater than 0.9; and there was no diagnostic value when the AUC was less than or equal to 0.5.

The mirDIP, TargetScan, and MiRTarBase databases can be accessed through the following links: mirDIP: http://ophid.utoronto.ca/mirDIP/; TargetScan: https://www.targetscan.org/vert_72/; MiRTarBase: https://mirtarbase.cuhk.edu.cn/~miRTarBase/miRTarBase_2022/php/index.php. The prediction of target mRNAs for diagnostic miRNAs was done using miRNA‒mRNA pairs screened from these three databases were then crossed. The Venn Diagram R package was utilized to create a Venn Diagram illustrating the intersection. To predict the diagnostic miRNAs’ target long noncoding RNAs (lncRNAs), the miRNet database (https://www.mirnet.ca/miRNet/upload/MirUploadView.xhtml), starBase database (https://rnasysu.com/encori/), and NPInter database (http://bigdata.ibp.ac.cn/npinter4/) were utilized. After that, the miRNA‒lncRNA relationship screened from these three databases was then crossed. The Venn Diagram R package was utilized to create a Venn Diagram illustrating the intersection.

A competitive endogenous RNA (ceRNA) network was constructed by utilizing the interplay among lncRNAs, miRNAs, and mRNAs. Next, the data were visualized using Cytoscape software (The Cytoscape Consortium, San Diego, CA, USA), and the analyze network tool was used to calculate the nodes’ degree. LncRNAs, miRNAs, and mRNAs with degrees in the ceRNA network were considered hub lncRNAs, miRNAs, and mRNAs.

First, we isolated plasma exosomes from 10 control and 10 EM patients. Exosomes were separated from the serum using ultracentrifugation. We employed TEM, NTA, and western blotting techniques to ascertain the successful isolation of the exosomes. Clear vesicle structures were visible under TEM (Fig. 1A). According to the NTA findings, the exosomes fulfilled the requirements in terms of their concentration (ranging from 3.5 $\times$ 10${}^7$ to 8.7 $\times$ 10${}^7$ particles/mL) and size (measuring between 112 and 130 nm), as shown in Fig. 1B.

Fig. 1.

Detecting exosomes derived from the endometrium. (A) Transmission electron microscopy (TEM) was used to observe the morphologies of exosomes derived from plasma. Scale bar = 100 nm; (B) Nanoparticle tracking analysis (NTA) was used to determine the dimensions and density of exosomes; (C) Exosomes were analyzed by western blotting using antibodies against exosomal markers (Calnexin, CD63, CD9 and TSG101).

We examined the presence of exosome-associated indicators (CD9, CD63, and TSG101) as well as a marker indicating the absence of exosomes (calnexin). The findings indicated that CD9, CD63, and TSG101 were expressed positively in the exosomes, while calnexin showed negative expression. On the contrary, the serum exhibited the opposite expression when exosomes were eliminated. Thus, we obtained qualified exosomes (Fig. 1C).

miRNA microarrays were used to extract and sequence total RNA from the eligible exosomes mentioned above. After quality control and preprocessing, we acquired the annotation information and expression values of the miRNAs. We utilized the R package ‘limma’ to conduct differential expression analysis to identify the distinctively expressed exosomal miRNAs among EM patients and control subjects. The DE-miRNAs were visualized by volcano plots (Fig. 2A). We obtained a total of 50 DE-miRNAs, consisting of 7 miRNAs that were upregulated and 43 miRNAs that were downregulated. The expression levels of the DE-miRNAs are shown in Fig. 2B.

Fig. 2.

Endometriosis is associated with differential expression of microRNAs (miRNAs) in exosomes. (A) Volcano maps of differentially expressed exosomal miRNAs. (B) Heatmaps displaying the miRNAs based on the magnitude of $|$logFC$|$ in the endometriosis (EM) and control groups.

WGCNA was conducted on all miRNAs in this research utilizing the WGCNA function in the R package. To begin with first, a total of 20 samples were subjected to hierarchical clustering using the Euclidean distance of expression. The purpose was to identify any outliers within the samples and subsequently remove them [16]. As depicted in Fig. 3A, there were no apparent samples that stood out.

Fig. 3.

Analysis of dysregulated miRNAs using a weighted gene coexpression network. (A) Sample hierarchical clustering by Euclidean distance of expressions. (B) Examining the scale-free fit index and average connectivity, we analyze the scale-free topology for different soft-threshold powers ($\beta$). (C) The dendrogram displays the clustering of differentially expressed miRNAs in the two comparison sets using a dissimilarity measure 1 ‒ topological overlap matrix (1 ‒ TOM) and includes assigned module colors. (D) Analyzing the connections between modules and traits in relation to endometriosis. A module eigengene is represented by each row, while each trait is represented by each column. Control and EM are indicated by groups. (E) A scatterplot displaying the correlation between gene significance and module membership in the blue coexpression modules.

In order to achieve a network that fulfills the scale-free topology requirement, we computed various network structures by employing diverse soft-thresholding powers spanning from 1 to 20. Distinct gene coexpression modules in EM were determined by the power value ($\beta$) and R2 value. The scale-free topological fitting index of the mRNA coexpression network reached 0.85 when the soft-thresholding power was 5 (Fig. 3B), which met the scale-free network criterion. Several modules were formed and colored differently in the scale-free networks, which were built using the screened soft thresholds for the miRNAs. Next, the mRNAs were grouped into multiple modules using hierarchical clustering based on their expression values, resulting in the identification of 27 modules (Fig. 3C).

Phenotypes and modules were correlated by utilizing samples from the EM and control groups as phenotypes. The matrix of correlation coefficients between module eigenvectors and phenotypic traits was computed, and heatmaps illustrating the correlations were generated [16]. According to the module-EM relationships in Fig. 3D, we determined that the blue module revealed the strongest correlation (module-trait weighted correlation = –0.47, p 0.05) with EM and was identified as a key module for EM, containing 320 miRNAs.

In order to further determine the key miRNAs, the internal module was examined to identify the genes exhibiting high GS (trait correlation with module gene, gene significance) and MM (gene correlation with module, module membership) [16]. The significant correlation between module membership (MM) in the blue module and gene significance (GS) for EM is presented in Fig. 3E. The coexpression modules that have been identified are displayed in various hues, with the blue module revealing a total of 38 core miRNAs.

The 50 differential miRNAs (DE-miRNAs) were intersected with 38 blue module core miRNAs screened by WGCNA. Thirty-one overlapping genes from DE-miRNAs and module miRNAs were retained as key DE-miRNAs (Fig. 4A).

Fig. 4.

Functional enrichment analysis of differentially expressed mRNAs. (A) A comparison of differentially expressed genes (DEGs) in the three sets using a Venn diagram. (B) Histogram of the top 10 Gene Ontology (GO) annotation terms of key differential miRNAs (DE-miRNAs). (C) Plotting a circular diagram illustrating the top 10 GO annotation terms associated with important DE-miRNAs. (D) Bubble plot of the top 10 Kyoto Encyclopedia of Genes and Genomes (KEGG) annotation terms of key DE-miRNAs. (E) Chord plot of the top 10 KEGG annotation terms of key DE-miRNAs in the enrichment analysis of dysregulated mRNAs between the EM and control groups. TNF/Th1, the tumor necrosis factor/T helper 1; CCR4-NOT, carbon catabolite repressor 4-Negative on TATA box (Goldberg–Hogness box).

To explore the biological functions and signaling pathways in which key DE-miRNAs were implicated, GO and KEGG enrichment analyses were conducted. A total of 3393 biological functions had 31 DE-miRNAs enriched when using p 0.05 as the screening criterion. The most important GO annotation terms, which are more significant with smaller p values, were represented visually. Histograms were used for this visualization (Fig. 4B), and circle plots (Fig. 4C) were drawn. The enrichment analysis revealed that the primary DE-miRNAs were predominantly enriched in ontology terms such as the ‘pathway of insulin receptor signaling’, ‘development of outflow tract’, ‘binding of translation initiation factor’, ‘binding of insulin-like growth factor receptor’, and ‘promotion of mesenchymal cell proliferation’.

A total of 223 signaling pathways had 31 DE-miRNAs enriched when using p 0.05 as the screening criterion. The most important 10 KEGG annotation terms, with smaller p values indicating higher significance, were represented through bubble plots (Fig. 4D) and chord plots (Fig. 4E). The main enrichment of the key DE-miRNAs occurred in toxoplasmosis, the tumor necrosis factor (TNF) signaling pathway, the phospholipase D signaling pathway, Chagas disease American trypanosomiasis, central carbon metabolism in cancer, the Toll-like receptor signaling pathway, and T helper 1 (Th1) and Th2 cell differentiation.

The Proc data package in R was utilized to plot the ROC curves of important DE-miRNAs in this research, and the validation of each crucial DE-miRNA’s diagnostic significance for endometriosis was conducted. AUC $>$0.7 was used as a screening index for ROC for the diagnosis of endometriosis. The analysis of the ROC revealed that out of the 31 important DE-miRNAs, one indicator had an AUC value below 0.7 (we removed this indicator). The remaining 30 DE-miRNAs had an AUC value above 0.7, thus making them suitable as diagnostic marker miRNAs (Fig. 5).

Fig. 5.

Diagnostic value of serum exosomal miRNAs. Receiver operating characteristic (ROC) curves were produced using each miRNA expression value for key DE-miRNAs. AUC, area under the curve.

We utilized three databases, namely, mirDIP, TargetScan, and MiRTarBase, to screen miRNA‒mRNA interaction pairs and forecast the target mRNAs for 30 diagnostic miRNAs. Fig. 6A shows that there were 362 pairs of miRNA-mRNA relationships, with 337 being mRNAs and 24 being miRNAs. To predict 30 diagnostic miRNAs for targeting lncRNAs, we screened miRNA‒lncRNA interaction pairs using three databases: miRNet, starBase and NPInter. Fig. 6B shows that there were 167 pairs of interactions between miRNA and lncRNA, comprising 24 different miRNAs and 37 different lncRNAs.

Fig. 6.

The networks of competitive endogenous RNA (ceRNA) regulation involving lncRNA, miRNA, and mRNA. (A) Venn diagram of miRNA‒mRNA intersections obtained from the mirDIP, TargetScan and MiRTarBase databases. (B) Venn diagram of miRNA‒lncRNA intersections obtained from the miRNet, starBase and NPInter databases. (C) mRNA‒miRNA regulatory network; red circles are miRNAs, and yellow triangles are mRNAs. (D) miRNA‒lncRNA regulatory network; red circles are miRNAs, and blue squares are lncRNAs. (E) Core mRNA‒miRNA‒lncRNA regulatory network; red circles are miRNAs, blue squares are lncRNAs, and yellow triangles are mRNAs.

Next, we utilized Cytoscape software to combine lncRNA‒miRNA and miRNA‒mRNA interaction pairs into a network consisting of mRNA‒miRNA‒lncRNA triplets. The regulatory network of lncRNA-miRNA-mRNA consists of 1796 pairs of relationships, encompassing 36 lncRNAs, 20 miRNAs, and 264 mRNAs (Fig. 6C).

The ‘Analyze Network’ tool in Cytoscape software was used to calculate the degree of each node in accordance with the ceRNA regulatory network. We screened the top 10 hub lncRNAs (GAS5, MALAT1, FGD5-AS1, HCG18, SNHG16, XIST, OIP5-AS1, NEAT1, KCNQ1OT1 and SNHG12), hub miRNAs (hsa-miR-361-5p, hsa-miR-19b-3p, hsa-let-7f-5p, hsa-miR-23a-3p, hsa-miR-199a-3p, hsa-miR-18a-5p, hsa-miR-221-3p, hsa-miR-17-5, hsa-miR-27a-3, and hsa-miR-25-3p), and hub mRNAs (GALC, ETNK1, RNF4, SOX4, ZBTB18, SPRY2, RUNX1, MYLIP, BTG2, and MAP2K4). The central network of mRNA-miRNA-lncRNA regulation was built, as depicted in Fig. 6D,E.