Deconvolution of DNA methylation data reveals cell type-specific epigenomic signatures in endometriosis and eutopic endometrium

A particularly important question is whether DNA methylation is altered between EE of EM cases and controls. We searched for previously published reports of DNA methylation generated to address this question. Although we could not find DNA sequencing-based studies, we were able to identify a large study ( Mortlock et al., 2023 ) in which DNA methylation was measured in aspirate biopsies of EE from 637 EM cases and 347 controls using Illumina’s Infinium Methylation EPIC array, which contains probes targeting around 800 000 CpGs. After adjusting for observed and latent confounds (e.g. cell type), the authors found that DNA methylation was not different between EE of cases and controls. That is, arrow 2 in Fig. 1A was absent. We confirmed this result (see Materials and methods ). We next assessed the evidence supporting the hypothesis that cell type proportions in the EE differ in EM cases and controls (arrow 1 in Fig. 1A ). This question was not addressed in the previously published analysis ( Mortlock et al., 2023 ). Since, to our knowledge, there are no published DNA methylation signatures for EE cell types that cover EPIC CpGs, we first used the reference-free approach BCconf ( McKennan and Nicolae, 2019 ) to recover latent factors from the DNA methylation data. A caveat to this approach is that, because factors are estimated without a cell type signature, we cannot be certain that factors correlated with EM represent cell type proportion. However, we have found that methylation-derived factors associated with disease states are almost always related to cell type ( Morin et al., 2020 , 2023 ; McKennan and Nicolae, 2022 ). Using this approach, we found that Factor 11 (out of 21) was highly correlated with EM status ( Fig. 1B ; Bonferroni-adjusted P -value = 7.7 × 10 −5 ), suggesting EE cell proportions vary between EM case patients and controls. Based on the difference in Factor 11’s values in case and control groups, we estimated that a minimum of 165 subjects from each group would be required to infer such a cell type factor as being EM-related with 80% power. This large sample size likely required to infer such differences between cell type proportions in EM cases and controls suggests that the difference is small and therefore unlikely to have diagnostic utility. This is confirmed in Fig. 1C , which shows that even though Factor 11 is significantly associated with EM, it does a poor job of predicting EM. We further explored the origin of Factor 11 by correlating it to EE cell proportion estimates obtained via cellular deconvolution using an in-house signature covering EPIC CpGs (see Materials and methods ). We identified five (out of 38) cell types (natural killer [NK], endothelium, breast basal epithelium, dermal fibroblast, and neuron) that were significantly correlated with factor 11 at a 5% family-wise error rate ( Supplementary Table S2 ). NK cells and endothelium exhibited by far the greatest correlation ( p adj = 6.2 × 10 −13 and p adj = 5.9 × 10 −7 ) and were both positively correlated with Factor 11 ( Fig. 1D ). These results confirm our hypothesis that Factor 11 reflects variation due to cell type and suggest NK cell and endothelium proportions are higher in EM cases. We next undertook targeted genome-wide bisulfite sequencing on DNA extracted from EE of surgically confirmed EM controls (n = 6), from EM cases (n = 6), and paired endometriotic lesions (n = 6). Solution-phase hybridization was used for genomic region targeting following enzymatic conversion to identify methylated CpG sites. This method delivers high sequencing read depth at single base resolution across 123 Mb of well-characterized regions of the human genome including promoters, exons, introns, CpG islands, CpG island shores and enhancers. The capture regions include 3.98 million CpG sites. A total of ~2.3 billion aligned read pairs across all samples were sequenced. Available sample specific clinical and demographic information is shown in Supplementary Table S3 . No statistically significant differences were identified between the ages of individuals providing paired EE/EM lesions (cases) and normal control EE (controls). To explore the molecular pathology of EM lesions at the epigenomic level, we undertook deconvolution analysis of our data to explore differences in cell type proportion estimates (see Materials and Methods ). A heat map of fractions estimated cell type proportions for different sample groups is shown in Fig. 2A . Cell types whose proportions differed between EE of cases and controls and EM lesions at a 10% FDR threshold are shown in Fig. 2B . Those whose fractions are significantly larger in EM lesions compared to EE included heart fibroblasts (FDR-adjusted P -value = 4.2 × 10 −4 ), vascular endothelial cells (FDR-adjusted P -value = 6.1 × 10 −3 ), smooth muscle (FDR-adjusted P -value = 0.024), breast basal epithelium (FDR-adjusted P -value = 0.031), and adipocytes (FDR-adjusted P -value = 0.032). Those whose fractions are significantly smaller in EM lesions compared to EE included dermal fibroblasts (FDR-adjusted P -value = 2.8 × 10 −6 ), ovarian/endometrial epithelium (FDR-adjusted P -value = 3.2 × 10 −3 ), NK cells (FDR-adjusted P -value = 0.031), and colon fibroblasts (FDR-adjusted P -value = 0.074). Identification of differentially methylated CpG sites (DMSs) focused on promoters, exons, introns, CpG islands (CGIs), CGI shores, and enhancers. We filtered for CpG sites that were significantly differentially methylated after controlling for false discovery (FDR-adjusted P -value ≤ 0.01) and that differed in methylation rate by at least 20% between EM lesions and paired EE, and between EM lesions and normal control EE (results summarized in Table 1 and Supplementary File S1 ). Clustering analysis highlighted broad differences between EM lesions, paired EE samples, and normal control EE. These analyses were carried out either using all differentially methylated CpGs ( Fig. 3A ), or just those located in CGIs and promoters ( Fig. 3B and C ). Notably, we found tighter clustering amongst EE samples from cases than normal controls, although both were clearly distinct from EM lesions. This may have been somewhat influenced by the inclusion of the control participant who provided sample EMBx20 ( Supplementary Table S1 ). Unlike the other controls, this individual presented with pelvic pain although was confirmed as not having EM via minimal surgery. It can be seen in Fig. 3A - C that this individual clustered separately from the other normal controls. PCA using only previously identified cell type-specific methylation biomarkers ( Loyfer et al., 2023 ) gave the clearest separation between EM lesions and EE regardless of the source of EE ( Fig. 3D ). We found that differentially methylated CpG sites were more methylated in EM lesions than normal control EE samples by a ratio of about 2:1 with ~62% of sites being more methylated in the former. This trend was also observed in introns (65%), CGI shores (67%), CGI shores located exclusively in promoters (63%), and enhancers (63%). CpG sites in promoters were slightly more likely to be more methylated in endometriotic lesions (55%) but this rose to 63% when promoters were considered that did not overlap with introns, exons, CGIs, CGI shores or enhancers ( Table 2 and Fig. 4A - E ). Similar observations were made when EM lesions were compared to paired EE samples, with ~69% of differentially methylated CpG sites being more methylated in the former, along with introns (~72%), CGI shores (~71%), promoters (~61%), CGI shores in promoters (67%) and enhancers (~63%) ( Table 2 ). In contrast to the genomic regions described above, CGIs showed the opposite trend, with decreased levels of CpG methylation in EM lesions compared to control EE. Approximately 32% of CGI sites were more methylated in EM lesions compared to EE from controls. This decreased to 24% for CGIs located in promoters and 20% for CGIs located in exons. The differential methylation of CGIs in introns was less skewed, with only 42% being hypermethylated in EM lesions compared to EE controls ( Table 3 , Fig. 4F - I ). Similar observations of reduced methylation were made when EM lesions were compared to paired EE samples from case patients, with CGI (35%), CGI in promoters (22%), CGI in exons (26%), and CGI in introns (45%) ( Table 3 ). To gain further insight into EM associated DMSs we used Ingenuity Pathways Analysis (IPA) software (Germantown, MD, USA) to identify biological pathways associated with genes that were differentially methylated between EM lesions and EE from control patients. We first considered DMSs located in CGIs and found that these CGIs were most likely to be located within genes associated with GABA receptor signaling ( P = 7.01 × 10 −5 ), axonal guidance signaling ( P = 7.44 × 10 −5 ) and CREB signaling in neurons ( P = 3.66 × 10 −4 ) ( Supplementary File S2 ). Further analysis revealed similar pathway associations for genes containing DMSs in introns, exons and CGI shores. For example, we found that DMSs in introns were commonly located within genes associated with axonal guidance signaling regardless of whether we focused on DMSs that were hypomethylated ( P = 2.92 × 10 −14 ) or hypermethylated ( P = 3.26 × 10 −25 ) in lesions when compared to normal eutopic endometrial tissue ( Supplementary File S2 ). Similarly, DMSs in CGI shores were most likely to be in genes associated with axonal guidance signaling ( P = 1.03 × 10 −16 ), and this was also true when only CGI shores located in promoters were considered ( P = 5.33 × 10 −8 ) ( Supplementary File S2 ). When we focused on DMSs located in the CGIs of promoters, identified by comparing EM lesions with EE from controls, we made several interesting observations with respect to neuronal biology. For example, we found that a CGI located in the promoter of neuronal marker ASCL2 displayed reduced methylation in EM lesions when compared to EE ( Fig. 5A ). Similar patterns of reduced methylation in EM lesions were observed for other neuronal markers including BARHL2 , CLDN11 , FARP2 , ( Fig. 5B - D ). Other notable genes with altered methylation patterns include GATA2 , HAND2 (increased methylation in EM lesions), and TCF21 (reduced methylation in EM lesions) ( Fig. 5E - G ). We next undertook an analysis of upstream regulatory mechanisms associated with genes in the ontological gene networks we identified. By far, the most prominent finding was that beta-estradiol was frequently found to be a major upstream regulator of network enriched genes. This was true for most of the genomic regions we considered. For example, beta-estradiol was the top predicted upstream regulator for genes containing DMSs in introns ( P = 1.75 × 10 −23 ), CGI shores ( P = 1.56 × 10 −42 ), promoters ( P = 1.61 × 10 −28 ), and exons ( P = 4.43 × 10 −15 ) ( Supplementary File S3 ). However, this was not the case for genes containing DMSs in CGIs, for which the top two predicted upstream regulators were CTCF ( P = 3.42 × 10 −13 ) and SOX2 ( P = 2.03 × 10 −12 ) ( Supplementary File S3 ). Included in the list of genes containing DMSs in CGIs that are predicted to be regulated by CTCF are those of the PCDH-GAMMA gene cluster ( PCDHG ) which is located on chromosome 5. Predicted CTCF-regulated genes containing DMSs within their CGIs include PCDHG1 , 2 , 3 , 5 , and 7 and PCDHGA5 , A6 and A10 . Further analysis revealed that this region contains a very large number of CGIs in multiple members of the protocadherin gene family that were hypomethylated in EM lesions when compared to EE from controls ( Table 4 and Supplementary File S4 ). This is notable because protocadherins encode neural cell adhesion proteins that function to mediate neuronal survival and dendrite patterning ( Ing-Esteves et al., 2018 ). This observation aligned with our finding that genes containing DMSs are frequently associated with ontological pathways related to neuronal biology (see above).

