Reconstruction of the Decidual Pathways in Human Endometrial Cells Using Single-Cell RNA-Seq

The fate of the human endometrium is determined during the mid-luteal window of implantation, coinciding with differentiation of endometrial stromal cells (EnSCs) into specialized decidual cells. In response to successful embryo implantation, differentiating EnSCs transform the endometrium into a decidua that maintains the placenta throughout gestation; whereas falling progesterone levels in the absence of pregnancy lead to tissue destruction and menstrual shedding. We used single-cell RNA sequencing (scRNA-seq) to map the temporal transcriptomic changes in cultured EnSCs along a decidual time-course and in response to withdrawal of differentiation signals. We demonstrate that decidual transformation of EnSCs first involves a continuous and largely synchronous transition through intermediate states before emerging as divergent subpopulations, representing mature decidual cells and stressed cells. Subsequent withdrawal of decidual signals imposes a second branching event, driven primarily by de-differentiation of a subset of cells. Further, scRNA-seq analysis of timed biopsies highlighted endometrial cellular complexity, confirmed the preponderance of different EnSC subpopulations upon progression from mid-to late-luteal phase, and uncovered genes with conserved branching dynamics in vivo and in vitro. Taken together, our single-cell analysis indicates that multiple EnSC states and dynamic subpopulations of decidual cells underpin endometrial fate decisions during the window of implantation.


Introduction
The human endometrium undergoes iterative cycles of menstrual shedding, regeneration, rapid growth and differentiation in response to the rise and fall in circulating estradiol and progesterone levels (Gellersen and Brosens, 2014). During the mid-luteal phase of the cycle the endometrium begins to remodel intensely, heralding the start of a transient implantation window. This remodeling process is characterized by synchronized transcriptional reprogramming of glandular and luminal epithelium, angiogenesis, influx of uterine natural killer (uNK) and other immune cells, and differentiation of endometrial stromal cells (EnSCs) into secretory decidual cells (Gellersen and Brosens, 2014). Whether or not implantation takes place determines the subsequent fate of the endometrium. In the absence of pregnancy, falling ovarian progesterone levels trigger an inflammatory tissue response resulting in neutrophil infiltration, breakdown of the extracellular matrix, and menstrual shedding (Evans andSalamonsen, 2012, 2014;Jabbour et al., 2006). However, upon implantation, decidualizing EnSCs first encapsulate the conceptus (Weimar et al., 2012;Weimar et al., 2013), and then form a tightly adherent matrix that controls trophoblast invasion and transforms the cycling endometrium into a semi-permanent tissue, which is maintained throughout gestation.
At a molecular level, decidual transformation of EnSCs involves genome-wide remodeling of the chromatin landscape (Vrljicak et al., 2018), wholesale reprogramming of several signaling pathways (Cloke et al., 2008;Leitao et al., 2010), and activation of decidual gene networks (Lynch et al., 2015;Takano et al., 2007). In culture, decidual reprogramming of EnSCs is a multistep process that starts with a burst of free radical production and secretion of inflammatory mediators in response to cyclic adenosine monophosphate (cAMP) and progesterone signaling (Al-Sabbagh et al., 2011;Brighton et al., 2017;Salker et al., 2012). Unlike the secretome produced by undifferentiated EnSCs, decidual inflammation supports blastocyst development in culture and promotes implantation in a mouse model (Peter Durairaj et al., 2017;Salker et al., 2012). After a lag period of several days, morphologically distinct decidual cells emerge, characterized by nuclear enlargement, abundant cytoplasm and expanded endoplasmic reticulum, reflecting acquisition of a mature secretory phenotype. A hallmark of decidual cells is their resistance to a variety of stress signals that trigger an inflammatory response in undifferentiated EnSCs. Several mechanisms underpin stress resistance of decidual cells, including silencing of c-Jun N-terminal kinase (JNK) and inositol trisphosphate (IP3) signaling (Leitao et al., 2010;Muter et al., 2016), downregulation of O-linked N-acetylglucosamine transferase (OGT), a key enzyme in metabolic stress responses (Muter et al., 2018), and induction of free radical scavengers (Kajihara et al., 2006). Decidualizing EnSCs also highly express HSD11B1, coding 11β-hydroxysteroid dehydrogenase type 1. This enzyme catalyzes the conversion of inert cortisone to active cortisol, thus increasing the local levels of anti-inflammatory glucocorticoids at the implantation site (Kuroda et al., 2013). Thus, compared to their fibroblastic progenitor cells, decidual cells are quasi-autonomous, anti-inflammatory and adapted to withstand the intense tissue inflammation caused by deep trophoblast invasion and formation of a hemochorial placenta.
We demonstrated recently that not all EnSCs differentiate into mature decidual cells (Brighton et al., 2017). Instead, a subpopulation of decidualizing EnSCs acquire a senescent phenotype. This subpopulation of stressed decidual cells is then cleared by activated uNK cells through granule exocytosis. We argued that acutely stressed decidual cells may play a pivotal role in kick-starting endometrial remodeling during the implantation window and that their subsequent elimination by uNK cells rejuvenates the endometrium and facilitates its transformation into a gestational tissue (Brighton et al., 2017). In the current study, we employed high-throughput singlecell droplet barcoding to profile individual EnSC transcriptomes across the decidual pathway in vitro, to characterize the emerging subpopulations, and to define cellular responses to withdrawal of the differentiation signal. We then extended our scRNA-seq analysis to timed biopsies and demonstrate that transition of the endometrium from the mid-luteal receptive phase to the lateluteal refractory phase intersects with branching of EnSCs into different subpopulations.

