Global RNA expression analysis of patient samples identified potential diagnostic biomarkers specific for peritoneal, ovarian and deep endometriosis

Table  1 shows the demographics and clinical characteristics of the participants ( n  = 81) included in RNAseq cohort ( n  = 42) as well as in ELISA validation cohort ( n  = 44) (5 participants were included in both cohorts). No significant differences were observed for age, BMI, menarche or cycle length across the groups. Additional information on participants’ pain perception was collected and rated on a numeric rating scale from 1 to 10. Pelvic pain perception was higher in women with all types of endometriosis, however only women with peritoneal endometriosis (1.3 vs. 3.5, p  = 0.009) and deep endometriosis (1.3 vs. 2.7, p  = 0.03) reported significantly higher pain scores, compared to women without diagnosed endometriosis (Table  1 ). Scores for dysmenorrhea (painful menstruation) and dyspareunia (painful intercourse) were also higher across all types of endometriosis, however reaching significance only for menstrual pain associated with deep endometriosis (4.1 vs. 7.6, p  = 0.0002) and ovarian endometriosis (4.1 vs. 6.3, p  = 0.04). Additionally, higher score for painful defecation and infertility was again only associated with deep endometriosis (dyschezia: 1 vs. 3.4, p  = 0.003; sterility 31% vs. 67%, p  = 0.02). Higher scores for infertility in the control population are related to other gynecological conditions, as the control women without endometriosis were recruited from patients scheduled for laparoscopy or hysteroscopy due to conditions such as fibroids, unexplained pelvic pain and dermoids, or other benign cysts. Table 1 The clinical and demographic characteristics of participants included in the study. Statistical significance between groups was determined using t-test and actual p-values are displayed. Statistically significant differences are highlighted in green. The clinical and demographic characteristics of participants included in the study. Statistical significance between groups was determined using t-test and actual p-values are displayed. Statistically significant differences are highlighted in green. To investigate differentially expressed genes in three subtypes of endometriosis (ovarian endometriosis – OE, peritoneal endometriosis – PE, deep endometriosis – DE) we performed RNA sequencing analysis capturing both mRNAs as well as limited amount of lncRNAs as it was shown that 40–70% of lncRNAs do contain poly-A tail 20 , 21 . We analyzed 11 samples of OE, 5 samples of PE, 11 samples of DE and 15 control normal endometrium samples, each sample obtained from different patient. First, we investigated overall variability of gene expression profiles of different groups of endometrial samples in comparison to control samples. PCA plot shows 3 discrete clusters representing control samples, OE samples and PE samples together with DE samples (Fig. 1 A). This clustering is in line with previously published gene expression profiles and reflect different pathology and emergence routes 5 , 7 . To further evaluate the major differences in expression profiles between evaluated groups, we constructed heat map including all significantly differentially expressed genes across the groups with p-value less than 0.05 (Fig. 1 B). The analysis revealed that samples from control endometrium and from ovarian endometriosis each displayed internally consistent gene expression profiles, while gene expression profiles for individual samples of extraovarian endometriosis (DE and PE) are similar to each other and did not separate into distinct groups, recapitulating the results from PCA (Fig. 1 A). Fig. 1 Principal Component Analysis and Heatmap of patient samples that were analyzed using RNA sequencing. ( A ): PCA plot represents the distribution of samples based on the first two principal components (PC1 and PC2). Each point represents an individual sample, colored according to tissue type (OE, DE, PE, control endometrium). ( B ): A hierarchical clustering analysis was performed on both genes and samples, and a heatmap was generated to depict the expression levels of all significantly differentially expressed genes ( p  < 0.05) across the different sample types (OE, DE, PE, and control endometrium). Each row represents a gene and each column represents a sample, with clustering applied in both dimensions to highlight expression patterns. Expression levels are normalized and displayed as z-scores, with colors indicating relative expression (red = upregulated, blue = downregulated). Dendrograms illustrate hierarchical clustering of genes and samples based on similarity. The color key represents the intensity of expression changes, and annotations indicate sample types. Alt text: Fig. 1 shows the results of RNA sequencing analysis through Principal Component Analysis (PCA) and a heatmap of gene expression. Panel A displays a PCA plot with each point representing a patient sample, color-coded by tissue type: ovarian endometriosis (OE), deep endometriosis (DE), peritoneal endometriosis (PE), or control endometrium. Panel B presents a heatmap of significantly differentially expressed genes ( p  < 0.05), with rows as genes and columns as samples. Principal Component Analysis and Heatmap of patient samples that were analyzed using RNA sequencing. ( A ): PCA plot represents the distribution of samples based on the first two principal components (PC1 and PC2). Each point represents an individual sample, colored according to tissue type (OE, DE, PE, control endometrium). ( B ): A hierarchical clustering analysis was performed on both genes and samples, and a heatmap was generated to depict the expression levels of all significantly differentially expressed genes ( p  < 0.05) across the different sample types (OE, DE, PE, and control endometrium). Each row represents a gene and each column represents a sample, with clustering applied in both dimensions to highlight expression patterns. Expression levels are normalized and displayed as z-scores, with colors indicating relative expression (red = upregulated, blue = downregulated). Dendrograms illustrate hierarchical clustering of genes and samples based on similarity. The color key represents the intensity of expression changes, and annotations indicate sample types. Alt text: Fig. 1 shows the results of RNA sequencing analysis through Principal Component Analysis (PCA) and a heatmap of gene expression. Panel A displays a PCA plot with each point representing a patient sample, color-coded by tissue type: ovarian endometriosis (OE), deep endometriosis (DE), peritoneal endometriosis (PE), or control endometrium. Panel B presents a heatmap of significantly differentially expressed genes ( p  < 0.05), with rows as genes and columns as samples. To further analyze gene expression profiles in more detailed manner, we performed pathway analysis of global gene expression alterations in individual subtypes of endometriosis. Our analysis identified several distinct pathways to be upregulated across all three types of endometriosis such as increased extracellular matrix deposition and signaling, increased contractility characteristic for fibrotic conditions, as well as increase chemokine and cytokine signaling characteristic for inflammatory conditions (Fig. 2 A-C, Supplementary Table 2 ). Additionally, we identified several pathways that are specific for individual type of endometriosis or pathways that can differentiate between ovarian and extraovarian (DE and PE) endometriosis. One of the prominent pathways identified in our sequencing analysis in deep endometriosis as well as in peritoneal endometriosis is cyclic phosphodiesterase signaling (Fig. 2 B and C, Supplementary Table 2). PDE1 class of phosphodiesterases represented by PDE1A, 1B and 1 C, as well as PDE2A are significantly upregulated in all types of endometriosis, while PDE3A and PDE5A are upregulated only in DE (Supplementary Table 3). PDE1-class phosphodiesterases are regulated through calmodulin binding which was also identified as enriched pathway in our dataset (Fig. 2 B and C). Increased PDE1 class activity was not previously shown in endometriosis, however it is shown to be important for pathology of other fibrotic diseases such as idiopathic pulmonary fibrosis or liver fibrosis 22 , 23 . On the other hand, pathways enriched in ovarian endometriosis are significantly different from the ones enriched in DE and PE. Major enriched pathway is MHC class II receptor pathway represented by several classes of HLA (HLA-DQ or HLA-DR) 24 , 25 . MHC class II receptors are mainly expressed on antigen presenting cells however can be also expressed by cells of epithelial origin 26 , 27 . Moreover, previous IHC-based as well as single cell RNA seq works identified increased expression of MHC class II receptors in ovarian endometriosis 24 , 25 . Additionally, we analyzed the expression data using xCell module 28 to investigate relative abundance of various cell types within the endometriotic lesions as well as in control endometrium. Our in silico analysis confirmed that control endometrial tissue was mostly composed of cells with epithelial origin and contained limited amount of stromal cells (Supplementary Fig. 1, Supplementary Table 4). While on the other hand, endometriotic tissues (OE, DE, PE) were strongly enriched in cells with mesenchymal origin (e.g. fibroblasts) which is in agreement with the pathophysiology of the endometriosis 29 , 30 as well as previous single cell RNA sequencing results 31 – 33 . Fig. 2 Pathway enrichment analysis. Bubble enrichment plot representing the top 20 curated significantly enriched Gene Ontology pathways (left: GO molecular function, right: GO cellular component) identified in each endometriosis subtype ( A ): Ovarian Endometriosis; ( B ): Deep endometriosis; ( C ): Peritoneal endometriosis). The x-axis represents the enrichment score, while the y-axis lists the enriched pathways. Color intensity indicates significance (FDR), and dot size represents gene count. Alt text: Bubble enrichment plots illustrating the top 20 significantly enriched Gene Ontology (GO) pathways for each endometriosis subtype. Pathway enrichment analysis. Bubble enrichment plot representing the top 20 curated significantly enriched Gene Ontology pathways (left: GO molecular function, right: GO cellular component) identified in each endometriosis subtype ( A ): Ovarian Endometriosis; ( B ): Deep endometriosis; ( C ): Peritoneal endometriosis). The x-axis represents the enrichment score, while the y-axis lists the enriched pathways. Color intensity indicates significance (FDR), and dot size represents gene count. Alt text: Bubble enrichment plots illustrating the top 20 significantly enriched Gene Ontology (GO) pathways for each endometriosis subtype. To identify the most differentially expressed genes in each subgroup that could be used as diagnostic biomarkers, we constructed volcano plots (Fig. 3 A-D). Interestingly, the analysis shows that the fold change of gene expression for particular genes was greatly variable with some of the genes showing up to 14-fold increase in certain conditions indicating major differences in signaling pathways governing the pathophysiology of endometriosis. (Fig. 3 A-D, Supplementary Table 1). First we investigated what are the most upregulated genes common for all types of endometriosis. In line with previously published results and pathophysiology of endometriosis, genes involved in fibrotic conditions and myofibroblast transition such (COL10A1, MYH1, MYH11) 34 – 36 and inflammation (CCL11, C7, PLA2G2A, PLA2G5) 37 – 39 . Interestingly, components of non-canonical Wnt signaling pathways are among the most overexpressed genes (SFRP2, TBX18, TCF21 or PDLIM3) indicating a strong, yet underexplored influence of this pathway on pathophysiology of endometriosis (Fig. 3 A-D; Table 2 ). Additionally, several of the most upregulated gene products are shown to be secretory which could be used for development of non-invasive blood-based test for endometriosis pending further validation – SFRP2, CCL11, LEP, C7, FIBIN, PLA2G5 40 – 45 (Fig. 3 A; Table 2 ). Another challenge in management of endometriosis is diagnosis of individual types of endometriosis, therefore we set out to identify genes, that are specifically upregulated in individual types (OE, DE, PE). Interestingly, a major transcription factor that is upregulated to the largest extend in DE is NK2 homeobox 3 (NKX2-3) which has not been implicated in endometriosis yet (Fig. 3 B; Table 2 ). However, given its role in differentiation and fibrotic signaling it could represent a viable target for treatment development 46 , 47 . Additionally, we see several secreted proteins (ANGPTL7, GDF6, ADIPOQ, BMP5) that are among the most upregulated in DE and PE (Fig. 3 B-C Table 2 ). On the contrary, the most upregulated genes associated with ovarian endometriosis are significantly different from DE and PE and are part of the steroid hormone synthesis – NR5A1, CYP11A1, CYP17A1, CYP19A1, DLK1 or STAR (Fig. 3 D; Table 2 ). Interestingly, one of the top upregulated genes was KLHDC8A, a multifunctional protein implicated in ciliogenesis which was not yet associated with endometriosis 48 . Fig. 3 Volcano plots of differentially expressed genes. Volcano plot depicting the distribution of differentially expressed genes between aggregated endometriosis samples and healthy control ( A ), individual endometriosis samples and healthy control ( B , C , D ) and between ovarian and extraovarian (DE/PE) endometriosis samples ( E ). The x-axis represents the log₂ fold change (log₂FC), while the y-axis represents the -log₁₀(p-value). Each point corresponds to a gene, with significantly upregulated genes shown in red, downregulated genes in blue, and non-significant genes in gray. The dashed horizontal line indicates the significance threshold (p-value  1 or < -1). Labeled points indicate key genes of interest. Alt text: Volcano plots displaying differentially expressed genes in various comparisons related to endometriosis. Volcano plots of differentially expressed genes. Volcano plot depicting the distribution of differentially expressed genes between aggregated endometriosis samples and healthy control ( A ), individual endometriosis samples and healthy control ( B , C , D ) and between ovarian and extraovarian (DE/PE) endometriosis samples ( E ). The x-axis represents the log₂ fold change (log₂FC), while the y-axis represents the -log₁₀(p-value). Each point corresponds to a gene, with significantly upregulated genes shown in red, downregulated genes in blue, and non-significant genes in gray. The dashed horizontal line indicates the significance threshold (p-value  1 or < -1). Labeled points indicate key genes of interest. Alt text: Volcano plots displaying differentially expressed genes in various comparisons related to endometriosis. Table 2 mRNA expression levels of differentially expressed genes across sample Groups. Table presents the mRNA expression levels of selected differentially expressed genes between control endometrium and individual types of endometriosis (OE, DE, PE). Gene symbols and description, fold change (log₂FC), secretory category and novelty are included. Positive log₂FC values indicate upregulation, while negative values indicate downregulation in displayed endometriosis subtypes compared to normal endometrium. Red color scale indicate the degree of expresion change.. mRNA expression levels of differentially expressed genes across sample Groups. Table presents the mRNA expression levels of selected differentially expressed genes between control endometrium and individual types of endometriosis (OE, DE, PE). Gene symbols and description, fold change (log₂FC), secretory category and novelty are included. Positive log₂FC values indicate upregulation, while negative values indicate downregulation in displayed endometriosis subtypes compared to normal endometrium. Red color scale indicate the degree of expresion change.. Additionally, we evaluated the expression of non-coding RNAs, such miRNA and lncRNAs as an emerging biomarker class for various conditions, including endometriosis 49 – 51 . Our sequencing efforts identified more than 1600 non-coding RNAs with significant expression (Base Mean > 2 5 ), of which approximately one third not yet described (Supplementary Table 1 ). Filtering the list based on the degree of overexpression (> 2 log2FC) yielded 27 non-coding RNAs that could be potentially used as biomarkers (Table 3 ). Most upregulated non-coding RNAs across all types of endometriosis was MIR202HG and LINC02381 both of which are associated with endometriosis 52 , 53 . Additionally, we identified several non-coding RNAs with increased expression in specific endometriosis type. FENDRR and PGM5-AS1 were highly upregulated in deep infiltration endometriosis (log2FC 4.67, and 4.5 respectively) (Table 3 ). These non-coding RNAs were not yet implicated in endometriosis, however FENDRR was shown to increase IL-6 expression and liver fibrosis, indicating a potential role in endometriosis as well 54 . On the other hand, there are 4 non-coding RNAs with specific and highly significant increase in ovarian endometriosis, namely AFAP1-AS1, HMGA2-AS1, CDH3-AS1 and LINC00607 out of which only AFAP1-AS1 was previously associated with endometriosis promotion through epithelial-to-mesenchymal transition 55 . Since non-coding RNAs are often released in exosomes and are relatively easy to detect, these highly upregulated RNAs could potentially be explored as endometriosis biomarkers. Table 3 LncRNA expression levels across sample Groups. Table presents expression levels of selected differentially expressed the LncRNA genes between control endometrium and individual types of endometriosis (OE, DE, PE). Gene symbols and description, fold change (log₂FC), and novelty are included. Positive log₂FC values indicate upregulation, while negative values indicate downregulation in displayed endometriosis subtypes compared to normal endometrium. Red color scale indicate the degree of expresion change. LncRNA expression levels across sample Groups. Table presents expression levels of selected differentially expressed the LncRNA genes between control endometrium and individual types of endometriosis (OE, DE, PE). Gene symbols and description, fold change (log₂FC), and novelty are included. Positive log₂FC values indicate upregulation, while negative values indicate downregulation in displayed endometriosis subtypes compared to normal endometrium. Red color scale indicate the degree of expresion change. It is well established that immune system plays a definite role in pathology of endometriosis, especially in DE and PE subtypes, however the complete understanding of the interaction between immune system and endometrial and stromal cell is lacking. To infer relative infiltration of different immune cell type within endometriosis lesions, we employed computational tool CIBERSORT which utilize gene expression data and machine learning to estimate the relative abundance of different immune cell types in mixed cell populations 56 , 57 . Our in silico analysis identified 3 major immune cell types that are differentially represented in endometriosis compared to healthy controls – macrophages, NK-cells and B-cells (Fig. 4 A, B). Most significantly enriched immune cell population in all three types of endometriosis were M2 macrophages that are tightly linked to fibrosis (Fig. 4 B, top right panel). There was also a significant decrease of M0 macrophages in DE and PE ( Fig. 4 B, top left panel). Similarly, B-naive cells were highly enriched in all three types of endometriosis, while plasma cells showed great increase of abundance only in OE which corresponds to increased expression of immune interaction genes seen in mRNA expression analysis (Fig. 4 B, bottom panel, Fig. 2 A). On the other hand, there was a significant decrease in NK-cell infiltration as seen on the proportion of both resting as well as activated NK-cells (Fig. 4 A and C, middle panel). Fig. 4 Immune cell infiltration analysis. ( A ): Bar plot showing the estimated proportions of immune cell types in different types of endometriosis and control tissue using CIBERSORT analysis. Each bar represents an individual sample, with different colors indicating the relative abundance of specific immune cell types. B : Box plots displaying the relative abundance of selected immune cell types (Top panel: M0 and M2 macrophages, Middle panel: NK resting and NK activated cells, Bottom panel: Naïve B cells and Plasma cells) across each endometriosis subtype and control endometrium. The y-axis represents the estimated immune cell proportion, while the x-axis represents the different groups. Statistical significance between groups was determined using Wilcox test and actual p-values are displayed. Alt text: Bar and box plots illustrating immune cell composition in endometriosis and control tissues based on CIBERSORT analysis. Immune cell infiltration analysis. ( A ): Bar plot showing the estimated proportions of immune cell types in different types of endometriosis and control tissue using CIBERSORT analysis. Each bar represents an individual sample, with different colors indicating the relative abundance of specific immune cell types. B : Box plots displaying the relative abundance of selected immune cell types (Top panel: M0 and M2 macrophages, Middle panel: NK resting and NK activated cells, Bottom panel: Naïve B cells and Plasma cells) across each endometriosis subtype and control endometrium. The y-axis represents the estimated immune cell proportion, while the x-axis represents the different groups. Statistical significance between groups was determined using Wilcox test and actual p-values are displayed. Alt text: Bar and box plots illustrating immune cell composition in endometriosis and control tissues based on CIBERSORT analysis. Our analysis identified several significantly upregulated genes in individual subtypes of endometriosis that were previously described such as SFRP2, FABP4, CCL11 or C7 58 – 61 further confirming the validity of our results. Moreover, we identified novel genes that were significantly overexpressed in endometriotic subtypes that could be potentially used for diagnostic purposes (Fig. 3 ; Table 2 ). From the list of upregulated proteins, we selected three secreted proteins for targeted validation using ELISA – PLA2G2A as a positive control that was previously shown to be upregulated in endometriotic tissues patients with endometriosis 62 and two novel potential biomarkers – PLA2G5 and ANGPTL7. We processed endometriotic tissues from 40 individuals (9 OE, 9 PE, 15 DE, 7 Control) and evaluated the protein expression of the shortlisted candidates. In line with previously published studies, our analysis showed that the expression of PLA2G2A was significantly increased in all three types of endometriosis, with the highest mean increase in samples from ovarian endometriosis (OE: 399 pg/ml, ctrl: 44pg/ml; p = 0.0015) (Fig. 5 A). On the other hand, protein expression of ANGPTL7 was significantly upregulated in DE (DE: 41ng/ml, ctrl: 22.4ng/ml; p = 0.0037) and PE (PE: 41.2ng/ml, ctrl = 22.4pg/ml; p = 0.00054) while we didn’t see any increase in tissues form ovarian endometriosis which indicates that ANGPTL7 could be used to diagnose specifically DE or PE ( Fig. 5 B). Finally, we analyzed the levels of PLA2G5 across different endometriosis types and found out that PLA2G5 protein levels are significantly increased only in DE subtype (DE: 0.55ng/ml, ctrl = 0.22ng/ml; p = 0.04) which suggest that PLA2G5 could be a specific biomarker for DE however further evaluation would be needed (Fig. 5 C). Fig. 5 ELISA analysis of selected secreted protein levels in tissue samples. Box plots representing the concentration of PLA2G2A ( A ), ANGPTL7 ( B ) and PLA2G5 ( C ) measured by ELISA in tissue samples across each endometriosis subtype and control endometrium. The y-axis represents protein concentration (as indicated as y-axis title), while the x-axis represents the sample groups. Statistical significance between groups was determined using Wilcox test and actual p-values are displayed. Alt text: Box plots showing ELISA-based quantification of selected secreted proteins in tissue samples from different endometriosis subtypes and control endometrium. ELISA analysis of selected secreted protein levels in tissue samples. Box plots representing the concentration of PLA2G2A ( A ), ANGPTL7 ( B ) and PLA2G5 ( C ) measured by ELISA in tissue samples across each endometriosis subtype and control endometrium. The y-axis represents protein concentration (as indicated as y-axis title), while the x-axis represents the sample groups. Statistical significance between groups was determined using Wilcox test and actual p-values are displayed. Alt text: Box plots showing ELISA-based quantification of selected secreted proteins in tissue samples from different endometriosis subtypes and control endometrium.

