Single-cell analysis of endometriosis reveals a coordinated transcriptional programme driving immunotolerance and angiogenesis across eutopic and ectopic tissues

This study was approved by the Ethics Committee of the Institutional Review Board at University of Connecticut Health Center (UCHC), The Jackson Laboratory, and the Human Research Protection Office of U.S Department of Defense and conducted according to all relevant ethical regulations regarding human participants. Written informed consent was obtained from all participants. All participants consented to share recoded information in public, unrestricted databases. Tissue samples were obtained from UCHC. Pre-menopausal female patients (aged 18 to 49 years old) pre-operatively diagnosed with stage II-IV endometriosis and scheduled for laparoscopic surgery were invited to participate in this study. Endometriosis staging was confirmed at the time of laparoscopy according to the revised American Society for Reproductive Medicine guidelines. The majority of the patients (in both control and endometriosis cohorts) were treated with similar hormonal treatments at the time of sample collection, and as detailed in Supplementary Table 1 . To recruit comparable controls, we chose patients under similar progestin/estrogen therapies and subject to the strict selection criteria we established: no history of inflammatory condition or cancer; presence of hormonal treatment; agematched to endometriosis cohort; visual inspection of absence of endometriosis (via laparoscopy with expert surgical evaluation of endometriosis absence). Matched eutopic endometrium and endometriosis tissues were collected from endometriosis patients ( Fig. 1a ; Supplementary Fig. 1a, b). Eutopic endometrium was obtained by performing an endometrial biopsy during hysteroscopy. Ectopic peritoneal endometriosis was obtained by resecting the entire endometriosis lesion and adjacent peritoneum ensuring the entire visible lesion was excised. For the control cohort, non-endometriosis eutopic endometrium biopsies were obtained from patients scheduled for surgery who were not suspected to have endometriosis. Complete patient demographic information is provided in Supplementary table 1 . Analysis of the endometrium histology was performed on H&E stained tissue sections by UCHC pathologists. Upon resection, fresh tissue was immediately stored in MACS tissue storage solution (Miltenyi, 130–100-008) and kept on ice until processing. Fresh tissues were immediately processed for scRNA-seq. Ectopic endometriosis lesions from the peritoneum were divided into ectopic lesion (EcP) and ectopic adjacent (EcPA) (Supplementary Fig. 1a). Viable single cells were obtained by mechanical and enzymatic digestion using cold active protease (CAP), following a modified version of the previously described protocol 75 . Briefly, minced tissue was transferred to GentleMACS C tubes (Miltenyi, 130–096-334) containing protease solution (10mg/ml Bacillus Licheniformis protease (CAP) (Sigma, P5380) in DPBS supplemented with 5mM CaCl2 and 125U/ml DNaseI (Stemcell, 07900) and incubated in cold water bath (6°C) for 7–10 minutes, performing trituration steps every 2 minutes. After incubation, sample were mechanically dissociated on a Miltenyi GentleMACS Dissociator for 1 minute, twice. Undigested tissue was allowed to settle by gravity for one minute. Single cells within the supernatant were transferred into a collection tube containing wash buffer PBS supplemented with 10% fetal bovine serum (FBS) (Gibco, 10082147), 2mM EDTA, and 2% bovine serum albumin (BSA, Miltenyi 130–091-376). Remaining undissociated tissue was incubated with fresh CAP protease for a total of 20 to 40 minutes, proceeding with a trituration step every 5 minutes and a Miltenyi gentleMACS Dissociator step every 15 minutes. After recovery of single cells, residual undissociated tissue was incubated with PBS supplemented with 1 mg/ml dispase on the Miltenyi gentleMACS Dissociator at 37°C for 15 minutes, and until complete tissue dissociation. Single cells were then pelleted, washed, and filtered through 70μm MACS Smartstrainer (Miltenyi, 130–098-462). Prior to FACS sorting, single cell suspension were stained with propidium iodide (PI) (BD Biosciences, 556364) and calcein violet (Invitrogen, C34858 ) in FACS buffer (PBS, 2mM EDTA, 2% BSA) and according to manufacturer protocols. Viable cells (propidium iodide negative and calcein violet positive) were sorted using the BD FACS Aria Fusion cell sorter, gated using FACS Diva (9.0.1) and recovered in Advanced DMEM/F12 (Gibco, 12634010) supplemented with 2mM GlutaMAX (Gibco, 35050061), 10mM HEPES (Gibco, 15630080), 20% FBS, 1% BSA. Sorted viable cells were then washed and resuspended with 0.04% BSA in PBS and assessed for viability using trypan blue staining for subsequent scRNA-seq experiments. The detailed protocol is available at Protocols.io 76 . Following tissue dissociation and single cell recovery, and after 10x chromium chip loading, remaining single cells were pelleted and resuspended in cold Matrigel (Corning, 356231). Fifty microliter (50 μL) droplet were plated onto 24-well plate wells (Greiner Bio-one, 662102) to generate endometrial epithelial organoids (EEOs). After Matrigel dome solidification, organoid media was added to cover each dome, as previously describe by Boretto et al. 77 . Organoid passaging was performed every 7–10 days and according to the established protocol from Turco et al. 78 . For scRNAseq experiments, organoid cultures between passage 3 and 5 and at day 7–11 after plating were collected, washed twice with wash media (Advanced DMEM/F12, 2mM GlutaMAX, 10mM HEPES, 0.1% BSA), and dissociated into single cells using TrypLE Express (Gibco, 12605010) for 3–5 minutes at 37°C. Cell suspensions were filtered with 40 μm mesh filter to remove debris and cell aggregates. Lastly, cells were washed and resuspended with cell staining buffer (Biolegend, 420201) for hashing with TotalSeq-A anti-human Hashtag reagents ( Supplementary Table 8 , Biolegend) for 30 min at 4°C and following previously published protocol 79 . After staining, cells were washed to remove excess antibody and resuspended in PBS/0.04% BSA for subsequent counting. Hashed cells were assessed for viability and sorted for viable cells as described below. Single cell suspensions were analyzed for viability and counted on a Countess II automated cell counter (Thermo Fisher). A total of 12,000 cells were loaded onto a channel of 10X Chromium microfluidic chips for a targeted cell recovery of 6,000 cells per lane. Single cell capture, barcoding, and library preparation were performed using 10X Chromium v3 chemistry according to manufacturer’s protocol (10x Genomics, CG000183 ). Sample cDNA and library quality controls were performed using the Agilent 4200 TapeStation instrument and quantified by qPCR (Kapa Biosystems/Roche). Libraries were sequenced on a NovaSeq 6000 (Illumina) with the S2 100 cycle kit targeting 100,000 reads per cell for tissues or 50,000 reads per cell for organoids. Illumina base call files for all libraries were demultiplexed and converted to FASTQs using bcl2fastq v2.20.0.422 (Illumina). The CellRanger pipeline (10x Genomics, version 3.1.0) was used to align reads to the human reference GRCh38.p13 (GRCh38 10x Genomics reference 3.0.0), deduplicate reads, call cells, and generate cell by gene digital counts matrices for each library. The resultant counts matrices were further processed with Scanpy package (version 1.7.1) 80 to exclude genes that are detected in less than 3 cells and to exclude cells with (1) fewer than 500 genes, (2) fewer than 1,000 UMIs, (3) maximum of 100,000 UMIs, and (4) maximum mitochondrial content of 25%. Doublet identification were performed using Scrublet 81 . Filtered matrices were then combined and normalized such that the number of UMI in each cell is equal to the median UMI across the dataset and log transformed. Scanpy was used to identify the top 2,000 highly variable genes from log transformed combined matrix. The (1) mitochondrial genes, (2) hemoglobin genes, (3) ribosomal genes, (4) cell cycle genes 82 , and (5) stress response genes were excluded from highly variable gene set 83 . Principal component analysis and neighborhood graph generation were performed based on highly variable genes set. Harmony (version 1.0) batch correction was performed to reduce variabilities introduced by inherent patient differences, tissue types, and endometriosis staging to enhance clustering by major cell type 84 . For subsequent clustering of each major cell type, batch correction was performed to account only for inherent patient differences and/or endometriosis staging to preserved tissue types specific expression 85 . Batched-corrected principal components were used for dimensionality reduction using Uniform Manifold Approximation and Projection (UMAP). Clustering was then performed with Leiden community detection algorithm 86 , 87 . Further doublet identification was calculated based on the median distance of a cell to the center of its respective cluster centroid in UMAP space and the coexpression of marker genes of two or more cell types. All suspected doublets were removed from the analysis. Marker genes of each cluster were identified using Wilcoxon Rank-Sum test in a one-versus-rest fashion, with (1) minimum 0.5 – 2 fold change between group, (2) expressed by at least 0.7 fraction of cells in the group, and (3) expressed by maximum 0.3 fraction of cells outside the group. Cell types were determined by matching the biomarkers with previously described cell types and cell states, and from biomarkers curated from the literature ( Supplementary Table 3 ). In order to assess possible biases in cellular diversity or transcript capture in the droplet-based scRNA-seq, we performed bulk RNA-seq on tissue-type-matched aliquots of 6 patients from the scRNA-seq cohort as well as 9 additional patients. In total, we collected 12 eutopic, 6 ovarian lesions, and 6 peritoneal lesions bulk RNA-seq samples ( Supplementary Table 2 ). Total RNA of endometrium and endometriotic lesions was extracted from snap frozen tissue or RNAlater (Invitrogen, #AM7020) stabilized tissue using QIAGEN RNeasy Mini Kit according to the manufacturer’s instructions. Library preparation was performed using KAPA mRNA Hyperprep kit (Roche) according to manufacturer’s instruction. Bulk RNAseq libraries were sequenced on NovaSeq 6000 (Illumina) with SP 100 cycles single-end reads kit resulting in an average of 38.9 million reads per sample. Reads were aligned to the GRCh38.p13 reference genome (GRCh38 10x Genomics reference 3.0.0), filtered, and quantified with nf-core/rnaseq (version 1.4.2) 88 utilizing the STAR aligner. Read counts were normalized to counts per millions (CPM) reads. scRNA-seq was compared to bulk RNAseq by utilizing pseudo-bulk transform (summing UMI counts for all cells in each sample and CPM normalization. For each tissue type, we compared these bulk transcriptomes with pseudo-bulk scRNA-seq profiles of the same type ( Extended data Figure 1 ). Briefly, we computed the spearman correlation for each pairwise combination (n=144 Eutopic, n=24 Ovarian, n=90 peritoneal) of bulk RNA-seq transcriptome and pseudo-bulk transcriptome, shown in Extended Data Figure 1a ; spearman correlation helps minimize unwanted biases derived from differences in total mRNA abundance and differences in normalization strategies between pseudo-bulk and bulk expression profiles. Then we computed the spearman correlation between the mean bulk and single cell pseudo-bulk expression profiles across samples sharing the same sample type (shown in Extended Data Figure 1b ). Differential gene expression between scRNA-seq and bulk RNA-seq data was analyzed with edgeR exactTest 89 . Differentially expressed genes (DEGs) were generated sequentially for eutopic endometrium (Control and EuE), ectopic peritoneal endometriosis (EcP and EcPA), and ectopic ovarian endometriosis (EcO) ( Supplementary Table 4 ). DEG analysis between tissue types within a population was performed on clusters with more than 500 cells. We utilized edgeR’s glmQLFTest function to compare each tissue types to Control samples. Significant DEGs were considered at FDR < 0.01 ( Supplementary Table 5 ). Gene Set Enrichment Analysis (GSEA) to GO Biological process (2018) was performed on significant DEG (FDR < 0.00001) with gseapy (version) prerank function for each cell subtype 90 , 91 . Resulting enriched gene ontology list was filtered at FDR < 0.10 ( Supplementary Table 7 ). Analyses were executed with functions implemented in Scanpy (1.7.1) package. Similarities between eutopic endometrium (Control and EuE) tissues were based on hierarchical clustering calculated from Pearson correlation using the Ward linkage algorithm. Cell cycle phase (G1, S or G2M) estimation was calculated following the protocol previously described in Satija et al. 92 and based on markers retrieved from Tirosh et al. 93 . The cell density was estimated with Gaussian kernel density estimation on major cell subtype within each tissue type. Read counts of spliced and unspliced RNA was computed with velocyto (0.17.17) 94 on all 10x libraries obtained from tissue biopsies. We utilized the run10x function which takes output from CellRanger pipeline. Reads were aligned to GRCh38 (10x Genomics reference 3.0.0) and GRCh38 repeat mask downloaded from UCSC Genome Browser as recommended. Projected stream and PAGA trajectory was calculated with scVelo (0.2.3) following recommended workflow previously described 95 , 96 . First, clusters of interest are isolated based on cell barcodes (e.g., myeloid cells in Control). Second, spliced and unspliced counts were log normalized and used for nearest-neighbors estimation. Then, RNA velocity was computed using scVelo’s dynamical model which infers the splicing trajectory for each gene and allows for differential kinetics across distinct lineages and functional states that may be present in the dataset. In order to summarize the mRNA velocity computations, we performed PAGA 96 (specifically regarding DC populations). The number of cells per cell type population vary, especially when subset within a single sample type. With the aim of increasing the robustness of the resulting PAGA graphs, we have utilized the following bootstrapping procedure to generate the edges for all PAGA networks presented in this manuscript: starting with the cells derived from one sample type, we generated 100 randomly subsampled datasets such that each cell type population contained 50 cells; populations with fewer than 13 cells were discarded and those with between 13 and 50 cells were supersampled to include 50 cells. We then recomputed the velocyto neighbor graph, moments distributions, and kinetics, as well as the resulting velocity and PAGA graphs for each randomized dataset 95 . The PAGA connectivities and transitions graphs were used to construct a distribution of linkages between each pairwise combinations of cells. The mean transition probability across these 100 bootstrap samples for each linkage is what is used to plot the PAGA graph, and this process was repeated independently for all sample types comprised by a given group of cell types. Ligand-receptor analysis was performed using CellPhoneDB (2.1.4) 97 on all 58 subclusters. We modified the protocol by running CellPhoneDB on each 10x library separately to reflect the interactions only within individual tissue sample. As such, we added additional parameters to obtain list of interactions that are (1) p-value < 0.01, (2) detected in at least 50% fraction of each tissue type, (3) is not self-interaction, and (4) is a unique cell-to-cell interactions (number of cell type pair is less than 150 counts, Supplementary Fig. 5). The database of ligand-receptor unique interactions obtained from this analysis is supplied in Supplementary Table 6 . The full list of interactions is available at https://github.com/TheJacksonLaboratory/endometriosis-scrnaseq . Formalin-fixed paraffin-embedded (FFPE) tissues were cut into 5-μm sections, mounted on slides, and stained for hematoxylin and eosin (H/E). The slides were then scanned with a Hamamatsu Nanozoomer slide scanner at 40x magnification for histopathological examination by a pathologist. Cell counting and cell type classification was performed with QuPath (0.3.0) using random forest classifier ( Source Data 1 ). Welch’s T-test was performed to obtain p-value. Immunofluorescence staining was performed on FFPE tissue sections. Slides were incubated for 10 minutes at 55°C in a dry oven, deparaffinized in fresh Histoclear (National Diagnostics, #HS-200), and rehydrated through a series of graded alcohols. Antigen retrieval was performed in a decloaking chamber (BioSB TintoRetriever) for 15 minutes at 95°C in neutral citrate buffer, pH 6.00 (Abcam, #ab93678). Tissue was blocked and permeabilized with 10% donkey serum/0.1% Triton X-100 in PBS for 30 minutes at room temperature, then incubated with primary antibodies MUC5B (1/1000, Novus Biologicals, #NBP1–92151) and E-cadherin (5 μg/ml, R&D Systems, #AF648) overnight. Tissue sections were subsequently incubated with secondary antibody Donkey anti-rabbit Alexa Fluor 647 (Invitrogen, #A-31573) and Donkey anti-Goat Alexa Fluor 488 (Invitrogen, #A-11055) for 1 hour at room temperature. DAPI (1 μg/ml, Sigma, MBD0015) was used to counterstain the nuclei, then mounted with ProLong Diamond (Thermo Fisher, # P36970 ). Images were taken using a Leica SP8 Confocal microscope at 40x magnification using Leica Application Suite X (LAS X) and processed with FIJI 98 . FFPE of each tissues were cut into 5-μm sections and mounted on slides. Slides were incubated for 15 minutes at 55°C in a dry oven, deparaffinized in fresh histoclear, and rehydrated through a series of graded alcohols. Antigen retrieval was performed in a decloaking chamber (BioSB TintoRetriever) for 15 minutes at 95°C in citrate buffer, pH 6.0. After blocking in buffer containing 3% BSA, slides were incubated overnight at 4°C with a cocktail of metal-conjugated IMC-validated primary antibodies and described in Supplementary table 9 . The following day, slides were washed twice in DPBS and counterstained with iridium intercalator (0.25 μmol/L) for 5 minutes at room temperature to visualize the DNA. After a final wash in ddH 2 0, the slides were air-dried for 20 minutes. The slides were then loaded on the Fluidigm Hyperion imaging mass cytometer. Regions of interest were selected using Fluidigm CyTOF Software (7.0) and ablated by the Hyperion. The resulting images were exported as 16-bit “.tiff” files using the Fluidigm MCDViewer software and analyzed using napari-imc (0.6.4) 99 or the open source Histocat++ (3.0.0) toolbox/ Histocat web 100 . Fluorescence Activated Cell Sorting (FACS) was performed to isolate MUC5B + and MUC5B − epithelial cells from endometrium dissociated tissue ( Extended Data Figure 10 ). Briefly, control eutopic tissue was dissociated as described above. Cells were then stained with PI and with antibodies marking immune cells (CD45 + ), endothelial cells (CD31 + ), epithelial cells (EpCAM + ), and MUC5B ( Supplementary Table 10 ). We sorted both MUC5B + and MUC5B − epithelial cells (BD Bioscience Symphony S6), gated using FACS Diva (9.0.1), and plated 2000 cells of each population in Matrigel domes (Corning, #356231). Growth of organoid was monitored every 4 hours using an Incuycte S5 (Sartorius) live microscope on brightfield imaging for 10 days. Organoid area and counts were analyzed directly in the onboard Incucyte software ( Source Data ED9 ) and a paired t-test was performed with GraphPad PRISM8 for each timepoint. All hypothesis tests were conducted with the Wilcoxon rank-sum test unless otherwise stated, and the Benjamini-Hochberg correction was used to correct for multiple simultaneous hypotheses tests where applicable. RNA–seq data that support the findings of this study have been deposited in the Gene Expression Omnibus (GEO) under accession codes GSE179640 . In order to further safeguard patients’ genomic identities, SNVs relative to the reference genome are masked in all bam files (BAMboozle v0.5.0) 101 . Moreover, we have made the final single cell datasets available for download and interactive exploration at https://singlecell.jax.org/datasets/endometriosis-2022 . For mapping of scRNA-seq and bulk RNA-seq data, GRCh38.p13 (Ensembl Release 93, https://www.ncbi.nlm.nih.gov/assembly/GCF_000001405.27 ) was used. All other data supporting the findings of this study are available from the corresponding author on reasonable request. Source data are provided with this paper. All code developed for and utilized in this study are available at https://github.com/TheJacksonLaboratory/endometriosis-scrnaseq (temporary link), including modified CellPhoneDB scripts developed to optimize data interpretation for this study.

scRNA-seq was performed on biopsies from 14 individuals. Control eutopic endometrium (Ctrl) represented samples from non-endometriosis patients. Eutopic endometrium (EuE), ectopic peritoneal lesions (EcP) and the adjacent regions to these (EcPA), and ectopic ovarian lesions (EcO) were collected from revised ASRM Stage II-IV patients, the majority under similar hormonal treatment ( Fig. 1a , Extended Data Figure 1a , Supplementary Table 1 – 2 ). EcPA was included in order to study the environment where lesions establish and evolve. In total, 108,497 cellular transcriptomes were generated from tissues, with a median 9,186 unique transcripts and 2,823 genes per cell ( Fig. 1b ). Cells were assigned into five overarching cell types: epithelial, stromal, endothelial, lymphocyte, and myeloid ( Fig. 1c , d , Extended Data Figure 1b ). Subsequent sub-clustering identified 58 subpopulations ( Fig. 1d , Supplementary Table 3 ), highlighting the cellular complexity of both the endometrium and ectopic lesions. We compared bulk transcriptomes from undissociated tissue to pseudo-bulk single cell transcriptomes to interrogate potential biases ( Extended Data Figure 1a ). These showed a strong correspondence within each tissue type ( Extended Data Figure 1c – d ). However, differential expression analysis indicates a few expected cell types are underrepresented in our single-cell dataset ( Extended Data Figure 1e , Supplementary Table 4 ). Nevertheless, similarities between bulk and single-cell transcriptomes indicate that our single-cell dataset reflects much of the original tissue composition and cellular complexity. The heterogeneity among the profiled tissues is evident in cell type composition changes ( Fig 1e , Extended Data Figure 2a ). In order to understand spatial organization and potential cell signaling pathways responsible for changes in cell-type proportions, we designed an antibody panel to spatially resolve cell types of interest with IMC ( Extended Data Figure 2b ). First, we observed endometrium composition in EuE differs dramatically relative to Ctrl, with much of the epithelial component replaced by stroma and lymphocytes in EuE ( Fig. 1e ). Consistent with this, we observed a smaller epithelial proportion in EuE compared to Ctrl ( Fig. 2a – b , Supplementary Table 2 ), together with increased expression of cell cycle-related genes and proliferation of endometrial fibroblasts in EuE ( Fig. 2c – e ). On a per-patient basis, analysis of the scRNA-seq data reveals EuE biopsies stratify into two groups distinguished by immune cell or fibroblast abundance, both distinct from Ctrl samples, and highlights heterogeneity in EuE across patients that exists independent of their hormonal treatment ( Fig. 2f , Extended Data Figure 2a ). Further, osteoglycin ( OGN ) expression is higher in EuE fibroblasts than Ctrl ( Fig. 2g – h , Supplementary Table 5 ), indicating EuE is transcriptomically and compositionally distinct from Ctrl. Second, EcP and EcPA are highly similar, particularly among epithelial cells, which suggests that lesions may extend beyond their macroscopic core and into the surrounding peritoneum ( Fig 1e , f ). Third, the two lesion types display markedly different cellular compositions, confirmed by IMC, where we observed scarce epithelial glands and predominance of stromal cells within EcO when compared to EcP ( Fig 1e – g ). We found endothelial cells (EC) to be markedly increased in EcP, hinting at angiogenesis. We identified cellular components of vasculature: four mural and seven endothelial cell subpopulations, by careful analysis of previously described marker genes 5 , 6 , 9 – 11 ( Fig. 3a – c ). Mural cells, which include vascular smooth muscle (VSMC) and perivascular (Prv) cells, are specialized cells directly interact with ECs providing support and promoting vasculature stabilization. Mural cells account for roughly 40% of stromal cells in EcPA ( Extended Data Figure 3a ) and, together with an increased proportion of EC, suggest a highly vascularized microenvironment ( Fig.1e ). Prv-STEAP4 and Prv-MYH11 were previously identified in the endometrium 6 , and we observed an unreported Prv-CCL19 subpopulation expressing both STEAP4 and MYH11 ( Fig. 3b ). This subpopulation accounts for the majority of perivascular cells in EcP and EcPA ( Extended Data Figure 3b ) and exhibits tissue-specific gene expression ( Extended Data Figure 3c , d ). Interestingly, SUSD2 , a marker for endometrial mesenchymal stem cells identified in endometriosis 12 , is specifically co-expressed in EcP and EcPA CCL19 + Prv cells ( Extended Data Figure 3d ). Prv-CCL19 cells are more abundant in and around peritoneal lesions together with an upregulation of CCL19 and other known angiogenesis regulators such as Synuclein-y ( SNCG) 13 and angiopoietins ( ANGPT1 , ANGPTL1 and ANGPT2 ) 14 . Similarly, Prv-CCL19 upregulate expression of ligands implicated in T-cell recruitment 15 ( CCL21 and FGF7 ) ( Fig. 3d – e , Extended Data Figure 3e ). Together, these data show the presence of an endometriosis-specific perivascular subpopulation, likely promoting angiogenesis and immune chemotaxis in peritoneal lesions ( Fig. 2f ). To further elucidate the interactions between ECs and Prv-CCL19, we performed ligand-receptor analysis to identify unique interactions involving Prv-CCL19 cells (See Methods and Extended Data Figure 4a ); our data indicate EC-tip cells likely respond to ANGPT1 produced by perivascular cells, an interaction known to induce tube formation and branching 16 ( Supplementary Table 6 ). In endometriosis, TEK expression is upregulated while expression of TIE1 , the anti-angiogenic receptor for angiopoietins 14 , 17 , is downregulated ( Fig. 3g ). Previous studies have shown that TEK pathway activation leads to EC proliferation and activation of a feedback loop through DLL4-NOTCH signaling to induce tip EC maturation 18 – 20 . We found significantly increased expression of cell cycle gene CCND1 and decreased expression of DLL4 in EuE, EcP, and EcPA relative to Ctrl ( Fig. 3g , Extended Data Figure 4b , Supplementary Table 5 ), suggesting higher proliferative capacity and sprouting of tip ECs. Contrasting this, the cell cycle arrest gene BTG2 is upregulated and DLL4 - NOTCH signaling is maintained in EcO, suggesting inhibition of sprouting tip ECs in ovarian microenvironment. Immune cell trafficking involves extravasation of immune cells from the blood stream, crossing the ECs barrier into interstitial tissue. Extravasation mainly occurs at the capillary and post-capillary venous (PCV) level 21 . We observed proportions of two PCV subpopulations 10 , 11 , activated PCV (EC-aPCV) and EC-PCV cells proportions, are remarkably increased in EcP and adjacent tissue ( Extended Data Figure 4c – d ). Genes which regulate immune cell attachment and monocyte trafficking 21 — PECAM1 , JAM2 , VCAM1 , ICAM1 , CD99 —and genes associated with EC permeability 22 — PLVAP, AQP1, CXCL12 —are among significantly upregulated genes in endometriosis EC-aPCV, while ICAM2 which encodes a tight junction protein responsible for endothelial-to-endothelial cell contact 21 is downregulated ( Fig. 3h , Supplementary Table 5 ). Expression of AQP1 , which is associated with angiogenesis and migration of endothelial cells 23 , is substantially increased in endometriosis EC-PCVs and EC-aPCVs ( Extended Data Figure 4e – f ). Myeloid and lymphocyte cells are abundant both within and surrounding blood vessels in EcP, pointing towards immune trafficking activity at this site ( Fig. 3i ). Together, these data suggest peritoneal lesions possess leaky PCV vasculature. Thus, our data describe an endometriosis-specific perivascular subtype and emphasize a dynamic orchestration of Prv-CCL19 and PCV endothelial subpopulations, likely promoting angiogenesis and immune cell trafficking in peritoneal endometriosis lesions. This analysis also highlights substantial differences between ovarian and peritoneal lesion microenvironments. scRNA-seq uncovered 15 myeloid cell subpopulations ( Fig 4a , Extended Data Figure 5a ) and 14 lymphocyte subpopulations ( Fig. 1d ). Myeloid cells, particularly macrophages (Mɸs), have been characterized as central components of the endometriosis ecosystem, playing a key role in the establishment of endometriosis 2 . This, together with our observations indicating an increased abundance of myeloid cells in peritoneal lesions ( Fig. 1e ), prompted us to investigate macrophage heterogeneity. We identified five Mɸ subpopulations, of which Mɸ1-LYVE1 and Mɸ3-APOE were previously identified by single-cell analysis in other systems 9 , 24 , 25 . Tissue-resident macrophage subpopulations (Mɸ1-LYVE1, Mɸ2-peritoneal and Mɸ3-APOE) are distinguished by expression of FOLR2— a gene associated with embryonic-derived tissue-resident macrophages 22 , 26 . Mɸ2-peritoneal cells are exclusive to peritoneal tissue and express ICAM2 , a known marker for peritoneal macrophages 27 . Mɸ4-infiltrated cells are present in all tissues, express CLEC5A , CCR2 , and VEGFA , all markers for blood infiltrated macrophages 28 , 29 , and present similar to monocytes ( Fig. 4b , c , Extended Data Figure 5b ). RNA velocity trajectory analyses suggest Mɸ3-APOE are more similar to tissue resident Mɸ1-LYVE1 and Mɸ2-peritoneal subpopulations ( Fig. 4d , Extended Data Figure 5b ). Interestingly, Mɸ5-activated cells appear to arise from both infiltrated and tissue-resident macrophages ( Fig. 4b – d ). Together, these data illustrate the presence of distinct tissue-resident and blood infiltrated macrophage populations in endometrial and lesion tissues. The relative abundance of macrophage subpopulations differ dramatically between control and endometriosis; for instance, Mɸ1-LYVE1 and Mɸ5-activated are enriched in endometriosis eutopic endometrium and ectopic tissues, respectively, and most macrophages present in EcO are Mɸ1-LYVE1 ( Fig. 4c , Extended Data Figure 5c ). Across patients, LYVE1 + macrophages are enriched in both eutopic and ectopic tissues compared to Ctrl ( Fig. 4e ). Tissue resident LYVE1 + macrophages have been previously associated with angiogenesis 22 , 24 , arterial stiffness 30 and anti-inflammatory phenotypes 26 . In agreement, we found that endometriosis Mɸ1-LYVE1 upregulate tolerogenic ( VSIG4 , RGS1 , IL10 , EGFL7, LYVE1 ) and angiogenesis-related genes ( THBS1 , HBEGF , PDGFB , PDGFC, IGF1 ) ( Fig. 4f ). Furthermore, inflammation and antigen presenting pathways of Mɸ1-LYVE1 are downregulated in endometriosis tissues while angiogenesis pathways are enriched in EcPA ( Supplementary Table 7 ). Mɸ1-LYVE1 localization along the vasculature—but not within—confirms the likely link between angiogenesis and this specific cell population ( Fig. 4g ). In addition, two genes ( IGF1 , EMB ) previously shown to promote neurogenesis sprouting in endometriosis 31 and neuromuscular junctions 32 were among the top upregulated genes in Mɸ1-LYVE1 in eutopic and ectopic tissues, suggesting their implication in pain-related mechanisms ( Fig. 4f ). Mɸ4-infiltrated cells, almost completely absent in EcO, presented with pro-tolerogenic features in endometrial tissue, starkly contrasting its pro-inflammatory presentation in control endometrium ( Fig. 4h ). Altogether, we identified endometriosis-associated changes in multiple macrophage subpopulations, promoting tolerogenic, pro-angiogenic, and pro-neurogenic microenvironments. Moreover, the altered macrophage landscape is shared in endometriosis lesions and eutopic endometrium, affecting both tissue resident and blood infiltrated macrophages. Among DCs, we classified three CD1C + populations as pre-cDC2, cDC2, and DC3 according to previously reported markers 25 , 33 , 34 ( Fig 5a ). DC proportions varied greatly across patients ( Extended Data Figure 6a ). However, CD1C + DCs consistently accounted for the majority of the DCs in all tissues ( Fig. 5b , c ). Despite studies reporting altered DC proportions in endometriosis 35 , 36 , the field is still lacking a comprehensive characterization of DC heterogeneity. Previous studies suggest DCs maintain themselves within tissue by proliferating under normal conditions but can be bolstered by an influx of blood-derived DCs during periods of heightened immune activity 37 . Our analysis suggests cDC2s derive from pre-DCs, consistent with the substantial number of proliferative pre-cDC2 cells observed ( Extended Data Figure 6b – c ). The relationship between cDC2s and DC3s appears tissue-specific; both populations seem to derive from an intermediate population in Ctrl (red arrows, Extended Data Figure 6b ) co-expressing FLT3, SIGLEC6, and AXL (DC progenitor- and blood-derived DC-markers 38 ). Meanwhile, cDC2 and DC3 subpopulations appear to derive from pre-cDC2s in EcP ( Extended Data Figure 6d ). As a possible reservoir for tissue-resident DCs, we further interrogated cDC2 diversity. A cDC2 subset in eutopic endometrium and peritoneal lesions, but not EcO, uniquely expresses CD207 , a gene expressed by Langerhans cells and immature DCs 39 ( Fig. 5d – f ). Further analysis revealed two cell states in EcP and EcPA cDC2 populations characterized by mutually exclusive expression of CD207 and MSR1 ( Fig. 5f , g ). Differential gene analysis highlighted CD207 + cDC2 cells express genes related to immunogenic DC maturation ( IL18 , GNLY , RUNX3, LTB ) 40 – 42 , whereas MSR1 + cDC2s express immunomodulatory genes ( MRC1 , VSIG4 , SGK1, and PECAM1 ) 43 – 45 ( Fig. 5f , h ). GSEA analysis of cDC2s across tissues indicated phagocytosis and cytokine-mediated signaling pathways were upregulated in endometriosis tissues ( Extended Data Figure 6e ). These data suggest disease-specific DC heterogeneity and highlights a potential immunomodulatory role for MSR1 + cDC2s in the peritoneum. Next, we interrogated lymphocyte diversity and their interactions with other immune subpopulations ( Fig. 6a and Extended Data Figure 7a , b ). Based on ligand-receptor analysis, we found numerous unique interactions, spatially supported by IMC, between T-cells and various immune subpopulations ( Extended Data Figure 7c ). Particularly, interactions between CD86 (expressed in Mɸ1-LYVE1) and CTLA4 (in T Reg ) is upregulated in endometriosis ( Fig. 6b , Extended Data Figure 7d , Supplementary Table 6 ). This interaction is important for T Reg suppression and homeostasis 46 , and suggests that macrophages’ cooperation with T reg s may be an additional mechanism through which an immunomodulatory microenvironment is promoted in endometriosis. Further, we found that genes associated with T Reg regulatory function are altered between control and endometriosis ( Fig. 6c ): Ctrl T Reg s express HAVCR2 , LAG3, ENTPD1, ICOS, TNFRSF4, and CTLA4 while TIGIT , PRDM1 and CD96 expression is prevalent in endometriosis tissues and EcO specifically. Also noteworthy, ENTPD1—which encodes an important regulator in uterine NK cells promoting immune tolerance and angiogenesis during pregnancy 47 —is upregulated in lesion NK1s ( Fig. 6d ). Collectively, these changes in gene expression indicate modulation of interactions between the various immune subpopulations in ectopic lesions, though the exact mechanism remains unclear. We interrogated the immune cell spatial localization among ectopic lesions using IMC. Unexpectedly, we observed a large clusters of immune cells in 2 out of 7 peritoneal lesions that fit the description of tertiary lymphoid structures (TLSs) ( Fig. 6e – g ). TLSs consist of a germinal center (GC) microarchitecture comprising a central B-cell (CD20 + ) population surrounded by T-cells (CD3 + ) and the additional presence of follicular dendritic cells (PDPN + ) and antigen-presenting cells (HLA-DRA + ), although TLS maturity can modify its composition and organization. TLSs are present in autoimmune disease, chronic inflammatory disease, and tumors but have never been described in endometriosis 48 , 49 . We did not observe similar structures in EcO or across all EcP ( Fig. 6h , i ), suggesting TLS formation may not be a driver of lesions but perhaps a consequence of a sustained inflammatory response. Gene expression analysis of B-cells shows subtle transcriptomic differences among genes related to GC B-cells ( BCL6, SEMA4A and CXCR5 49 , 50 ) ( Extended Data Figure 7e ), suggesting this phenomenon is variable among peritoneal lesions and between patients. Altogether, these data emphasize the diversity among immune cells co-existing within endometriosis lesions. We identified ten epithelial populations comprising the endometrial glands and mesothelium ( Fig. 7a – c ); some have previously been observed in healthy endometrium 5 , 6 whereas other populations differ substantially or have not been previously observed. Namely, we observed previously unreported mesothelial cells found in ectopic tissue and MUC5B+ epithelial cells ( Fig. 7a , Extended Data Figure 8a ). Epithelial cell composition in peritoneal lesions reflects that found in the endometrium while EcO epithelial populations are smaller and less diverse ( Fig. 7b ). This suggests endometrial-like epithelial cells in ovarian and peritoneal lesions may differ in their ability to respond to hormonal or differentiation signals. The MUC5B+ population is present in both eutopic (4–10% of epithelial cells) and ectopic tissues (< 1%) and uniquely expresses RUNX3 , TFF3 and SAA1 ( Fig. 7b and Extended Data Figure 8a ). We confirmed their presence in eutopic endometrium through immunohistochemistry and IMC ( Fig. 7d and Extended Data Figure 8b ). Both trefoil factor 3 (encoded by TFF3 gene) and serum Amyloid A (encoded by SAA1 and SAA2 ) have reported involvement in epithelial restitution, a process initiating mucosal epithelial repair and immune cell recruitment, although details of this mechanism are still unclear 51 – 53 . SAA is a major modulator of inflammation 54 and known promotor of phagocyte chemotaxis through interaction with its receptor FPR2 55 . TFF3 upregulation has been linked to endometriosis and inflammation 56 . This prompted us to look for potential interactions involving MUC5B+ cells through ligand-receptor analysis. We found FPR2 to be uniquely expressed by myeloid cells, particularly monocytes and Mɸ4-infiltrated cells ( Fig. 7e , Supplementary Table 6 ) suggesting a potential interaction between MUC5B+ cells and blood-derived myeloid cells. Additionally, we noted the co-expression of PROM1 and SIX1 - progenitor cell markers 57 , 58 - in in vivo MUC5B+ cells ( Fig. 7b ). To investigate the presence of MUC5B+ cells in in vitro cultures, we established endometrial epithelial organoids (EEO) from primary tissue, starting with unselected single-cell suspensions ( Fig. 7f ). EEOs were maintained in proliferative conditions and subsequently profiled by scRNA-seq, yielding data on 13,326 cells ( Extended Data Figure 8d ). We combined the EEO dataset with the epithelial single-cell transcriptomes from primary tissue and analyzed them jointly ( Fig. 7g ). Of the populations identified in the primary tissue, EEO cells distribute along ciliated (14%) and glandular proliferative epithelial cells (11%), with few differences related to their tissue of origin ( Extended Data Figure 8e ). The largest population of EEO cells (70%) cluster closely between glandular proliferative and MUC5B + cells from primary tissue ( Fig. 7g , arrow) and expresses markers similar to MUC5B + cells in vivo , such as RUNX3, TFF3 and SAA1 ( Extended Data Figure 8f ). The elevated proportion of MUC5B + cells in organoid culture led us to further interrogate the role of MUC5B+ cells. Thus, we isolated MUC5B + and MUC5B − epithelial populations from primary tissue and performed live imaging during organoid derivation ( Extended Data Figure 9a ). Organoids derived from MUC5B + cells grow significantly larger and in higher numbers than those of MUC5B − cells ( Extended Data Figure 9b – c ). Remarkably, MUC5B expression was confirmed in organoids derived from both populations ( Extended Data Figure 9d ). Altogether these data suggest that MUC5B+ epithelial cells may potentially represent a progenitor-like population.

We report a comprehensive description of peritoneal and ovarian endometriosis lesions at single-cell resolution and compare this data to control endometrium, endometrium from endometriosis patients, and organoids derived from these tissues. Our approach was holistic, capturing all cell types (or at least those that survive dissociation) that comprise lesions and their adjacent surroundings and thus providing a view of cellular composition and communication within the niche where lesions establish and evolve. We utilized histopathological imaging and hyperplex antibody-based imaging, with selection of antibodies guided by scRNA-seq, to provide a spatial context. Our data is generated from patients undergoing lesion excision for endometriosis symptom relief and under hormonal treatment, thus representative of the vast majority of patients, as hormonal therapy is the most frequent management strategy for the condition. This continuous low-dose estrogen and/or progestin hormonal treatment comprise the bulk of our control cohort as well. Hormonal treatment induces systemic histological and molecular changes that vary from patient to patient and differ substantially from the cyclic changes observed during the normal menstrual cycle. Our study was specifically designed to interrogate transcriptional and compositional differences between control endometrium and endometriosis, irrespective of hormonal treatment. Indeed, despite the inherent patient-to-patient variability in treatment and clinical history common to human studies, we confirmed robust differences between the two tissues, such as OGN upregulation and increased stromal cell presence. As lesions are described as a piece of endometrial tissue identical to the eutopic endometrium, it was not surprising we identified extensive similarities in cell type composition between the eutopic endometrium and peritoneal lesions. We also detected profound dysregulation of the innate immune and vascular systems in peritoneal lesions. Ovarian lesions, however, display extensive and distinct compositional and gene expression differences relative to peritoneal lesions. Our study provides important clues on the interconnected cellular networks where myeloid, endothelial, epithelial, and perivascular subpopulations influence the formation of the endometriosis-favoring microenvironment ( Extended Data Figure 10 ). Among myeloid cell subpopulations, macrophages and DCs have been described as key players in endometriosis pathology 1 , 2 , 4 with reports showing endometriosis-related alterations of macrophages 29 , 59 and DCs 35 , 36 . However, a comprehensive description of myeloid sub-types was previously lacking. Here, we present a precise characterization of immunomodulatory macrophage and DC populations in peritoneal endometriosis that adopt a coordinated immunotolerant phenotype in the endometriosis microenvironment, where DCs expressing MRC1 and VSIG4 , potentially promoting immunosurveillance escape, and thus benefiting lesion establishment 60 , 61 . Such a phenotype was reported in decidual macrophages and associated with fetal tolerance during pregnancy 62 , 63 . Thus, the present dataset constitutes an ideal starting point to understand how endometriosis may hijack a naturally occurring immunotolerant process to sustain lesion formation and evolution. A deeper understanding of this myeloid compartment in endometriosis is critical, as therapeutics targeting the immune system have been proposed as treatment strategies 2 , 64 , 65 . The spatial analysis provided by IMC provides valuable information to understand the full dynamic of cellular interactions. One such example is the intriguing discovery of TLSs in peritoneal lesions. Their role in endometriosis remains to be determined. A functional understanding of each myeloid and lymphoid subpopulation will determine if TLSs constitute key drivers of the disease, and therefore key therapeutic targets, or are simply a byproduct of the continuous inflammation provoked by lesion settlement. The accumulation of myeloid cells in lesions, together with the presence of LYVE1 -expressing macrophages near the vasculature, is likely linked with increased vascularization, a distinctive trait of peritoneal lesions, and accentuated in the adjacent tissue surrounding these lesions. CCL19 + (and CCL21 + ) perivascular cells may play a role in this as such cells in primary and secondary lymphoid organs have been shown to play a role in immune cell chemoattraction 66 , 67 . This population has not been previously described in endometriosis. Supporting our findings, inhibition of SNCG, a marker uniquely expressed by CCL19 + Prv, prevents endometriosis vascularization and growth 68 , 69 . The peritoneal angiogenic setting contrasts from the ovarian lesion microenvironment where CCL19 + perivascular cells are absent. Thus, while ovarian and peritoneal lesions are currently binned under a common disease name and treatment, we uncover fundamental differences in lesion type that may assist in tailoring lesion-specific therapeutic strategy design such as vascular targeting 70 , 71 for peritoneal lesions. Endometrial epithelial glands form integral components for both eutopic endometrium and endometriotic lesions. The characterization of endometrial epithelial stem cells has been challenging due to the dynamic nature of the regenerative endometrium. Recent single-cell driven descriptions of endometrial epithelial cells from healthy endometrium provide important insights into epithelial subpopulations and the associated hormone responses across the menstrual cycle 5 , 6 , 72 . The field, however, is still lacking a precise characterization of stem or progenitor epithelial cell populations that could explain epithelial gland establishment and initial lesion formation in ectopic tissues. We describe a previously uncharacterized epithelial cell population expressing MUC5B among other specific markers, present in both eutopic and, to a lesser extent, in ectopic tissues. Our success in capturing these cells may be a combination of an optimized tissue dissociation protocol, the use of a hormonally-treated cohort, and/or the speed at which we process these from surgery. However, it remains unclear how this cell subset contributes to endometrial regeneration or the genesis of lesions. Further functional studies will be key to define the precise role of these MUC5B+ cells. While recent studies have begun to describe endometriosis at single cell resolution 73 , 74 , here we have generated the most comprehensive data set, inclusive of spatial organization, describing the eutopic endometrium and ectopic peritoneal and ovarian endometriosis lesions. This atlas represents a unique tool to understand the key players and the dynamic interplay that constitutes the endometriosis niche in hormonally-treated patients. We believe this dataset will be instrumental for designing effective therapeutic strategies or diagnostic biomarkers to provide some relief to the large group of underserved endometriosis patients.

Endometriosis is an inflammatory gynecologic condition that affects 10% of women of reproductive age worldwide 1 , 2 , with symptoms including pelvic pain and infertility. It is characterized by the presence of endometrium-like tissue outside of the uterine cavity (termed lesions), commonly found within the peritoneal cavity, as superficial peritoneal or ovarian lesions. Despite the first description of endometriosis occurring almost a century ago, the exact etiology and molecular drivers of the disease remain largely unknown. Limited non-invasive diagnostic tools impede early detection, resulting in up to seven years’ delay from onset of symptoms to definitive diagnosis, which currently necessitates invasive surgical biopsies of lesions. Treatment of endometriosis remains similarly challenging, relying on hormonal therapy often in conjunction with surgery. Oral contraceptives aim to reduce symptoms but do not necessarily promote lesion clearance. Even after excision, lesions often recur, and repeated surgery is frequent 3 . Challenges in diagnoses and treatment is, at least in part, due to a poor understanding of the pathophysiology of—and heterogeneity within—endometriosis. The tissue microenvironment, including immune cells, has been highlighted as a critical factor for normal endometrium development and endometriosis progression 4 – 8 . Advancements in single-cell RNA sequencing (scRNA-seq) and organoid culture systems enable interrogation of the dynamic interactions within the endometrial microenvironment and the cellular complexity and heterogeneity present in endometriosis. Recent studies have demonstrated the power of such cutting-edge technologies to characterize the dynamic changes of human endometrium through the menstrual cycle and pregnancy 5 – 8 . In this study, we profiled the transcriptomes of endometrium and endometriotic lesions using scRNA-seq and hyperplexed antibody imaging. Our cohorts included patients undergoing oral contraceptive treatment, the most common medical treatment for endometriosis. Consequently, we sought to understand changes within the endometrium and in endometriotic lesions during treatment. Through profiling eutopic endometrium, peritoneal and ovarian lesions, and patient-derived organoids, we uncover distinct cellular changes in endometriosis endometrium as well as specific subsets of immunomodulatory macrophages, immunotolerant dendritic cells (DCs), and vascular changes unique to endometriosis. Our data highlight an unreported endometriosis-specific perivascular population, the presence of tertiary lymphoid structures in some lesions, and a progenitor-like epithelial cell population which may be critical for deeper understanding of this disease.

a , Experimental workflow. b , UMAP showing distribution of cell based on tissue types, PID, and endometriosis stage, before and after batch correction with Harmony. c , Box plot showing Spearman’s correlation rank (ρ) between bulkRNA-seq and pseudobulk from scRNA-seq in Eutopic (Ctrl & EuE, n = 144 ), Peritoneal (EcP & EcPA, n =90 ), or Ovary (EcO, n =24 ). Each dot represents a sample pair. The box represents the interquartile range with median and minimum/maximum represented by box centerline and whiskers, respectively. d , Scatterplot showing distribution of average gene expression between bulk RNA-seq and scRNA-seq (Spearman ρ). Each dot represents a gene. e , Volcano plots representing DEGs between scRNA-seq pseudo bulk (red) and bulk RNA-seq from undissociated tissue (blue) (edgeR, FDR 3). The genes highlighted are exclusively expressed in bulk RNA-seq and associated with erythrocytes (orange), neuronal projections (green), adipocytes (brown), and muscle cells (purple). Related to Fig 1 . a , Major cell types were determined based on UMAP. The mean distribution for all 5 major cell populations is represented for each tissue type Ctrl, EuE, EcP, EcPA, and EcO (left of the line). Each pie chart represents major cell type proportions for each replicate (right of the line). b , Each antibody was selected according to the cell types identified by the scRNA-seq data analysis. Representative images show single channels for each metal-conjugated antibody in a EuE biopsy. A total of 26 antibodies was used to identify cellular heterogeneity within stromal, endothelial, epithelial, lymphocyte, and myeloid major cell types. Additional antibodies (in “Others”) were used to identify cell proliferation (Ki67), active metabolism (pS6), extracellular matrix (Collagen1), and nuclei (DNA). A complete list of cell subpopulations identified through this panel of markers is listed on Supplementary Table 8b . Related to Fig 1 . a , Bar plot representing the proportion of stromal cell types in control endometrium and endometriosis lesions. Endometrial fibroblasts were found in all lesions. Fibroblast C7 is the predominant fibroblast type in EcO. b , Density plot showing distribution of mural cells for each tissue. Arrows points to Prv-CCL19. c , Heatmap of markers genes for mural cell subtypes. d , Track plot representing gene expression pattern for selected DEG in Prv-CCL19 subpopulations. GGT5 and ABCC9 are pan-markers for this cell subtype. e , Box plot showing the proportion of CCL19 -expressing cells in Prv-CCL19 subpopulation within each tissue type. Each dot represents the percentage of CCL19 + cells in a tissue biopsy (Ctrl n = 3, EuE n = 9, EcP n = 8, EcPA n = 6, EcO n = 4). The box represents the interquartile range with median and minimum/maximum represented by box centerline and whiskers, respectively. Related to Fig. 2 . a , Unique cell-to-cell interaction counts obtained from a modified CellPhoneDB procedure. To recover meaningful interactions, we analyzed ligand-receptor interaction in each sample independently. Unique interactions in each tissue type are counted as follows; each ligand-receptor pair observed in a specific cell type pair is counted as one interaction; this is tabulated for all possible pairwise cell type combinations (up to 58 subpopulations in this study) within a sample (n). The total count (Σ, n_celltype_pairs) represents the commonality of the ligand-receptor interaction of interest. The more common interactions (observed in multiple cell type pairs and in all individual samples) will have higher counts while restricted interactions (observed in specific cell type pairs) will have lower counts. We arbitrarily restricted our analysis to interactions observed fewer than 150 times to narrow the scope of analysis and focus on potentially uncovering unique cell-to-cell interactions. b , Box plot showing the proportion of DLL4 -expressing cells in EC-tip subpopulation within each tissue type. c , Density plot showing distribution of endothelial cells for each tissue. d , EC proportions by sample type. e , AQP1 + cell abundance is substantially increased in peritoneal lesions (EcP and EcPA). f , (top) Proportion of aPCV among ECs across tissue types. (bottom) Swarm plot showing AQP1 expression per cell. Horizontal lines represent the median value. For box plots, each dot represents percentage of DLL4+ cells in EC-tip cluster (b) or AQP1+ cells in EC-aPCV cluster (e), in a tissue biopsy (Ctrl n = 3, EuE n = 9, EcP n = 8, EcPA n = 6, EcO n = 4). The box represents the interquartile range with median and minimum/maximum represented by box centerline and whiskers, respectively. Related to Fig. 3 . a, Heatmap representing marker genes for each myeloid subpopulation. b, Dendrogram showing the hierarchical clustering (Pearson correlation) for the myeloid cell clusters. c, Bar plot showing the representation of each myeloid subtype across tissue types. Related to Fig. 4 . a , Bar plot represents the proportion of DCs among all myeloid cells for each patient (Ctrl n = 3, EuE n = 9, EcP n = 8, EcPA n = 6, EcO n = 4). Patient-to-patient variability was observed in DC proportions within the myeloid population and across different tissue types. The box represents the interquartile range with median and minimum/maximum represented by box centerline and whiskers, respectively. b , PAGA and RNA velocity trajectory analyses suggest that pre-cDC2 differentiate towards cDC2 and DC3 in Ctrl and EuE. Red arrows indicate that some cDC2 and DC3 cells derive from a smaller intermediate cell population. c , Cell cycle analysis for pre-cDC2, cDC2 and DC3 populations. d , Expression of DC progenitor markers FLT3 , AXL , and SIGLEC6 . e , Phagocytosis pathway is enriched in cDC2 subpopulations of peritoneal lesions. Bar plot shows the Normalized Enrichment Score (NES) for the top 10- GSEA pathways in cDC2 cells of EuE and EcPA (FDR < 0.1). Related to Fig. 5 . a, Density plot showing distribution of lymphocyte cells for each tissue. b, Dot plot representing marker genes for each lymphocyte subpopulation, including four natural killer cell (NK) clusters, innate lymphoid cells (ILCs), effector memory T-cells (T EM ), cytotoxic T-lymphocytes (CTL), naïve/central memory T-cells (T N /T CM ), T regulatory cells (T Reg ), CD4- and CD8- tissue resident T cell (CD4-T RM and CD8-T RM , respectively), CD8− mucosal-associated invariant T cells (CD8-MAIT), plasma cells, and B cells. c, Representative IMC images showing the presence and proximity of myeloid cells labelled with CD68 (yellow) with T cells labelled with CD3 (cyan), and T Reg labelled with FOXP3 (magenta) in EcO (n = 5); nuclei are marked with DNA intercalation (blue). Scale bar = 100 μm. d, Proportion bar plot of CTLA4 expressing cells from the total T Reg subpopulation. e, Proportion box plot of BCL6 , SEMA4A , CXCR5 expressing cells from the total B cells within each sample type. For box plots, each dot represents a unique patient (Ctrl n = 3, EuE n = 9, EcP n = 8, EcPA n = 6, EcO n = 4). The box represents the interquartile range with median and minimum/maximum represented by box centerline and whiskers, respectively. Related to Fig. 6 . a , Proportions of epithelial subpopulations per sample type. b, Representative IMC images of MUC5B+ epithelial cells in eutopic endometrium (Ctrl: C07, EuE: E12, E06) from multiple tissues. Epithelial cells are marked with PanCK, EpCAM, E-cadherin (green); MUC5B (magenta); nuclei (white). Scale bar = 100 μm. c, Proportion box plot of SAA1 expressing cells from the total MUC5B+ cells within each sample type. Each dot represents a unique patient (Ctrl n = 3, EuE n=9, EcP n=8, EcPA n=6, EcO n=4). The box represents the interquartile range with median and minimum/maximum represented by box centerline and whiskers, respectively. d , Sequencing metrics from EEO scRNA-seq; UMIs and unique genes counts are shown for Control (C) and endometriosis (E) patients and across tissue type. Undetermined (UD) group represents single cells which could not be assigned due to the lack of multiplexing hashtag but otherwise passed QC. e , Density plot showing distribution of EEO cells derived from Ctrl, EuE, EcP, and EcPA (UD cells were not included). f , UMAP showing the co-expression of MUC5B, SAA1, TFF3, and RUNX3 in the MUC5B+ population comprising in vivo epithelial cells and EEO. Related to Fig. 7 . A, Schematic and FACS sorting gating strategy to isolate MUC5B+ and MUC5B- epithelial cells from eutopic tissue for organoid generation. B, Representative brightfield images showing the progression of organoid generation from sorted single cells at day 2, 6 and 10. MUC5B+ cells formed EEO faster than MUC5B- cells. Each panel shows a whole Matrigel dome and magnified organoids are shown in the inset. Inset scale bar = 100 μm. c, Line graph showing area (top) and number (bottom) of EEO generated from MUC5B+ (dark blue, n=1) and MUC5B- (sky blue, n=1) cells over time. Area and Count of EEO is significantly higher in MUC5B+ compared to MUC5B- (paired t-test, two-tailed p < 0.0001). d, IF staining of EEO generated from MUC5B+ (n=1) and MUC5B- (n=1) sorted cells showing the co-localization of endometrial epithelial (E-Cadherin, in green) and MUC5B+ (magenta) staining. Nuclei were counterstained with DAPI (gray). Scale bar = 100 μm. Related to Fig. 7 . In peritoneal lesion (left), the proportion of myeloid and endothelial is increased, and endometrial-like epithelial population is reduced. CCL19-expressing perivascular cells mediate immune cell recruitment, such as macrophages and T cells, which contributes to the immunomodulatory microenvironment. We observe the presence of MSR1-expressing dendritic cells contributing to immunomodulation. TLS is also observed in some lesions. In addition, Mø1-LYVE1 and perivascular cells contribute to angiogenesis by regulating endothelial tip proliferation. In contrast, ovarian ectopic lesions (right) show a striking increase in the proportion of stromal cell and a reduced endometrial-like-epithelial cell presence. The immunomodulatory microenvironment is mainly driven by Mø1-LYVE1 expressing IL10. In ovary lesions, the regulation of angiogenesis is marked by endothelial cell arrest, resulting in mature vasculature. Created with Biorender.com .