Single-cell analysis of the decidual pathway in vitro
To construct a temporal transcriptomic map along the decidual pathway, confluent primary EnSCs (passage 2) were treated with 8-bromo-cAMP and MPA (medroxyprogesterone acetate, a progestin) for 8 days followed by withdrawal of the deciduogenic signals for an additional 2 days.
Cells were recovered every 48 h and subjected to single-cell analysis using nanoliter droplet barcoding and high throughput sequencing of RNA (Drop-Seq) (Macosko et al., 2015).
Approximately 800 cells were captured and sequenced per time-point ( Figure S1), with an average of 2,267 unique transcripts per cell, corresponding to an average of 1,282 genes per cell ( Figure   S2). After initial quality filtering, a total of 4,580 single cells were computationally assigned to 7 transcriptional cell states using a combination of Shared Nearest Neighbour (SNN) and t-Distributed Stochastic Neighbor Embedding (t-SNE) methods. Figure 1A and 1B show t-SNE plots color-coded according to the day of treatment (D0-8) or withdrawal (WD) and transcriptional states (i-vii), respectively. The expression pattern of the top 10 differential genes in each of the 7 cellular states is presented as a heatmap ( Figure 1C and Table S1).
The time-course analysis yielded a number of unexpected observations. First, we identified a discrete population of cells (state i) that highly express genes involved in cell cycle progression. This population, confined to undifferentiated (D0) and D2 cultures, were all in G2/M phase of the cell cycle ( Figure S3A) and abundantly expressed ANLN ( Figure 1C), coding the evolutionarily conserved actin binding protein anillin that serves as a key mediator of cytokinesis (Piekny and Maddox, 2010). Immunofluorescence microscopy confirmed the presence of anillin-positive cells in independent EnSC cultures maintained over 4 passages ( Figure S3B). A total of 896 genes were differentially expressed between undifferentiated (state ii) and D2 cells (state iii) (Table S2), representing the most profound transcriptional change across the decidual pathway, in keeping with the conjecture that decidual initiation is akin to an acute stress response (Al-Sabbagh et al., 2011;Salker et al., 2012). Subsequent progression to state iv (D4) involved 289 differentially expressed genes (Table S3). Although heterogeneity in cellular states was observed on D2 and D4, the overall pattern infers a continuum of differentiation-associated gene expression during this initial phase of the decidual process. After 4 days of treatment, however, the temporal pattern in gene expression was lost and cells diverged into two subpopulations, represented by states v and vi ( Figure 1B). We identified 326 differentially expressed genes between these subpopulations of decidual cells (Table   S4). Manual mining of the data revealed that state v cells highly express genes that confer decidual autonomy (e.g. DUSP1, PLCL1) and resistance to stress signals (e.g. HSD11B1, IL1RL1, GLRX) (Kuroda et al., 2013;Leitao et al., 2010;Muter et al., 2018;Salker et al., 2016); whereas state vi cells are enriched in genes implicated in oxidative apoptosis [e.g. C10orf10 (also known as decidual protein induced by progesterone, DEPP), FOXO3, NOX4] and cellular senescence (e.g. CLU, ABI3BP) ( Figure S4) (Guo and Chen, 2015;Kajihara et al., 2006;Latini et al., 2011;Salcher et al., 2017;Trougakos, 2013). The emergence of two transcriptionally distinct cell states after 4 days of differentiation substantiates our recent report demonstrating that primary EnSCs give rise to both mature and senescent decidual subpopulations (Brighton et al., 2017). Withdrawal of the deciduogenic signal for 48 h resulted in the emergence of state vii cells, although the response was heterogeneous with multiple cells at this time-point assigned to state ii or vi ( Figure 1B). Compared to either mature (state v) or stressed (state vi) decidual cells, 341 and 361 genes, respectively, were differentially expressed upon withdrawal of 8-bromo-cAMP and MPA (Tables S5 and S6). Finally, two additional primary cultures from different women were analyzed by scRNA-seq across a partial time-course (D0, D2 and D8) to test the robustness of our findings. Aggregation of data from three biological repeat experiments demonstrated that cells cluster on D0 and D2 by state, rather than origin of the culture, and confirmed branching of the decidual trajectory into transcriptionally distinct subpopulations by D8 ( Figure S5). By contrast, expression and secretion of IL-8 was incongruent; with some undifferentiated EnSCs (state ii) expressing but not secreting IL-8 whereas the transient rise in IL-8 secretion on D2 (state iii) coincided with loss of IL8 transcripts ( Figure 2). This pattern, which also applied to IL-6 ( Figure S6), is in keeping with previous observations demonstrating that active release of stored pro-inflammatory mediators plays a role in kick-starting the decidual process (Brighton et al., 2017). Interestingly, a second rise in IL-8 secretion occurred at later time-points, coinciding with the emergence of the two decidual subpopulations, and peaked in response to withdrawal (state vii) despite a concurrent fall in mRNA levels. Discordance between expression and secretion was also apparent for IL-15, a pivotal cytokine for uNK cell activation (Brighton et al., 2017). IL15 mRNA expression increased markedly during the initial decidual phase but secretion of this cytokine also coincided with emergence of distinct decidual subpopulations after D4. Further, IL15 transcripts were more abundantly expressed in mature (state v) decidual cells when compared to stressed (state vi) cells ( Figure 2). Taken together, the data show that decidual factors are secreted in a dynamic but temporally restricted manner. Further, the kinetics of these secretory responses are not necessarily aligned to transcriptional activation of corresponding genes.