This study was approved by local Ethics committee of General University Hospital, Department of Gyneacology, Obstetrics and Neonatology (Ethical approval ID: Ref. No. 3/23 Grant GIP) on 16/02/2023. The recruitment of the patients for the study started on 01/04/2023. The clinical trial was registered at Clinicaltrials.gov under ID number NCT06020482 ( https://clinicaltrials.gov/study/NCT06020482 ), date of registration: 28.8.2023. Clinical and demographic characteristics of the participants in the study are presented in Table 1 . A total of 81 individuals were enrolled in the study. In the RNA-seq cohort, we recruited 11 patients with ovarian endometriosis (OE), 11 patients with deep endometriosis (DE), 5 patients with peritoneal endometriosis (PE), and 15 controls withouth diagnosed endometriosis. In the ELISA validation cohort, we recruited 10 patients with OE, 14 patients with DE, 9 patients with PE, and 11 controls withouth diagnosed endometriosis. Additionally, 4 individuals withouth diagnosed endometriosis and 1 OE patient from the RNA-seq cohort were included in the ELISA validation cohort. All patients as well as control individuals were recruited at the General University Hospital, Department of Gynecology, Obstetrics, and Neonatology, and provided written informed consent. Endometriosis was diagnosed according to international and Czech guidelines for the diagnosis of endometriosis. The surgical staging of endometriosis was performed according to the American Society for Reproductive Medicine (ASRM) classification system 77 . During laparoscopy, lesion size, depth of infiltration, location, and the extent of adhesions were systematically recorded, and each parameter was assigned a weighted score. The cumulative score was then used to stratify patients into disease stages I–IV: minimal (I), mild (II), moderate (III), or severe (IV). Samples of ectopic as well as normal eutopic endometrium were excised and immediately flash-frozen in liquid nitrogen and stored at -80 °C until processing. Frozen tissue was transferred into the pre-chilled mortar and liquid nitrogen was added. The tissue was ground into a fine powder using a pestle, while liquid nitrogen was replenished as necessary to maintain the frozen state. The homogenized powdered tissue was transferred into sterile, pre-labeled microcentrifuge tubes and stored at -80 °C until further processing. Frozen powdered tissue samples were homogenized in TRIzol™ reagent (Invitrogen) by vigorous pipetting to ensure complete disruption, and total RNA was extracted following the manufacturer’s protocol. The concentration and purity of the RNA were assessed using a NanoDrop™ spectrophotometer (Thermo Fisher Scientific), and integrity was verified by agarose gel electrophoresis. The extracted RNA was stored at -80 °C until further use. All procedures were carried out under RNase-free conditions to prevent RNA degradation. RNA library preparation was performed using the NEBNext ® Ultra™ II Directional RNA Library Prep Kit (NEB#E7760L) in conjunction with the NEBNext Poly(A) mRNA Magnetic Isolation Module (NEB #E7490), following the manufacturer’s instructions. Briefly, poly(A)-containing mRNA was isolated from total RNA using magnetic beads conjugated with oligo(dT). The captured mRNA was then fragmented at 94 °C for 15 min in first-strand synthesis buffer. First-strand cDNA synthesis was performed using random primers and reverse transcriptase, followed by second-strand synthesis incorporating dUTP to maintain strand specificity. The resulting cDNA was purified using SPRISelect beads (Beckman Coulter) and subjected to end repair, adapter ligation, and USER enzyme treatment. Library amplification was performed using index primers, and the final libraries were purified and assessed for quality and concentration using a Bioanalyzer 2100 (Agilent). The libraries were quantified using a Qubit™ fluorometer (Thermo Fisher Scientific) and pooled for sequencing. RNA sequencing was performed using the Illumina NovaSeq 6000 platform with NovaSeq 6000 S1 Reagent Kit version 1.5 (total 138 sequencing cycles). The prepared RNA libraries were pooled at 2.5 μm and loaded onto the NovaSeq 6000 sequencing system following the manufacturer’s guidelines. Cluster generation and sequencing were conducted using the patterned flow cell technology of the NovaSeq platform. Paired-end sequencing (P7 index: 8 cycles, P5 index: 8 cycles, Read1: 70 cycles, Read2: 52 cycles) was aimed at an average 30 million reads per library to ensure sufficient read depth and coverage. Real-time data processing and base calling were carried out using the Illumina Real-Time Analysis (RTA) software. The raw sequencing data were demultiplexed and converted into FASTQ format using Illumina’s bcl2fastq conversion software locally. Quality control of the raw sequencing reads was assessed using FastQC, and adapter trimming and filtering of low-quality reads were performed using cutadapt. The processed sequencing reads were then aligned to the human reference genome GRCh38 using STAR aligner (version 2.7.10a) 78 . Subsequent downstream analyses including read count filtering, normalization, dispersion estimation, variance stabilization and identification of differentially expressed genes (DEGs) were carried out in R (4.4.2) using Bioconductor package DEseq2 79 . DESeq2 modeled the counts using a negative binomial distribution and expression of normalized data was adjusted utilizing a Bayesian shrinkage estimator package apeglm 80 . Identification of DEGs was preceded by the Wald test and the final P values were adjusted using the Benjamini and Hochberg method. Genes with log2 fold change > 1 and p-adjusted value < 0.05 were considered significantly differentially regulated and for pathway analyses DEGs with p-adjusted value < 0.05 were included. For ELISA sample preparation, previously frozen and ground tissue samples were lysed in phosphate-buffered saline (PBS) using sonication followed by a single freeze-thaw cycle. Briefly, the powdered tissue was suspended in cold PBS at a ratio of 1:10 (w/v) and subjected to probe sonication on ice in short bursts to prevent overheating. Following sonication, the lysates were subjected to one freeze-thaw cycle by freezing at -80 °C and rapidly thawing on ice to ensure complete lysis. The lysed samples were then centrifuged at 12,000 × g for 10 min at 4 °C to remove debris, and the supernatants were collected and stored at -80 °C until analysis. The ELISA assay was performed according to the manufacturer’s instructions for the specific target of interest (PLA2G2A – Invitrogen, EH369RB; PLA2G5 – CUSABIO, CSB-EL018103HU; ANGPTL7 – Invitrogen, EH30RB). Briefly, samples and standards in diluent buffer were added to 96-well plates pre-coated with capture antibodies and incubated overnight at 4 °C or two hours at 37 °C. Afterwards, samples were removed and incubated with biotinylated detection antibodies incubated for 1 h, followed by washing steps ( n  = 5) and incubation with streptavidin-HRP conjugate. The colorimetric reaction was developed using TMB substrate and stopped with STOP solution. Absorbance was measured at 450 nm using a microplate reader, and the protein concentration was determined by interpolation from a standard curve. Pathway enrichment was performed using Cytoscape software 81 , 82 . Gene sets of interest were imported, and enrichment was assessed using standard functional annotation and network analysis workflows. The results were visualized and further illustrated using the SRPlot web tool 83 , which provides customizable graphics for enrichment data. The relative abundance of tumor-infiltrating immune cell populations was estimated using the CIBERSORT algorithm 84 . Gene expression profiles were analyzed against the LM22 leukocyte gene signature matrix, which defines 22 distinct immune cell subsets. The analysis was run with 1,000 permutations to improve robustness. The results were further processed and visualized in excel where stacked bar plots were generated to illustrate the distribution and relationships of immune cell subsets across samples. Additionally, enrichment scores for selected cell populations were visualized using SRPlot webtool 83 and statistically evaluated using Wilcox test. We applied xCell method 28 to estimate the relative enrichment of immune and stromal cell populations from bulk transcriptomic data. xCell is a gene signature–based method that infers 64 immune and stromal cell types by integrating single-sample gene set enrichment analysis with spillover compensation to reduce dependencies between closely related cell types. Normalized gene expression matrices (log2-transformed transcripts per million, TPM) were used as input. Expression data were uploaded to the xCell web server ( https://xcell.ucsf.edu/ ) and enrichment scores for each cell type were calculated. Scores are non-absolutely scaled and reflect relative enrichment within each sample, rather than absolute cell fractions. For data interpretation we only retained non-immune cell types with enrichment scores above 0.05 present at least in two sample types (Eutopic endometrium, Ovarian endometriosis, Deep endometriosis, Peritoneal Endometriosis). xCell scores were subsequently used for comparative analyses between sample types. Exploratory data analyses, including heatmaps, principal component analysis (PCA), and volcano plots, were generated in R (R Foundation for Statistical Computing, Vienna, Austria). Heatmaps were created to visualize relative expression patterns across samples, PCA was applied to assess sample clustering and variation, and volcano plots were used to highlight differentially expressed genes/proteins based on fold change and statistical significance. Tables were prepared in excel using conditional formatting rules to indicate fold changes of highlighted genes, box plots were prepared suing SRPlot web tool 83 . Group comparisons scores for patient characteristics were treated as continuous variables, consistent with prior practice, and analyzed using Student’s t-test. Unless otherwise specified, a two-tailed test was applied. A p-value < 0.05 was considered statistically significant. Group comparisons of cell population enrichment and ELISA protein concentrations were performed using the non-parametric Wilcoxon rank-sum test and p-values are reported.

In conclusion, our comprehensive, patient samples-centric RNA sequencing analysis identified several known as well as novel genes that are overexpressed in endometriotic lesions. We identified increased expression of members of phosphodiesterase family (PDEs) in endometriotic lesions, which has not been previously published. Involvement of these enzymes in fibrotic signaling indicates their potential as therapeutic targets for endometriosis. Furthermore, we identified increased expression of ANGPTL7 in extra-ovarian endometriosis and confirmed its elevated protein levels in both deep infiltrating and peritoneal lesions. These results highlight its potential as a subtype-specific biomarker for endometriosis. Based on these findings, future studies will focus on further validating the diagnostic and therapeutic utility of the identified targets.

Despite being described more than a century ago, endometriosis still presents significant challenges in its management, starting with accurate and timely diagnosis and extending to curative treatment. The current diagnostic algorithm is cumbersome, involving multiple approaches, often including invasive procedures such as laparoscopic surgery 63 . Moreover, on average, it takes approximately 8 years to reach a conclusive diagnosis, which often results in advanced disease that is much more difficult to manage compared to disease in its early stages 64 . Although basic mechanisms associated with endometriosis such as inflammation and fibrotic changes are characterized, underlying cause and risk factors are largely underexplored 7 . Furthermore, endometriosis is not homogeneous condition, but at least 3 major types of endometriosis are postulated: ovarian, peritoneal and deep infiltration endometriosis 1 . These subtypes differ in anatomical localization, clinical presentation, and possibly in their pathogenesis and molecular signaling pathways 7 . Given this complexity, there is a pressing and unmet need to develop non-invasive diagnostic approach based on circulating biomarkers that could enable not only earlier detection of endometriosis but also the discrimination between its subtypes. Recently this area of biomedical research is gaining traction and several groups published reports on developing diagnostic tests based on blood biomarkers 14 , 18 . Here we aimed to analyze the expression profiles across subtypes of endometriosis in order to identify potential biomarkers that could aid in development of improved diagnostic algorithms. We analyzed 26 samples of ectopic endometrium along with 14 samples of control healthy endometrium using unbiased RNA sequencing. Overall, we identified major changes in the gene expression profiles among different types of endometriosis in comparison to healthy endometrium as well as between various types of endometriosis. Upregulation of genes belonging to established pathways governing endometriosis pathophysiology such as inflammation, fibrosis or steroid hormone synthesis were used as validation of our approach and results. Additionally, we identified novel pathways and upregulated genes in individual types of endometriosis. One of the key undescribed pathways that was identified in our dataset is phosphodiesterase signaling that was greatly enriched in DE as well as in PE, represented by PDE1A, 1B and 1 C as significantly upregulated in both DE and PE, while other classes (PDE2, PDE3 and PDE5) are upregulated only in DE. Although it is established that phosphodiesterase activity is crucial for stimulation of fibrosis and contractility 65 , there is very limited research conducted in the role of PDEs in endometriosis. Another upregulated gene cluster included members of the calcium dependent phospholipase A2 family, such as PLA2G2A, PLA2G2F or PLA2G5. These are secretory phospholipases that are crucial for prostaglandins production and inflammatory signaling 66 – 68 . Given their role, it is likely that these phospholipases could play an important role in mediation of inflammatory signaling spread thus facilitating recruitment of immune cells as well as fibroblasts characteristic for endometriosis. On the other hand, there were upregulated pathways highly specific for OE such as increased MHC class II receptor signaling, increased vesicular transport or steroid synthesis pathway. Interestingly, small nucleotide polymorphisms in MHC II receptors (HLA -DQB1, HLA-DRB1) 69 were associated with increased risk of endometriosis, and our analysis showed that expression of the same HLA genes (HLA -DQB1, HLA-DRB1) is significantly upregulated in ovarian endometriosis. These expression changes are well correlated with immune cell infiltrations. Increased expression of PLA2G-family of phospholipases is correlated with increased polarization of macrophages towards M2 phenotype which corresponds to described role of PLA2G5 in M2 polarization of macrophages 40 . Moreover, we can see that plasma B-cells are significantly enriched specifically in ovarian endometriosis which is aligned with increased MHC class II expression. On the contrary, we see high degree of NK-cell depletion across all types of endometriosis, which could be attributed to increased M2 macrophage polarization and increased prostaglandin production through action of PLA2G-family of phospholipases 70 , 71 . In conclusion, we identified several potential molecular and cellular targets such as PDEs, PLA2Gs or M2 macrophages that could be used to develop targeted therapies for endometriosis. Additionally, we aimed to identify potential secreted biomarkers that could facilitate development of improved diagnostic algorithm for endometriosis in general as well as for individual subtypes. First, we see increased expression of several previously identified genes (SFRP2, CCL11, LEP, C7, PLA2G5, PLA2GA) whose protein products are secreted thus could be used for endometriosis diagnosis. SFRP2, PLA2G5 and C7 expression was shown to be increased in endometriosis however their presence was not conclusively investigated in the blood of patients 58 , 61 , 72 . CCL11, LEP and PLA2GA were also shown to be consistently overexpressed in endometriosis, however increased presence of these proteins in serum of patients was not confirmed 19 , 62 , 73 . More importantly, we have identified several novel secreted proteins whose expression is greatly upregulated across all subtypes (FIBIN) or specifically in extraovarian endometriotic lesions (ANGPTL7, GDF6, ADIPOQ, BMP5). These proteins were not yet directly linked with endometriosis however play important role in endometriosis-associated pathways such as contractility regulation, ECM remodelation or TGFβ signaling 74 – 76 which suggests potential for further exploration as biomarkers for endometriosis. Although our findings are consistent with previous studies, several limitations should be acknowledged. First, the use of bulk RNA-seq may include residual host tissue admixture despite careful surgical excision and histological validation; to address this, we performed in silico deconvolution analyses, which confirmed enrichment of cell populations characteristic of diseased tissue or control epithelium. Second, the exploratory nature of the study restricted validation to a focused ELISA panel without targeted immunohistochemistry. Finally, confirmation in independent cohorts, including evaluation of circulating biomarkers in blood or other body fluids, will be essential to establish the translational potential of the identified candidates.

Endometriosis is a chronic, progressive gynecological condition characterized by the abnormal presence of endometrial glands and stroma outside the uterus, resulting in chronic inflammation within the peritoneal cavity. It primarily affects pelvic organs such as the ovaries, fallopian tubes, bladder, intestines, and peritoneum and may occasionally present in extrapelvic sites such as the diaphragm, pleura, abdominal wall, or even the central and peripheral nervous systems 1 . Globally, endometriosis affects approximately 10% of reproductive-aged women (~ 190 million) 2 . Despite its prevalence, the diagnosis is frequently delayed by 8–10 years 3 . The condition also imposes a significant economic burden, with estimated direct medical costs ranging from US $1459 to US $20,239 per patient annually, and indirect costs reaching US $4572 to US $14,079 (2022) 4 . Endometriosis presents with a wide range of manifestations, from asymptomatic lesions discovered incidentally to severe cases. Symptoms do not always correspond with the size of the lesions and typically appear before the age of 20 5 . The primary symptoms include chronic pelvic pain, dysmenorrhea, dyspareunia, dysuria, dyschezia, abdominal bloating, and constipation. Additionally, endometriosis may increase the risk of mental health issues such as anxiety and depression. The disease significantly impacts patients’ quality of life and social well-being. Due to pain and other symptoms like fatigue, heavy bleeding, or mood swings, women may miss school or work and may also avoid sexual activity 6 . Endometriosis is primarily classified based on its localization and histopathology, with three main subtypes: superficial peritoneal endometriosis, ovarian endometriosis, and deep endometriosis. Superficial peritoneal endometriosis is the most common form, accounting for 70–80% of diagnosed cases and it occurs on the surface of pelvic organs and often adheres to the peritoneum. Ovarian endometriosis affects the ovaries, forming cystic structures known as endometriomas or “chocolate cysts.” This subtype is frequently linked to infertility and carries an elevated risk of ovarian malignancy. Deep endometriosis invades visceral organs to a depth of 5 mm or more, either within or beyond the pelvic cavity, leading to anatomical distortion and severe symptoms that frequently require surgical intervention 7 . Although DE is present in only ~ 1% of women with endometriosis in the general population, its prevalence increases to 21.8% among women with subfertility 8 . Moreover, different subtypes often coexist. It is estimated that 20–50% of affected individuals exhibit multiple lesion types; for instance, up to 50% of women with ovarian endometriosis also have deep lesions 8 . Although endometriosis was first described over a century ago 9 , 10 , diagnostic methods remain limited, relying on medical history, physical examination, imaging, laparoscopy, and histopathological confirmation. Diagnostic accuracy varies significantly depending on lesion location, disease severity, and clinician expertise 11 . Laparoscopy remains the gold standard for diagnosis, but it is an invasive, costly procedure with potential risks, including nerve damage, injury to pelvic organs or major blood vessels, postoperative complications, and adhesion formation. In developed countries, expert ultrasound may be sufficient for the diagnosis of endometriosis; however, globally, where such diagnostic methods may not be available, laparoscopy remains the diagnostic standard 12 . Imaging techniques like transvaginal ultrasound and magnetic resonance imaging (MRI) can detect ovarian endometriosis and deep endometriosis lesions, yet non-invasive identification and diagnosis of endometriosis – particularly its subtypes - remains a significant challenge 12 . Non-invasive diagnostic methods, such as blood biomarker analysis, are well-established for other conditions. However, despite extensive research, no reliable blood-based test for endometriosis diagnosis has been developed 13 . To this date, there are biomarkers whose levels has been shown to be elevated in serum of patients with endometriosis such as CA-125, BDNF, IL-6, IL-8, TNF-a, or IL-32 14 – 17 . However, there is a significant variability in the outcome of the studies investigating the utility of these biomarkers for diagnosis of endometriosis and specificity and sensitivity of individual or combined biomarkers ranges between 60% and 80% 18 . Moreover, majority of these biomarkers are associated with other inflammatory conditions which further impacts their clinical utility for endometriosis diagnosis 19 . To complement the current landscape of potential endometriosis biomarkers we analyzed 26 tissue samples from women suffering from various types of endometriosis (Ovarian endometriosis – 11, Deep endometriosis – 11, Peritoneal endometriosis – 5) along with endometrial 15 samples of healthy individuals using RNAseq followed by ELISA validation of selected biomarkers. Our analysis revealed significant differences in expression profiles between normal eutopic endometrium and various types of endometriosis. Moreover, we confirmed differential immune cell landscape associated with endometriosis with most significant differences in presence of NK cells and M2 macrophages. Finally, we postulate PLA2G2A, ANGPTL7 and PL2G5 as potential secreted biomarkers of endometriosis.

Below is the link to the electronic supplementary material. Supplementary Material 1 Supplementary Material 1 Supplementary Material 2 Supplementary Material 2 Supplementary Material 3 Supplementary Material 3 Supplementary Material 4 Supplementary Material 4 Supplementary Material 5 Supplementary Material 5