Integrating endometrial proteomic and single cell transcriptomic pipelines reveals distinct menstrual cycle and endometriosis-associated molecular profiles

2.1. Study approval

All participants provided informed consent in accordance with a protocol approved by the MassGeneralBrigham Human Research Committee and the Massachusetts Institute of Technology Committee on the Use of Humans as Experimental Subjects (Protocol number IRB-P001994). Endometrial pipelle biopsies were obtained from reproductive age women (N=17, ages 18-45) undergoing laparoscopic surgery for non-malignant gynecologic indications. Study enrollment was limited to pre-menopausal women and excluded patients with an irregular or ambiguous cycle history or a history of hormone use in the 3 months prior to surgery. A complete list of patients can be found in Supplemental Table 1.

2.2. Tissue Processing for scRNA-seq analysis

A single-cell suspension was generated from fresh endometrial Pipelle biopsies using enzymatic digestion and filter separation as previously described (Osteen et al., 1989). Briefly, the tissue was first minced and washed in 1X PBS to remove contaminating red blood cells (RBCs). Tissues were then incubated in a digest solution composed of collagenase type I (Worthington Biochem, Lakewood, NJ) solution containing DNAse I (SKU# 10104159001, MilliporeSigma, Burlington, MA) for 30 minutes in a 37°C water bath with frequent mechanical dissolution to separate the stromal cells from the epithelial glands. Isolated cells and fragments were then centrifuged and the pellet was resuspended in 5mL of Tryp-LE express (Cat# 12604013, Thermofisher, Waltham, MA) and DNAse I in a 37°C water bath for 15 minutes to generate a single cell suspension. Mechanical disruption of the cells was performed, and visual inspection of single cell populations was assessed using a microscope. Single-cell suspensions of endometrial cells were then either immediately processed for scRNA-seq or cryopreserved at -80°C in a cryopreservation (10% DMSO, 90% FBS) solution.

2.3. Tissue Preparation for Proteomics and Running the LC-MS/MS

Before digestion of the biopsy tissue for single cell analysis, a roughly 5 mm slice was isolated for proteomic analysis. Following collection, 50 μl of RapiGest (Product number 186001861; Waters Corporation, Milford, MA) reconstituted in 50 mM NH₄HCO₃ at a concentration of 2mg/ml was used to digest the biopsy sample for 24 hours at 37°C. Next, urea (Cat # 29700, Thermofisher, Waltham, MA) and dithiothreitol (Cat # R086, Thermofisher, Waltham, MA) were added to the sample to bring the mixture of a final concentration of 4 mM and 5 mM, respectively. This was left at 37°C for 24 hours in an incubator with a stir bar spinning at 1400 rpm. Finally, iodoacetamide (CAS 144-48-9, MilliporeSigma, Burlington, MA) was added to the sample, for a final concentration of 15 mM and kept in the dark, at 37°C, for 30 minutes. Samples were frozen at -80°C until all patient samples had been collected and processed.

Samples were lyophilized overnight and then reconstituted in NH₄HCO₃. A standard mixture was created by mixing equal parts from all of the patient samples into one tube. For each proteomics run, one aliquot of this standard mixture was used to normalize the sample readings between runs. Samples were digested with Trypsin Gold, mass spectrometry grade, from Promega (Madison, WI) using a ratio of 1 ug trypsin to 50 ug sample protein as quantified using a Pierce BCA Protein Assay Kit (Cat # 23225, Thermofisher, Waltham, MA). Samples were digested overnight at room temperature using a rotor for mixing and then placed in a speed vacuum to reduce sample volume. Dried down peptide samples were labeled with TMT10plex (Pierce, Rockford, IL, USA). After labeling, the samples were mixed and purified with peptide desalting spin columns (Pierce, Rockford, IL, USA). 1/5 of the sample was used for one LC/MS/MS analysis. The dried peptide mix was reconstituted in a solution of 2 % formic acid (FA) for MS analysis. Peptides were loaded with the autosampler directly onto a 50cm EASY-Spray C18 column (ES803a, Thermo Scientific). Peptides were eluted from the column using a Dionex Ultimate 3000 Nano LC system with a 5 min gradient from 1% buffer B to 10 % buffer B (100 % acetonitrile, 0.1 % formic acid), followed by a 44.8 min gradient to 25%, and a 9.2 min gradient to 36% B, followed by a 0.5 min gradient to 80% B, and held constant for 4.5 min. Finally, the gradient was changed from 80 % buffer B to 99 % buffer A (100% water, 0.1% formic acid) over 0.1 min, and then held constant at 99 % buffer A for 19.9 more minutes.

The application of a 2.2 kV distal voltage electrosprayed the eluting peptides directly into the Thermo Exploris480 mass spectrometer equipped with a FAIMS and an EASY-Spray source (Thermo Scientific). Mass spectrometer-scanning functions and HPLC gradients were controlled by the Xcalibur data system (Thermo Scientific). MS1 scans parameters were 120,000 resolution, scan range m/z 375-1600, AGC at 300%, IT at 50ms. MS2 scan parameters were either at 30,000 or 45,000 resolution, isolation width at 0.7, HCD collision energy at 36%, AGC target at 300% and IT set to auto. Cycle time for MS2 was 2sec for each MS1 scan.

All MS/MS samples were analyzed using Sequest (Thermo Fisher Scientific, San Jose, CA, USA; version IseNode in Proteome Discoverer 2.3.0.523). Sequest was set up to search uniprot_human_reviewed_032120.fasta (version March 21, 2020) assuming the digestion enzyme trypsin. Sequest was searched with a fragment ion mass tolerance of 0.020 Da and a parent ion tolerance of 10.0 PPM. Carbamidomethyl of cysteine and TMT6plex of lysine were specified in Sequest as fixed modifications. Oxidation of methionine and acetyl of the n- terminus were specified in Sequest as variable modifications.

Scaffold Q+ (version Scaffold_5.0.1, Proteome Software Inc., Portland, OR) was used to quantitate Label Based Quantitation peptide and protein identifications. Peptide identifications were accepted if they could be established at greater than 10.0% probability to achieve an FDR less than 1.0% by the Percolator posterior error probability calculation (Käll et al., 2008). Protein identifications were accepted if they could be established at greater than 99.0% probability to achieve an FDR less than 1.0% and contained at least 2 identified peptides. Protein probabilities were assigned by the Protein Prophet algorithm (Nesvizhskii et al., 2003). Proteins that contained similar peptides and could not be differentiated based on MS/MS analysis alone were grouped to satisfy the principles of parsimony. Channels were corrected by the matrix found in Supplemental Table 2 in all samples according to the algorithm described in i-Tracker (Shadforth et al., 2005). Normalization was performed iteratively (across samples and spectra) on intensities, as described in Statistical Analysis of Relative Labeled Mass Spectrometry Data from Complex Samples Using ANOVA (Oberg et al., 2008). Medians were used for averaging. Spectra data were log-transformed, pruned of those matched to multiple proteins and those missing a reference value, and weighted by an adaptive intensity weighting algorithm.

2.4. Linear Mixed Model for Proteomics Data

Linear mixed modeling (LMM) was used on the list of proteins discovered using isobaric tagging with LC/MS-MS. A total of 668 overlapping proteins were found and for each protein (Supplemental Table 3), a null model and an alternative model was built using the R package “lmer”. Protein expression was modeled by assuming that an individual protein response consisted of patient identity (Patient.ID) as a random variable and two vectors of fixed effect variables: DiseaseState - having endometriosis or not - or cyclePhase – the patient falling into either the proliferative or secretory phase. The two separate null models were built to assess the impact of cycle phase or disease state.

The each of the null models (1) were compared against the alternative model (2) using the likelihood ratio test (LRT).

To determine the influence of cycle phase, the model was run once with disease state as the sole null model fixed variable and the alterative model including both cycle phase and disease state as fixed variables and, to determine the influence of disease state, the null model including cycle phase as the fixed variable. Patient identity was assumed to be a variable in all cases. As a way to generate hypothesis to test, significant proteins (P <0.05) were identified using the LRT. When using cycle phase as the null model fixed effect variable, 115 proteins were found to be significant. This was interpreted as those 115 proteins being significantly associated with disease state variance. Conversely, when disease state was the null model fixed effect variable, 43 proteins were identified as significant. Both proteins lists can be found in Supplemental Table 4 and scripts used for analysis are available for download here github.com/lauren-baugh in the LLM-Prot-Data repository. These identified protein lists were then used as inputs to look for significantly altered pathways using gProfiler (Raudvere et al., 2019) (Supplemental Table 5).

2.5. Generation of scRNA-seq data with Seq-Well S³

Seq-Well S³ was performed as described previously.(Gierahn et al., 2017; Goods et al., 2021; Hughes et al., 2020) Briefly, at least 15,000 cells were loaded onto each array preloaded with barcoded mRNA capture beads (ChemGenes). Arrays were washed with protein-free RPMI media, then sealed with plasma-treated polycarbonate membranes. Arrays were incubated at 37°C for 30 minutes to allow membranes to seal, then transferred through a series of buffer exchanges to allow for cell lysis, transcript hybridization, bead washing, and bead recovery from arrays after membrane removal. Reverse transcription was performed with Maxima H Minus Reverse Transcriptase (ThermoFisher) and excess primers were removed using an Exonuclease I digestion (New England Biolabs). Second strand synthesis was performed and whole transcriptome amplification (WTA) by PCR was performed using KAPA Hifi PCR Mastermix (Kapa Biosystems). WTA product was purified using Agencourt Ampure beads (Beckman Coulter) and 3’ digital gene expression (DGE) sequencing libraries were prepared using Nextera XT (Illumina). Libraries were sequenced using custom primers on a Novaseq S4.

2.6. Computational analyses of scRNA-seq data

2.7. Tissue Sectioning, Staining, and Quantification

Generation of scRNA-seq and proteomics data in healthy and endometriosis patients

In order to develop a more complete proteomic and transcriptomic picture of the healthy and diseased uterus, we established a workflow to profile human endometrial samples. Endometrial biopsy samples were collected from patients undergoing hysterectomy. Prior to uterine tissue removal, pipelle biopsies are taken, primarily from the fundus region of the uterus. These biopsies were then digested as described above to prepare samples for scRNA-seq and proteomic (Figure 1A) analyses. Patients without a history or diagnosis of endometriosis were used as control subjects while diseased patients presented with symptoms of endometriosis with an official diagnosis made during or after surgery. In total, we profiled 11 total patients using scRNA-seq (n=6) and proteomics (n=8). For our proteomic data, we first performed principal component analysis across 668 proteins. Displaying this data along the first and second dimensions revealed little global separation on the basis of cycle phase and of disease state of the patient samples (Figure 1B). There is minor separation of patient disease state along principal component (PC) 1, driven by differences in extracellular matrix composition since COL3A1, COL6A1, and RACK1 were in the top 10 of PC1 contributors.

• Download figure
• Open in new tab

Figure 1. Sample collection workflow and study overview.

A. Sample collection workflow. Endometrial biopsies are collected via pipelle from women undergoing uterine surgery. The sample is split, with part used to isolate single cells for scRNA-seq analysis using Seq-Well, and the remainder digested and analyzed for proteomics using LC-MS/MS. B. A PCA plot of the proteomics samples show poor clustering based on disease state with cycle phase as a confounding factor (n = 8 donors, Supplementary Table 1). C. ScRNA-seq data (n = 6 donors, 6,659 cells) shown as a UMAP projection colored by cell identities (left), cycle stage (middle) or disease status (right).

We generated high quality scRNA-seq data on a total of six patients (3 control, 3 endometriosis) and 6,659 single cells, and identified major epithelial, stromal, and immune populations across all six donors (Supplemental Figure 1). Dimensionality reduction of the scRNA-seq data using PCA followed by uniform manifold approximation and projection (UMAP) on dimensions 1 and 2 highlights populations of cells that cluster strongly by both cycle phase (center, Figure 1C) and disease state (left, Figure 1C). Taken together, our data suggest that both cycle and disease state drive global variations in both data types, and suggests that analytical approaches that leverage both these variables are essential for deeper analyses.

Identification of cycle phase and disease associated proteins

Initial unsupervised clustering methods for both the proteomic and scRNA-seq data highlighted the difficulties of having both menstrual cycle phase and disease state as variables of interest. We therefore used linear mixed effect modeling (LMM) in order to tease apart the effects of these two biological variables in our proteomics results. Using a LMM, we attempt to capture the majority of variance in the system using either fixed or random variables and then by creating a null model to capture for comparison we isolate the variance attributed to the variable of interest. This allows us to assess individual factors while avoiding the pitfalls of trying to understand the contribution of each source of system variance. The likelihood ratio test was used between a null model and an alternative model for each protein as described (see Methods). The model was run twice to generate a list of proteins that showed significant association with cycle phase or disease state. Of the 669 proteins identified from the biopsy tissue, 115 were significantly affected by disease state and 43 were significantly affected by cycle phase, with 28 overlapping proteins. Volcano plots of the identified protein lists show more proteins are upregulated in disease state compared to those that are downregulated (or proteins upregulated in control patients) and more proteins are significantly upregulated in the secretory phase as opposed to the proliferative phase (Figure 2A). PCA of the refined set of proteins from the patient proteomic data separates groups based on disease state (Figure 2B, left) and cycle phase (Figure 2B, right). This suggests that the use of LMM was successful in deconvolving the impacts of cycle phase and disease state, revealing a set of proteins associated with each to investigate further. While numerous protein differences among samples occur as a result of the natural menstrual cycle, an independent set of proteins are associated with the disease state, and a subset of proteins is significantly affected by both cycle phase and disease state.

• Download figure
• Open in new tab

Figure 2: Linear mixed effect modeling (LMM) to identify variable proteins across cycle stage or disease status.

A. Volcano plots of the proteins identified in the LMM using either disease state or cycle phase only in the null model. Green points are above the log₂ fold change cut, blue points are above the -log₁₀ of the p-value cut off and red points are above both the - log₁₀ of the p-value and log₂ fold change cut offs. B. LMM is used to identify proteins that are variable either in disease state or with cycle phase. In the LMM for each protein, the null model excludes the variable of interest (disease or cycle phase) and is compared against the full model to identify proteins with a significant impact on the alternative model variable. PCA plots using the proteins identified in the two different models show better separation of both groups. C. A heatmap of functional enrichment results on protein lists identified as significant in disease or cycle phase.

In order to determine which signaling pathways might be influenced by either cycle stage or disease status, we performed pathway analysis on the LMM protein lists using the software gProfiler (Raudvere et al., 2019). We identified hundreds of pathways that were altered across both protein lists (Supplemental Table 4). To better understand the broader implications of these results, we classified these pathways based on each GO terms description: immune response, extracellular vesicles, cellular transport, cellular binding, and cellular metabolism. We show the top ten scoring pathways in each broad category (Figure 2C). These data suggest that pathways significant to cycle phase skew significantly towards cellular transport pathways while disease-related pathways are notably higher in the immune response, extracellular vesicles, cellular binding, and cellular metabolism categories (Figure 2C). Taken together, our analyses suggest that we are able to effectively tease apart differential proteins and associated biological functions that differ as a function of cycle stage or disease status using this LMM approach.

Multiple cell types are involved in expression of proteomics-identified genes

We next sought to leverage our scRNA-seq data set to better understand our proteomics results by determining which cell types were expressing cycle or disease-defined transcripts. We found that cycle-associated genes, including TPM2, TPM1, COL4A2, and TNS1, were predominantly expressed in smooth muscle cells (Figure 3A). We also found that immune related genes, such as HLA-A, were expressed in macrophages and NK cells. For those proteins that were identified as uniquely associated with disease, we found that many different cell types expressed the corresponding genes (Figure 3B). These included expression of IDH2 across several different clusters of fibroblasts and epithelial cells, high expression of SOD2 in epithelial cells, and high expression of LCP1 in macrophages and NK cells. We additionally use a MAST (Finak et al., 2015) modeling approach to identify genes associated with both cycle and disease in epithelial cells and stromal cells since these were the largest cell populations in our dataset (Figure 3C, Supplemental Table 7 and 8). We identified many genes that were significantly differentially expressed across both disease and cycle, including many that were found in our proteomi data. Comparing our disease-associated proteins with the corresponding lists in our scRNA-seq data, we found there were 71 common elements in epithelial cells and 69 in stromal cells. Similarly, for cycle-associated elements, we found 15 common for epithelial cells and 20 common for stromal cells. Overall, this suggest that many different cell types are contributing to disease and cycle associated alterations in the uterus, and by combining both scRNA-seq and proteomics, we can identify high-confidence lists of genes of interest.

• Download figure
• Open in new tab

Figure 3: Expression of cycle-associated or disease-associated genes in scRNA-seq data.

A. Dotplots of genes uniquely associated with cycle stage as identified in Figure 2 as shown across all cell identities. B. Dotplots of genes uniquely associated with disease status as identified in Figure 2 shown across cell identities. C. Scatter plots of MAST modeling coefficients where each dot is a gene and scaled according to the magnitude of the adjusted p value.

In order to look at these lists more wholistically, we plotted the MAST regression coefficients for cycle stage and disease for either epithelial cells or stromal cells (Figure 3C). These plots allow us to identify genes that significantly contribute to disease or cycle and the magnitude of contributions. For example, we see that for epithelial cells, UBC, SOD2, VCAN, and PAEP were found to be significantly positively associated with cycle stage and disease, while COL1A2 and PGR were both negatively associated. We additionally can see that there are fewer genes associated positively with disease and cycle stage in stromal cells, but they include several heat shock protein genes and LGALS1. The long non-coding RNA H19 was found to be associated positively with cycle stage in stromal cells. Taken together, these plots allow us to further identify genes associated with cycle stage and disease in the uterus and serve as a roadmap for deconvolving the complex interplay between genes associated significantly with both.

Enrichment of proteomics-identified pathways across scRNA-seq data

Next, we wanted to determine which cell types in our scRNA-seq data were enriched for top-scoring pathways identified in our proteomics data. To do this, we created module scores for proteomics-identified pathways and performed PCA on the resulting scores (Figure 4A). We found that we could separate cell types along PC1 or PC2 based on each cell type’s high- scoring pathways. We found that macrophages separated from the rest of the cells along PC2, where the top pathways driving this separation were immune in nature, including myeloid leukocyte mediated immunity (Figure 4B), neutrophil activation, and leukocyte mediated immunity. We additionally found that there was a spectrum of epithelial and stromal cells, with the epithelial 4 cluster and PAEP+ epithelial 2 cluster separating along PC1. This was predominantly driven by selanoamino acid metabolism and regulation of protein transport (Figure 4B). This suggests that there might be different functions of these epithelial and stromal subsets, and in light of other published reports (Tan et al., 2021; Wang et al., 2020), suggests that proteomics coupled with scRNA-seq can help to refine these phenotypes. Taken together, our data suggest that this approach can provide useful context for interpreting long lists of pathway enrichment results from proteomics data and also identify subsets of cells, like stromal, epithelial, and macrophages that drive these results.

• Download figure
• Open in new tab

Figure 4. Module scoring identifies cell types that may drive pathway enrichments from proteomics data.

A. PCA plot of module scores calculated for each top scoring pathway from Figure 2D. PCA plot is colored by cell type. B. Violin plots of the top variable loading pathways along PC1 and PC2.

Mining interactions between stromal cells or epithelial cells and macrophages suggests Axl-mediated signaling is altered across disease status

Given that our proteomic and scRNA-seq analyses suggested that macrophages, stromal cells, and epithelial cells have different roles in cycle and disease-related uterine biology, we wanted to determine how cell-cell interactions might influence these processes. To do this, we used NicheNet, (Browaeys et al., 2020) a cell-to-cell ligand-to-receptor interaction inference approach, to identify receptors and ligands involved in specific interactions identified strongly with disease state (see Methods). We performed this on a patient-by-patient basis in order to more accurately identify interactions between macrophages and stromal or epithelial cells in a subset of patients that had enough macrophages present in the dataset (see Methods). We identified several major receptor and ligand interactions that may be altered in the context of disease, health, and as a function of patient. For interactions between macrophages and epithelial cells, we found several receptor- ligand pairs with high interaction potentials that were conserved across all patients, including IL1B-IL1R, and several that were distinct across patients, such as APP interactions with various receptors (Supplemental Figure 2). Interestingly, for macrophage and stromal cell interactions (Figure 5), we found that Axl-GAS6 receptor-ligand pair scored highly across all patients. We previously found in our MAST modeling results that Axl was differentially expressed in both cycle and disease conditions in stromal cells as well (Supplemental Table 8). Axl is part of group of tyrosine kinase transmembrane receptors which also includes MERTK and Tyro3, collectively known as “TAM” receptors. They are known for their roles in efferocytosis and immune checkpoint blockades (Lemke, 2019; Lemke and Rothlin, 2008; Vago et al., 2021). Interestingly, we also see Gas6-Tyro3 binding between macrophages (Gas6) and stromal cells (Tyro3) as important cell to cell interactions, though this ligand receptor pair is less well studied compared to Axl. Taken together, this data suggests that macrophage interactions with stromal cells, as mediated by Axl and Tyro3 signaling may play a key role in uterine cycle and disease.

• Download figure
• Open in new tab

Figure 5. NicheNet modeling results suggest key receptor and ligand pairs that mediate immune-stromal interactions.

For each patient, top predicted ligands and receptors are shown for macrophages interacting with stromal cells.

Axl and Tyro3 mediated signaling in the healthy and diseased uterus

We next wanted to further investigate the role of Axl, Gas6 and Tyro3-mediated signaling in the healthy and diseased uterus. Further analysis of Axl (data not shown), Gas6 and Tyro3 expression in our scRNA-seq data of all patients shows that Tyro3 is downregulated in patients with disease while Axl and Gas6 do not follow this trend (Figure 6A, B). Following up on Tyro3, and one of its primary ligands found to be expressed in macrophages (from Figure 4) Gas6, expression of Gas6 does not seem to have a notable difference between disease states in patients with a measurable macrophage population (data not shown), while Tyro3 expression in stromal cells, and a lesser degree epithelial cells, is downregulated significantly in diseased patients (Figure 6A and B). To confirm these results and better delineate the role of Tyro3, we used immunohistochemistry on slides of patient endometrial biopsy tissue (representative image in Figure 6C) to stain for Tyro3. Tyro3 staining was found throughout all samples in controls, but was rare in endometriosis patient specimens (Figure 6D). We used CD45 to quantify immune cell populations between patient groups, capturing macrophages and other immune cell types. CD45 expression is comparable to values previously reported (Jones et al., 1996; Vallvé- Juanico et al., 2019) and does not differ between populations as shown in Figure 6E and quantified in Figure 6F. This, combined with our single-cell transcriptomic results, suggests that dysregulation of Tyro3 resides in endometrial stromal cells; however, more detailed analysis is required to confirm this finding.

• Download figure
• Open in new tab

Figure 6. Expression of CD45 and Tyro3 in endometrial tissue samples show no difference in immune cell populations, but Tyro3 is downregulated in disease patients.

A. Gas6 and Tyro3 expression across all epithelial cells. P values were derived from MAST modeling (Supplemental Table 7). B. Gas6 and Tyro3 expression across all stromal cells. P values were derived from MAST modeling (Supplemental Table 8). C. Representative images of histological staining for Tyro3, DAPI, and an overlay with brightfield (bottom) in patient endometrial biopsy samples are quantified in D by comparing the area of Tyro3 staining normalized by the entire tissue area. There is a significant decrease in Tyro3 expression in diseased patients. E. Histological staining for CD45 (some cells marked with yellow arrows), quantified in (F), show no significant difference in CD45 expression