Co-regulated gene networks
To provide further insight into the genes that drive decidualization (state ii to vi), 1749 genes with variable expression were used for K-means cluster analysis (Table S7). This analysis yielded 7 coexpressed gene networks, which we ordered in 3 main categories ( Figure 3A-C). Category A encompassed 3 gene networks, totaling 867 genes, that were down-regulated during the initial decidual phase (states iii & iv) and either remained repressed or were re-expressed upon the emergence of the two decidual subpopulations (state v and vi cells). Gene ontology (GO) analysis revealed that these attenuated networks were significantly enriched in genes involved in translation, mRNA splicing and cell-cell adhesion ( Figure 3D and Table S8). Category B, totaling 460 genes, denotes gene networks most highly expressed during either the initial decidual phase (states iii and iv), whereas category C networks, totaling 421 genes, are prominent in mature (state v) and stressed (state vi) decidual cells. Notable GO terms enriched in Category B networks included protein folding and cell redox homeostasis. Further, the most notable GO terms in Category C were signal peptide/secretion and type I interferon signaling ( Figure 3D).

Cell-fate transitions in decidualizing EnSCs
While aggregation analyses categorize heterogeneous cells in different states or populations, they do so by imposing a discrete framework. Pseudotime measures changes in cell state as a function of progress along a continuous trajectory. We used this approach to reconstruct a decidual pathway with a tree-like structure, reflecting divergence in cell-fate trajectories. However, as the initial Withdrawal of differentiation signals on D8 of the decidual time-course triggered further cell-fate divergence (branch point 2) and the emergence of cluster 4 and 5 cells ( Figure 4B).
Modeling of gene dynamics at this branching point revealed that cluster 4 cells upregulate numerous genes that are prominently expressed in undifferentiated EnSCs (state ii), indicating dedifferentiation of a subpopulation of decidual cells ( Figure S7).