We reanalyzed data from Mortlock et al. (2023) , who used Illumina’s Infinium Methylation EPIC array (San Diego, CA, USA) to measure the DNA methylation at 759 345 CpG dinucleotides from the aspirate biopsies of the eutopic endometrium from 637 EM cases and 347 controls. Since there is no published cell type-specific signature matrix for EE cell types covering EPIC CpGs, we took a reference free approach to deconvolving these data. We used the R package CorrConf ( McKennan and Nicolae, 2022 ) to estimate there were 21 latent factors after accounting for disease (EM case or healthy control) and the covariates batch number, recruitment site, methylation plate, ancestry (European or non-European), and menstral cycle phase. We used the R package BCconf ( McKennan and Nicolae, 2019 ) to recover latent factors while protecting disease state and adjusting for covariates, used the R function BCconf::Sig.Omega to test the null hypothesis that latent factors were independent of EM (after adjusting for covariates) against the 2-sided alternative, and the R function BCconf:: Correction to regress methylation levels onto disease state while adjusting for covariates and latent confounds (arrow 2 in Fig. 1A ). Significant CpGs were determined by controlling the false discovery rate (FDR) at 0.05 using the Benjamini–Hochberg procedure (as implemented by the R function p.adjust). To interpret latent factors related to EM, we used publicly available whole genome DNA methylation data ( GSE186458 ) to build a 38-cell type signature matrix covering 99457 EPIC CpGs. We excluded the cell type “Adipocytes” because it was estimated to make a trivial contribution in our deconvolution of whole genome EE data ( Fig. 2 ) and we could not identify any EPIC CpGs capable of distinguishing adipocyte and colon fibroblasts. We used the R package nnls to regress EPIC array beta values ( Du et al., 2010 ) onto the signature via non-negative least squares to estimate cell type proportions. Estimates were subsequently normalized to sum to one. We identifyied cell types correlated with EM-related latent factors by using robust regression (Huber’s loss using the R function MASS::rlm) to regress latent factors onto cell proportions while controlling for the above covariates. Robust regression was used to mitigate the influence of extreme or innacurately estimated proportions. We adjusted for multiple testing by controlling the FDR at 0.05 via the Benjamini–Hochberg procedure (as implemented by the R function p.adjust). All analyses were done using R version 4.3.1, R Foundation for Statistical Computing, Vienna, Austria. Tissue specimen collection and procurement was approved by the University of Pittsburgh Institutional Review Board (IRB) (Protocols 19110186 for controls and 19070028 for cases). All study subjects were consented at University of Pittsburgh Medical Center (UPMC) Magee-Womens Hospital at routine preoperative surgical consult office visits with a minimally invasive gynecologic surgeon in the Division of Gynecologic Specialties. Recruited cases included patients going to the operating room (OR) with an indication for pelvic pain. Recruited controls included patients going to the OR for benign gynecologic procedures without a chief complaint of pelvic pain nor a history of pelvic pain/endometriosis, with the exception of the participant who provided sample EMBx20 ( Supplementary Table S1 ) who presented with pelvic pain but was confirmed as not having EM via minimal surgery. Patients were excluded if they were < aged 18years or menopausal, used systemic immunosuppressants within the past 3months, took antibiotics within the past 1 month, or had a history of pelvic inflammatory disease/tuboovarian abcess (PID/TOA) or cancer. Cases were excluded if no endometriosis was found at time of surgery. Similarly controls were excluded if endometriotic lesions were found at time of surgical intervention. Endometrial biopsy aspirates were collected following induction of anesthesia and prior to the start of the planned surgical case. EM lesions were specimens resected by the surgeon for usual surgical care/intervention. Tissue samples were frozen and stored at −80°C until further analysis. DNA was extracted from tissue samples using the Nucleospin Tissue XS Kit (Macherey-Nagel, Duren, Germany). Extracted DNA was quantified by using the KAPA hgDNA Quantification and QC Kit (Roche). DNA was sheared with a Covaris S220 sonicator (Woburn, MA, USA) to a size of ~180bp. Libraries were prepared using the Enzymatic Methyl-Seq Library Preparation Protocol (NEB, Ipswich, MA, USA) and amplified 8–10 cycles. Converted DNA libraries were pooled and underwent targeted capture using the Twist Human Methylome Panel (Twist Biosciences, San Francisco, CA, USA), hybridized for ~16 hours and amplified 6 cycles. Samples were sequenced on a NovoSeq 6000 (Illumina) with 150bp paired-end reads to a mean targeted read depth of ~84x. DNA sequence reads were quality trimmed and adaptor sequences were removed using Trim-Galore ( https://www.bioinformatics.babraham.ac.uk/projects/trim_galore/ ). The reads were aligned to the human reference sequence (GRCh38/hg38) using Bismark Bisulfite Read Mapper10 ( Krueger and Andrews, 2011 ) in paired-end, bowtie2 mode. Any unmapped reads were aligned in single-end, bowtie2 mode. Read duplicates were removed using Bismark. Methylation was called on paired-end and single-end files and then merged to determine the methylation status for each CpG site. Sequencing data were aligned to the reference human genome GRCh38. We then extracted the methylation information using software Bismark Bisulfite Read Mapper ( Krueger and Andrews, 2011 ). To identify differentially methylated CpGs between groups of samples, we used negative binomial test, as implemented in the R package methylKit. Data were then filtered to determine CpG location with respect to specific types of genomic element (promoters, CpG islands, introns exons etc). Principal component analysis (PCA) was performed on the CpGs located within the promoter regions (defined as upstream 1500bp to downstream 500bp of a transcriptional start site [TSS]). Specifically, we calculated the standard deviation of the CpG methylation rates across all 18 samples for each CpG, and only CpGs with a standard deviation above the median were selected. Pairwise plots of the first 3 principal components were generated to examine the similarity among the samples. For the comparison between EM lesions and paired EE of EM patients, the FDR was controlled at the level of 0.01 ( Akalin et al., 2012 ). Available cell type-specific whole genome DNA methylation data ( GSE186458 ) was used to estimate cell type proportions in each sample. Briefly, we used the cell type reference atlas from Loyfer et al. (2023) . This data set consists of 205 tissue samples for 39 cell types, with 250 CpG marker regions per cell type. Markers with fewer than 30 counts in any sample were excluded from the analysis, which resulted in 2800 markers. Non-negative least squares were used to regress observed methylation fractions at cell type maker regions onto the cell type signature matrix. Resulting estimates were normalized to sum to one. The reanalysis of EE methylation data from the Mortlock et al. (2023) study (see above) indicated any differences in EE cell type proportions between cases and controls were too small to be detected with our sample size (see Results ). This was confirmed by an unpaired t-test comparing the EE of cases and controls, which revealed no cell types with significantly different cell type proportions (smallest FDR-adjusted P -value = 0.73). We therefore pooled EE case and control samples to compare the difference in proportions between the EE and the lesions of EM cases. To account for the many proportions equal to zero, proportions were modeled as a mixture of a point mass at zero and a log-normal distribution with a region- (EE or lesion) specific location and constant scale. The mixture weight, i.e., the probability a proportion was zero, was also region-specific. The correlation between proportions from the same subject was modeled using a random-effects probit model for binary variables indicating whether a proportion was zero, as well as a random-effects model for non-zero proportions. A P -value for the null hypothesis that the distributions of a cell type’s proportions were the same in each region type was computed using a likelihood ratio test. The FDR was used to adjust for multiple testing. Results generated using this approach are shown in Fig. 2 .

Biomarker discovery is an important step towards identifying means of early EM detection, the development of novel therapeutics, and the overall improvement of clinical care. Given the lack of diagnostic methods for EM detection, and the appeal of identifying biomarkers that can be detected either non-, or minimally invasively, we were motivated to identify epigenomic features of endometrial aspirate biopsies that might be associated with, and therefore biomarkers of, EM. Our reanalysis of existing DNA methylation data ( Mortlock et al., 2023 ) collected from eutopic endometrium and analyzed on Illuminma EPIC arrays identified one latent factor whose levels were significantly higher in EM cases. This factor was highly positively correlated with, and therefore reflects variation in, NK and endothelial cell proportions. Taken together, these results suggest NK and endothelial cells are more abundant in the EE of EM cases compared to controls, which is consistent with the observation that NK cells are increased in the perinatal fluid of women with EM ( Gomez-Torres et al., 2002 ). The increase in endothelial cells could be due to the vascularization required by EM lesions ( Laschke and Menger, 2018 ). Despite the statistically significant difference in the levels of this cell type-related factor between EM cases and controls, the absolute difference was likely too small to have diagnostic utility. Deconvolution analysis of our DNA sequencing data using the same publicly available cell type-specific DNA reference libraries atlas ( Loyfer et al., 2023 ) revealed several interesting alterations in cell type proportions across the three sample types. Surprisingly, in EM lesions, we found almost no evidence of cellular DNA containing endometrial epithelium-specific methylation biomarkers whereas these were, not surprisingly, abundant in EE samples from both case and control patients. EM lesions are thought to contain EE-derived squamous and glandular epithelial cells. Squamous EE cells do not express EpCAM, but it is expressed in the glandular epithelial cells of the endometrium. DNA methylation reference data used to identify markers of endometrial epithelium was, in fact, generated from a combination of 4 whole genome bisulfite sequencing libraries generated via EpCAM selection of EE samples (n = 3 donors) and ovarian surface epithelium (n = 1 donor) ( Loyfer et al., 2023 ). However, the latter is thought not to express EpCAM ( Szotek et al., 2008 ) so the reference data likely only contain genomic DNA from EpCAM expressing EE epithelial cells. Our failure to detect DNA methylation biomarkers for endometrial epithelial cells in EM lesions is surprising. It could be explained by there being a relatively low cell type proportion of genomic equivalents from endometrial epithelial cells in our EM tissue samples, making their detection difficult. This seems unlikely, however, because the proportion of these cells in EM lesions we studied was at or near zero ( Fig. 2B ). Also reduced in EM lesions were cell type proportions for NK cells. NK cells have been widely studied in EM with conflicting findings about their cell number and activity ( Yang et al., 2024 ). However, their reduced tissue number may correlate with lower infiltration into endometriotic tissue and impaired elimination of endometriotic cells ( Sciezynska et al., 2019 ). Not surprisingly, data from EM lesions showed relatively high levels of specific methylation biomarkers for vascular endothelium, with lower signal in EE. These findings are consistent with the known angiogenic phenotype of EM ( Laschke et al., 2011 ). An angiogenic molecular phenotype is supported by the elevated levels of smooth muscle cells observed in EE lesions. Also elevated, and possibly consistent with angiogenesis, were methylation biomarkers for heart fibroblasts. Notably, the suppression of HAND2 expression, which occurs in EE lesions relative to EE, has been shown to induce fibroblast growth factor 1 (FGF1), FGF2, and FGF9 in endometriotic stromal cells with a consequent increase in the migration and invasion capacity of endometriotic stromal cells ( Kato et al., 2019 ). As noted above (and Fig. 5 ), HAND2 DNA methylation was elevated in EM lesions relative to EE. Whether the published reference data used ( Loyfer et al., 2023 ) for deconvolving our data is specific for heart fibroblasts is unclear, but a fibrotic molecular phenotype is well-described in EM lesions ( Vissers et al., 2024 ). Furthermore, elevated cell type proportions for heart fibroblasts were mirrored by reduced levels of methylation biomarkers for other fibroblast cell types, specifically those matching dermal and colon fibroblasts. It is possible that these are surrogates for endometrial fibroblasts, which were not specifically represented in the corpus of DNA methylation data used for deconvolution ( Loyfer et al., 2023 ). When comparing EM to EE, differentially methylated CpG sites were found in genes whose functions are enriched within biological pathways associated with central nervous system biology and which are regulated by beta-estradiol-mediated signaling cascades. We identified several genes containing CGIs that were hypomethylated in lesions compared to paired case patient EE. Specifically, ASCL2 is involved in neural stem cell proliferation ( Liu et al., 2019 ). It is also involved in the control of intestinal stem cell fate and regeneration ( van der Flier et al., 2009 ; Murata et al., 2020 ) and trophoblast lineage determination during placental development ( Varberg et al., 2021 ). This is of interest not only because of the significance of ASCL2 in neuronal biology but also, perhaps, because of broad similarities between the intestine, placenta and endometrium. All are stem cell rich and undergo differentiation and proliferation. Also of note is BARHL2 , which we found to be hypomethylated in lesions. BARHL2 is a transcription factor whose expression is associated with neuronal differentiation ( Teratani-Ota et al., 2016 ). Similarly, CLDN11 , which we found to be hypomethylated in lesions, is a component of tight junctions within central nervous system myelin ( Denninger et al., 2015 ), FARP2 which was also hypomethylated in lesions, is involved in dendrite morphogenesis ( Danelon et al., 2020 ) and FOXD4L1 (also hypomethylated in lesions) which is essential for establishing neural cell fate and for neuronal differentiation ( Sherman et al., 2017 ). GATA2 contains a promoter CGI that was hypermethylated in lesions compared to paired eutopic endometrium. GATA2 has previously been shown to be hypermethylated and repressed in endometriotic cells ( Dyson et al., 2014 ). Also of note is HAND2 which contains a CGI shore and 3’ UTR (untranslated region) that was hypermethylated in EM lesions relative to EE. HAND2 DNA methylation has previously been shown to be altered in endometriosis and its hypermethylation may be predictive of endometrial cancer ( Multinu et al., 2020 ). The TCF21 gene, which is a tumor suppressor gene whose product regulates mesenchymal cell transition into epithelial cells ( Smith et al., 2006 ) and has previously been shown to display increased expression in EM lesions ( Ganieva et al., 2020 ), was also hypomethylated at DMSs within its promoter CGI. The identification of DNA methylation in genes associated with neuronal biology is notable because of the fact that one of the most prevalent symptoms of EM is chronic pelvic pain. Numerous reports have identified increased nerve density and altered neuronal phenotypes in EM lesions. For example, Yadav et al. (2021) used immunohistochemical methods to demonstrate that the majority of EM lesions they examined, but not EE aspirate biospies, were innervated. Similarly, Zhang et al. (2024) demonstrated an increased immunoreactive nerve fiber density located in abdominal wall EM lesions. Given that we also identified beta-estradiol as a major upstream regulator of the pathways enriched for genes containing differentially methylated CpGs, it is notable that estrogen receptor-β levels are elevated in EM tissues compared to the EE ( Lin et al., 2017 ). Estrogen plays a significant role in the secretion of chemokines from peripheral nerves, thereby facilitating the recruitment and polarization of macrophages within endometriotic tissue ( Montagna et al., 2008 ). Estrogen is also a key regulator of interactions between macrophages and neurons, driving macrophage recruitment and altered neurogenesis ( Greaves et al., 2015 ), the latter being involved in the EM-associated pain ( Liang and Yao, 2016 ). Abnormal neurogenesis is an important feature of endometriosis ( Tokushige et al., 2006 ; Demir et al., 2012 ; Morotti et al., 2014 ), and secretion of multiple nerve growth factors from macrophages can be induced by estrogen ( Kanda and Watanabe, 2003 ). Persistent activation by inflammatory cytokines from macrophages acting upon the nociceptors of peripheral nerves exacerbates neuroinflammation through the release of inflammatory neurotransmitters. This estrogen-regulated neuro-immune interaction sensitizes peripheral nerves, culminating in neuropathic pain manifestation in endometriosis ( Liang et al., 2018 ). The main limitation of this study is that cell type-specific DNA methylation reference data are sparse, especially for EE. Therefore, there is potential for cell type specific DNA methylation reference data to be insufficient for complete cell type specificity when undertaking deconvolution of bulk sample data. An example of this lack of specificity is likely demonstrated by our finding that cell type proportions for breast basal epithelium were altered both in the array based analysis of EE and sequencing analysis of EM/EE. It is likely that this signal represents an epithelial cell type (or types) that is not specifically breast-derived but that is distinct from the ovarian/endometrial epithelial signal that we found to be virtually undetectable in EM lesions. It is also likely that a significant number of the differentially methylated CpG sites we identified are the result of differences in cell type proportion between tissue samples ( Jaffe and Irizarry, 2014 ). However, we do not believe this detracts from their potential to illuminate our understanding of the molecular pathology of EM. Further investigations will improved DNA methylation reference data that is perhaps obtained via emerging methods for single cell DNA methylation analysis ( Acharya et al., 2024 ).

Endometriosis (EM), characterized by the presence of endometrial tissue outside of its normal anatomic location, is a common gynecologic disorder impacting nearly 15% of reproductive age women, and resulting in more than 50000 hospitalizations annually ( Eskenazi and Warner, 1997 ; McLeod and Retzloff, 2010 ). Commonly causing pelvic pain and infertility ( Hadfield et al., 1996 ; Eskenazi and Warner, 1997 ; McLeod and Retzloff, 2010 ; Nnoaham et al., 2011 ), recent studies now suggest that EM may also predispose women to endometrioid, clear-cell, and low-grade serous ovarian cancers—termed EM-associated ovarian cancer ( Edwards et al., 2015 ). Thus, in addition to quality-of-life issues ( Liu et al., 2023 ), EM may pose significant risks to a woman’s long-term health and well-being. Despite its prevalence, the molecular basis of EM remains incompletely understood, rendering diagnosis and treatment challenging. An emerging body of evidence points to epigenetic modifications as key contributors to the onset and progression of the disease. Among the epigenetic changes, DNA methylation—whereby a methyl group is added to the cytosine of a cytosine–phosphate–guanine (CpG) dinucleotide—has garnered substantial attention. Perturbations in this methylation equilibrium may have profound implications for genomic stability, gene expression, and chromatin structure. There is accumulating evidence that epigenetic aberrations play a role in the pathogenesis and/or the pathophysiology in EM and several differentially methylated genomic regions have been identified ( Dyson et al., 2014 ; Koukoura et al., 2016 ; Grimstad and Decherney, 2017 ; Yotova et al., 2017 ; Izawa et al., 2019 ). Numerous studies have reported altered DNA methylation in cells of eutopic endometrium (EE) and EM lesions. For instance, Dyson et al. (2014) suggested that DNA methylation is an integral component of EM, and identified a novel role for the GATA family in progesterone resistance and disease progression ( Dyson et al., 2014 ). Meanwhile, a number of groups have reported altered methylation of the homeobox gene HOXA10, which could influence endometrial receptivity and fertility ( Wu et al., 2005 ; Signorile et al., 2018 ; Elias et al., 2023 ). A major unmet need in EM is methods for routine noninvasive screening and phenotyping/classification. Because the only current way to definitively diagnose EM is via laparoscopic surgery, patients often suffer for many years (perhaps 8–10years) before their disease is identified and treated ( Ghai et al., 2020 ). Furthermore, the absence of screening methods has confounded efforts to develop successful non-surgical treatments. Alterations in the epigenomic molecular phenotype at the level of DNA methylation are intimately associated with the pathobiology of disease and exposure to environmental stressors ( Abdolmaleky et al., 2006 ; van Vliet et al., 2007 ; Ball et al., 2009 ; Good et al., 2020 ). Recent developments in liquid biopsy suggest that it may be possible to identify cell-free DNA (cfDNA) methylation signatures that are associated with EM. A potential source of these nucleic acids is EM lesions that exist outside of the uterus. These have the potential to release cfDNA into the peripheral circulation. Therefore, the identification of epigenomic dysregulation in EM is important because the further characterization of these altered molecular signatures may not only advance our understanding of the underlying pathobiology of endometriosis but may also enable novel diagnostic and phenotyping strategies in which these are targeted as biomarkers via liquid biopsy of plasma ( Bunce et al., 2012 ; Yu et al., 2015 ; Lehmann-Werman et al., 2016 ; Shen et al., 2018 ; Kisiel et al., 2019 ). Additionally, analysis of DNA methylation in samples of EE may identify biomarkers with utility for EM screening and epigenomic analysis of EM lesions may improve our understanding of the molecular pathology of EM and provide insight into inter-individual differences in disease phenotype ( Zhang et al., 2023 ). Given this, we undertook deconvolution analysis of DNA methylation data previously generated from endometrial aspirate biopsies obtained from 637 EM cases and 347 controls using Illumina’s Infinium Methylation EPIC array. Our goal was to identify cell type proportion differences between cases and controls with potential to form the basis of an EM screening assay. We further analyzed DNA methylation in EM lesions and paired EE samples via solution phase hybridization and DNA sequencing with the goal of identifying epigenomic dysregulation in EM and further advancing our understanding of the molecular pathology of EM lesions.

Supplementary data are available at Molecular Human Reproduction online.