Secreted markers of decidual subpopulations
Modeling of gene dynamics at branch point 1 in the decidual pathway suggested that both mature and senescent decidual cells contribute to the peri-implantation secretome. We postulated that secretion of clusterin, a stress-induced molecular chaperone molecule, could serve as a selective marker for senescent decidual cells (state vi), whereas IL15 and CXCL14 secretion is likely restricted to the mature subpopulation (state v). We first employed MAGIC (Markov Affinitybased Graph Imputation of Cells) to impute missing values in our data and then established the gene-gene relationship between CLU, IL15, and CXCL14 in decidualizing cells. As shown in Figure 5A, high CLU and high IL15 or CXCL14 expression were mutually exclusive in decidualizing EnSCs. To test our hypothesis further, three independent primary EnSC cultures were first treated with 8-bromo-cAMP and MPA for 6 days and then co-cultured with or without uNK cells isolated from luteal phase endometrial biopsies. uNK cells selectively target senescent decidual cells, perforate their membrane with perforin, and then trigger apoptosis by releasing granzyme into the cytoplasm. Analysis of culture supernatant revealed that uNK cells do not impact on decidual IL15 secretion. As anticipated, secreted clusterin levels started to rise sharply on D6 of 1 1 the time-course; and this response was markedly blunted upon co-culture with uNK cells (P < 0.01). Interestingly, although secreted CXCL14 levels plateaued around D6, uNK cell-mediated clearance of senescent decidual cells further enhanced secretion of this chemokine (P < 0.05).
Taken together, the data suggest that the decidual microenvironment is determined by the relative abundance of, and crosstalk between, decidual subpopulations; and further dynamically modified by uNK cells.

Single-cell analysis of timed endometrial biopsies
Based on the single-cell analysis in culture, we postulated that dynamic progression of EnSCs along an equivalent cell-state trajectory underpins the transition of endometrium from a receptive to a refractory state. To test this hypothesis, freshly isolated cells from seven endometrial biopsies were subjected to scRNA-seq. The biopsies were timed relatively to the pre-ovulatory luteinizing hormone (LH) surge to coincide with the mid-luteal implantation window (LH+8; n=3) or the start of the refractory late-luteal phase (LH+10; n=4). Out of a total of 4,797 cells, 3,304 passed our quality filtering with an average of 2,774 unique transcripts per cell, corresponding to an average of 1,284 genes per cell. The cells were computationally separated into stromal (n=2,373), epithelial (n=386), endothelial (n=159), and immune cells (n=386). The relative underrepresentation of epithelial cells is accounted for by the resistance of the glands to mild enzymatic digestion used to isolate single cells from the stromal fraction of the biopsies. The expression pattern of the top 10 marker genes for each of the four major cell populations in the endometrium is depicted as a heatmap ( Figure S8A). Figures 6A and 6B show t-SNE plots color-coded according to cell population and the relative contribution of different biopsy samples to each population, respectively. As was the case for primary cultures, we identified a discrete, highly proliferative (HP) subpopulation of EnSCs, representing 1.3% of all stromal cells. Next, we performed pseudotime analysis to test our hypothesis that transition of endometrium from a receptive to a refractory state coincides with progression of EnSCs along a cell-state trajectory. In addition to two minor branch points (designated 1 and 2) ( Figure S8B and S8C), pseudotime analysis yielded one major branch point (designated 3) as highlighted in Figure 6C. The top 50 dynamic genes at this major in vivo branch point are depicted in a modified heatmap ( Figure S8D). Strikingly, the timing of the biopsy relative to the LH surge mapped closely to the progress of EnSCs along the cell-fate trajectory in all but one sample. This sample, designated P1, was obtained at LH+10 but aligned in pseudotime closer to LH+8 samples, indicating that the biopsy was likely mistimed.
We speculated that the relationship between genes that drive branching of cellular trajectories to alternative cell fates would, at least partly, be conserved in vivo and in vitro. To test this hypothesis, we selected the top 50 genes that distinguished branch point 1 in the decidual timecourse ( Figure 4C) and calculated the Pearson correlations for each gene pair across stromal cells in the biopsy data ( Figure 6D, left panel). We added those to correlation coefficients in the in vitro data generating a sum of coefficients ranging from 2, for positively correlated, to -2, for negatively correlated gene pairs. The reversed analysis was also carried out by calculating the conservation in the in vitro data of the top 50 gene-gene interactions that drive divergence of EnSCs at the major in vivo branch point ( Figure 6D, right panel). Congruency was defined as the sum of correlation coefficients of >1 or <-1 for positively and negatively co-regulated genes, respectively. Using this criterion, 27% of the top 50 gene pairs in the in vitro data (branch point 1) were congruent with the in vivo analysis whereas 19% of the top 50 gene pairs in the in vivo data (branch point 3) conformed to the in vitro data. For example, gene pairs that mark the divergence of differentiating EnSC into mature (e.g. SCARA5, IL1RL1, and GLRX) and senescent (e.g. ABI3BP, CLU and IGFBP1) subpopulations were maintained in vivo. Other gene-gene interactions were less conserved. For instance, SCARA5 and PGRMC2 are co-regulated in vitro but this association was not maintained in EnSCs in vivo. Hypergeometric testing revealed that genes involved in branching dynamics in vitro (branch point 1) and in vivo (branch point 3) are significantly more co-regulated than expected by chance alone (P < 10 -82 and P < 10 -33 , respectively; Figure S9). Taken together, the data show that multiple genes involved in the segregation of differentiating EnSCs into mature and stressed subpopulations in culture are also implicated in endometrial fate decisions during the window of implantation.

Materials and Methods
Ethical approval and sample collection. Cell aggregation analysis. Analysis of DGE data was performed with Seurat v2 (Satija et al., 2015). To select high quality data for analysis, cells were included when at least 200 genes were detected, while genes were included if they were detected in at least 3 cells. Cells which had more than 4500 genes were excluded from the analysis as were cells with more than 5% mitochondrial gene transcripts to minimize doublets and low-quality (broken or damaged) cells, respectively.
After scaling and normalization of the raw counts in the DGE matrix, cell-cycle regression was applied. For cell aggregation, a set of highly variable genes was first identified, with an average   Expression for each gene is scaled (z-scored) across cells. Yellow, purple and black represents high, low or undetectable expression of a given marker gene, respectively.           Heatmap top 50 genes distinguishing branch-point 3.