Evolutionary transcriptomics implicates new genes and pathways in human pregnancy and adverse pregnancy outcomes

Evolutionary changes in the anatomy and physiology of the female reproductive system underlie the origins and diversification of pregnancy in Eutherian (‘placental’) mammals. This developmental and evolutionary history constrains normal physiological functions and biases the ways in which dysfunction contributes to reproductive trait diseases and adverse pregnancy outcomes. Here, we show that gene expression changes in the human endometrium during pregnancy are associated with the evolution of human-specific traits and pathologies of pregnancy. We found that hundreds of genes gained or lost endometrial expression in the human lineage. Among these are genes that may contribute to human-specific maternal–fetal communication (HTR2B) and maternal–fetal immunotolerance (PDCD1LG2) systems, as well as vascular remodeling and deep placental invasion (CORIN). These data suggest that explicit evolutionary studies of anatomical systems complement traditional methods for characterizing the genetic architecture of disease. We also anticipate our results will advance the emerging synthesis of evolution and medicine (‘evolutionary medicine’) and be a starting point for more sophisticated studies of the maternal–fetal interface. Furthermore, the gene expression changes we identified may contribute to the development of diagnostics and interventions for adverse pregnancy outcomes.


Introduction
Evolutionary changes in the ontogeny of anatomical systems are ultimately responsible for their functional conservation and transformation into new tissue and organ systems (novelties) with new physiological functions that are outside of the range of the ancestral ones (innovations). These same evolutionary and developmental histories limit (constrain) the range of genetic and environmental perturbations those physiological functions can accommodate before leading to dysfunction and disease (i.e., their reaction norms). Evolution of the structures and functions of female reproductive system and extraembryonic fetal membranes, for example, underlie the evolution of pregnancy (Armstrong et al., 2017;Hou et al., 2009;Kin et al., 2015;Lynch et al., 2015;Lynch et al., 2008) and likely adverse pregnancy outcomes such as infertility (Cummins, 1999), recurrent spontaneous Gene expression changes ultimately underlie the evolution of anatomical structures, suggesting that gene expression changes at the maternal-fetal interface underlie these primate-and humanspecific pregnancy traits. Therefore, we used comparative transcriptomics to reconstruct the evolutionary history of gene expression in the pregnant endometrium and identify genes that gained ('recruited genes') or lost endometrial expression in the primate and human lineages. We found genes that evolved to be expressed at the maternal-fetal interface in the human lineage were enriched for immune functions and diseases such as preterm birth and preeclampsia, as well as other pathways not previously implicated in pregnancy. We explored the function of three recruited genes in greater detail, which implicates them in a novel signaling system at the maternal-fetal interface (HTR2B), maternal-fetal immunotolerance (PDCD1LG2), and remodeling of uterine spiral arteries and deep placental invasion (CORIN). These data indicate that explicit evolutionary studies can identify genes and pathways essential for the normal healthy functions of cells, tissues, and organs, and that likely underlie the (dys)function of those tissue and organ systems.

Endometrial gene expression profiling and ancestral transcriptome reconstruction
To identify gene expression gains and losses in the endometrium that are phylogenetically associated with derived pregnancy traits in humans and catarrhine primates, we assembled a collection of Source data 1. Species names (common and binomial), genome annotations used for RNA-Seq analysis, and parity mode.
Source data 2. Gene expression matrix and ancestral reconstruction results.
Source data 3. Binary encoded matrix of gene expression from extant and ancestral reconstructions used to generate Figure 1B. transcriptomes from the pregnant or gravid endometrium of 20 Eutherian mammals, including human (Homo sapiens), baboon (Papio anubis), Rhesus monkey (Macaca mulatta), and Pig-Tailed macaque (Macaca nemestrina), three Marsupials, platypus, three birds, and six lizard species, including species that are both oviparous and viviparous ( Figure 1 and Figure 1-source data 1). The complete dataset includes expression information for 21,750 genes and 33 species, which were collected at different gestational times from early-to midpregnancy, by multiple labs, and sequencing methods. Thus, differences in transcript abundance between samples may reflect biological differences in mRNA abundances between gestational ages or species, differences in sequencing protocols, or other technical factors unrelated to the biology of pregnancy (i.e., batch effects). Therefore, we transformed quantitative gene expression values coded as transcripts per million (TPM) into discrete character states such that genes with TPM ≥2.0 were considered expressed (state = 1), genes with TPM <2.0 were considered not expressed (state = 0), and missing genes coded as unknown (?; Box 1). Consistent with significant noise reduction, multidimensional scaling (MDS) of species based on gene expression levels (TPMs) was essentially random (Figure 1-figure supplement 1A), whereas MDS of the binary encoded dataset grouped species by phylogenetic relatedness (Figure 1-figure supplement 1B). Next, we used the binary encoded dataset to reconstruct ancestral transcriptomes and trace the evolution of gene expression gains (0 → 1) and losses (1 → 0). Ancestral states were inferred with the empirical Bayesian method implemented in IQ-TREE 2 (Minh et al., 2020;Nguyen et al., 2015) using the species phylogeny ( Figure 1A) and the GTR2 + FO + R4 model of character (Soubrier et al., 2012). Interested readers are referred to other publications for more information about ancestral reconstruction methods (Joy et al., 2016;Pauling et al., 1963). Internal branch lengths of the gene expression tree were generally very short while terminal branches were much longer, indicating pronounced species-specific divergence in endometrial gene expression ( Figure 1A). MDS of extant and ancestral transcriptomes ( Figure 1B) generally grouped species by phylogenetic relationships, parity mode, and degree of placental invasiveness. For example, grouping platypus, birds, and reptiles (cluster 1), viviparous mammals with noninvasive placentas such as opossum, wallaby, and horse, pig, and cow (clusters 2 and 3), and Eutherians with placentas such as mouse, rabbit, and armadillo (cluster 4). Human, baboon, Rhesus monkey, and Pig-Tailed macaque formed a distinct group from other Eutherians (cluster 5), indicating that catarrhine primates have an endometrial gene expression profile during pregnancy that is distinct even from other Eutherians ( Figure 1B).

Gain and loss of signaling and immune regulatory genes in humans
We identified 923 genes that gained endometrial expression in the human lineage with Bayesian posterior probabilities (BPPs) ≥0.80 (Figure 2-source data 2; Box 2). These genes are enriched in 54 pathways, 102 biological process Gene Ontology (GO) terms, and 91 disease ontologies at a false discovery rate (FDR) ≤0.10 ( Figure 2). Among enriched pathways were 'GPCRs, Class A Rhodopsinlike', 'Signaling by GPCR', 'Cytokine-cytokine receptor interaction', 'Allograft Rejection', and 'Graftversus-host disease'. The majority of enriched GO terms were related to signaling processes, such as 'cAMP-mediated signaling' and 'serotonin receptor signaling pathway' or to the immune system, such as 'acute inflammatory response' and 'regulation of immune system process'. The majority of enriched disease ontologies were related to the immune system, such as 'Autoimmune Diseases', 'Immune System Disease', 'Inflammation', 'Asthma', 'Rheumatic Diseases', 'Dermatitis', 'Celiac Diseases', and 'Organ Transplantation', as well as 'Pregnancy', 'Pregnancy, First Trimester', 'Infertility', 'Habitual Abortion', 'Chorioamnionitis', 'Pre-Eclampsia', and 'Preterm Birth', consistent with observations that women with systemic autoimmune diseases have an elevated risk of delivering preterm (Kolstad et al., 2020).
Seven hundred and seventy-one genes lost endometrial expression in the lineage with BPP ≥0.80 ( Figure 2-source data 3). These genes were enriched in 48 pathways, 42 biological process GO terms, and 3 disease ontologies at FDR ≤0.10 ( Figure 2). Enriched pathways included 'immune system ', 'pregnancy', 'pregnancy first trimester', 'infertility', 'habitual abortion', 'preeclampsia', and 'preterm birth'. Unlike genes that gained endometrial expression in the human lineage, those that lost endometrial expression were enriched in disease ontologies unrelated to the immune system, but did include 'Preterm Birth', as well as 'Selenocysteine Synthesis' and 'Selenoamino Acid Metabolism', the latter two which have been previously implicated in preterm birth by genome-wide association study (GWAS;Zhang et al., 2017). In stark contrast, genes that gained (+) or lost (−) endometrial expression during pregnancy in the stem lineage of primates (+63/−34) did not include terms related to the immune system or pregnancy. Thus, genes that gained or lost endometrial expression in the human lineage are uniquely related to immune regulatory process, autoimmunity, inflammation, and allograft rejection, signaling processes such as cAMP-mediated and serotonin receptor signaling, and well as adverse pregnancy outcomes.

Box 1. Classification of genes into not/expressed categories.
A challenge with quantitative gene expression metrics such as RNA-Seq data is defining an expression level that corresponds to functionally active (expressed) genes. Previous studies, however, have shown that an empirically informed operational criterion based on transcript abundance distributions reasonably approximate gene expression categories (Hebenstreit et al., 2011;Kin et al., 2015;Wagner et al., 2013;Wagner et al., 2012). Hebenstreit et al., 2011, for example, showed that genes can be separated into two distinct groups based on their expression levels: the majority of genes follow a normal distribution and are associated with active chromatin marks at their promoters and thus are likely actively expressed, whereas the remaining genes form a shoulder to the left of this main distribution and are unlikely to be actively expressed. Similarly, Wagner et al., 2013;Wagner et al., 2012 found that gene expression data could be modeled as a mixture of two distributions corresponding to inactive and actively transcribed genes. Based on this mixture model, they proposed an operational criterion for classifying genes into expressed and nonexpressed sets: genes with transcripts per million (TPM) ≥2-4 are likely to be actively transcribed, while genes with TPM <2 are unlikely to be actively transcribed. Furthermore, Wagner et al., 2013 suggest that the expression cutoff should be chosen depending on the goal of the study. If it is important to reduce false positives (classifying genes as expressed when they are not), then a conservative criterion of TPM ≥4 could be used. In contrast, if it is more important to reduce false-negative gene expression calls (classifying genes as not expressed when they are), then a liberal criterion such as TPM ≥1 could be used. Both Hebenstreit et al., 2011 andWagner et al., 2013 suggest that genes with TPM ~2 are likely to be actively transcribed. We found that gene expression data (Log 2 TPM) from human decidual stromal cells (DSCs) generally followed a normal distribution with a distinct shoulder to the left of the main distribution, which could be modeled as a mixture of two Gaussian distributions with means of TPM ~0.11 and TPM ~19 (Box 1 - figure 1A). An empirical cumulative distribution fit (ECDF) to the Gaussian mixture model suggests that genes with TPM = 0.01-1 have less than 50 % probability of active expressin, whereas genes with TPM ≥2 have greater than 75 % probability of active expression (Box 1- figure 1A). Next, we grouped genes into three categories, TPM = 0, TPM = 0.01-1, and TPM ≥2, and explored the correlation between these categories and histone marks that are associated with active promoters (H3K4me3) and enhancers (H3K27ac), regions of open chromatin (DNaseI-and FAIRE-Seq), and regions of active transcription (RNA polymerase binding to gene). We found that genes with TPM = 0.01-2 and TPM = 0 were nearly indistinguishable with respect to H3K4me3 marked promoters, H3K27ac marked enhancers, regions of open chromatin assessed by FAIRE-Seq (but not DNaseI-Seq), and most importantly, regions of active transcription as assessed by RNA polymerase binding to gene bodies (Box 1- figure 1B). The promoters of genes with TPM ≥2 were also more enriched for binding sites for the progesterone receptor (PGR) and its cofactor GATA2 than genes with TPM <1 (Box 1 -figure 1C). These data suggest that genes with TPM <1 are unlikely to be actively expressed while genes with TPM ≥2 have hallmarks of active expression. Therefore, we used the TPM ≥2.0 cutoff to define a gene as expressed. We note, however, that other cutoffs could be used that either increase or decrease the probability that genes are actively expressed. Human recruited genes predominantly remodeled the transcriptome of endometrial stromal cells The maternal-fetal interface is composed of numerous maternal and fetal cell types including endometrial stromal lineage cells (perivascular, EFSs, and DSCs), uterine natural killer cells (uNKs), decidual macrophage (uMP), dendritic cells (DCs), T helper cells (Th cells), regulatory T cells (Tregs), various innate lymphoid cells (ILCs), and multiple trophoblast cell types (Suryawanshi et al., 2018;Vento-Tormo et al., 2018;Wang et al., 2020). To infer if genes recruited into endometrial expression in the human lineage are enriched in specific cell types, we used a previously published single-cell RNA-Seq (scRNA-Seq) dataset generated from the first trimester human decidua (Vento-Tormo et al., 2018) to identify cell types at the maternal-fetal interface ( Figure 3A; Figure 3-figure supplement 1). Next, we determined the observed fraction of human recruited genes expressed in each cell type compared to the expected fraction and used a two-way Fisher exact test to identify cell types that were significantly enriched in human recruited genes. Remarkably, human recruited genes were enriched in five of six endometrial stromal lineage cells, including perivascular endometrial mesenchymal stem cells (pvEMSCs) and four populations of DSCs, as well as plasmocytes, endothelial cells (ECs), DCs, and extravillus cytotrophoblasts ( Figure 3B). Consistent with these findings, the expression of human recruited genes defines distinct cell types at the maternal-fetal interface ( Figure 3-figure supplement 2). Our observation that human recruited genes have predominantly remodeled the transcriptome of endometrial stromal lineage cells prompted us to explore the development and gene expression evolution of these cell types in greater detail. Pseudotime single-cell trajectory analysis of endometrial stromal lineage cells identified six distinct populations corresponding to a perivascular mesenchymal stem cell like endometrial stromal population (pvEMSC) population and five populations of DSCs (DSC1-5), as well as cells between pvEMSCs and DSCs that likely represent nondecidualized ESFs and ESFs that have initiated decidualization ( Figure 4A and B). In addition, ESFs that decidualize branch into two distinct lineages, which we term lineage 1 DSCs (DSC1-DSC3) and lineage 2 DSCs (DSC4 Box 2. Genomic features of human recruited genes. The expression of human recruited genes is enriched in 25 tissues at false discovery rate (FDR) <0.05 (Box 1 -figure 1 inset), suggesting these genes were predominately recruited into endometrial expression from those tissues. The expression of human recruited genes in human gestation week 9-22 decidua followed a normal distribution that could be modeled as a mixture of two Gaussian distributions with means of transcripts per million (TPM) ~2.8 and ~ 10.5 (Box 2- figure 1B), grouping genes into low and high expression sets around TPM ~4.2 (Box 2 - figure 1B). An empirical cumulative distribution fit (ECDF) to the gene expression data also suggests a cutoff at TPM ~4.2, which defines an expression level at which 50 % of genes are binned into either the high or low expression sets. The promoters of human recruited genes with TPM <4.2 and ≥4.2 were indistinguishable with respect to H3K4me3 and H3K27ac signal at promoters and enhancers, DNaseI hypersensitive sites, progesterone receptor (PGR), and GATA2 binding, and RNA polymerase binding to gene bodies; both expression sets were generally enriched in these signals compared to genes with TPM = 0 or random genomic locations (Box 2 - figure 1C). In contrast, the promoters of human recruited genes with TPM ≥4.2 are in regions of chromatin with greater nucleosome depletion than recruited genes with TPM <4.2 as assessed by FAIRE-Seq. This observation is consistent with previous studies which found the promoters of highly transcribed genes are preferentially isolated by FAIRE-Seq (Giresi et al., 2007;Nagy et al., 2003). A particularly noteworthy human recruited gene is PRL, which we previously showed evolved endometrial expression in primates (Emera et al., 2012a) and is the most highly expressed human recruited gene in our dataset (Box 2 - figure 1B). Remarkably, PRL gene has independently coopted transposable elements (TEs) into decidual promoters in multiple Eutherian lineages (Emera et al., 2012a;Gerlo et al., 2006;Lynch et al., 2015), suggesting that TE cooption into decidual promoters may be a widespread phenomenon. Consistent with this hypothesis, human recruited genes with TEs in their promoters and 5′-UTRs had greater H3K4me3 signal and nucleosome depletion as assessed by FAIRE-and DNaseI-Seq (Box 2 - figure 1D). The majority of TEs within the promoters and 5′-UTRs of (40%) were primate specific, suggesting they may have played a role in recruiting these genes into endometrial expression (Box 2 - figure 1E).  continued on next page levels as transcripts per million (TPM) for human RefSeq genes (genes with TPM < 2 are classified as "not expressed" and are not shown). Expectation-maximization-based Gaussian mixture curve fits of expression data are shown in purple (low expressed genes) and magenta (high expressed genes), the TPM 4.2 cutoff for defining low and high expressed genes is shown as a black circle. Empirical cumulative distribution fit (ECDF) to the gene expression data is shown in blue, the point of the ECDF at which 50% of genes are binned into either the high or low expression sets is indicated with a blue circle.(C) Correlation of gene expression categories (random genomic locations, TPM = 0, TPM < 4.2, and ≥ 4.2) with histone marks that characterize active promoters (H3K4me3) Figure 2. Enriched pathways, gene ontologies, and disease ontologies among genes that gained or lost endometrial expression in the Hominoid (human) lineage. Data shown as a WordCloud, with term size proportional to −log10 hypergeometric p value (see inset scale) and colored according to enrichment ratio for genes that gained (purple) or lost (green) endometrial expression.
The online version of this article includes the following figure supplement(s) for figure 2: Source data 1. Custom gmt file used for enrichment tests related to preterm birth. and DSC5) ( Figure 4A and B). These cell populations differentially express human recruited genes ( Figure 4C), which are dynamically expressed during differentiation of perivascular cells (PVCs) into lineage 1 and 2 DSCs ( Figure 4D).

Co-option of serotonin signaling in human endometrial cells
Genes that were recruited into endometrial expression in the human lineage are enriched the serotonin signaling pathway (Figure 2), but a role for serotonin signaling in the endometrium has not  previously been reported. Among the recruited genes in this pathway is the serotonin receptor HTR2B.
To explore the history of HTR2B expression in the endometrium in greater detail, we plotted extant and ancestral gene expression probabilities on tetrapod phylogeny and found that it independently evolved endometrial expression at least seven times, including in the human lineage ( Figure 5A). To investigate which cell types express HTR2B, we used the scRNA-Seq dataset from the first trimester maternal-fetal interface and found that HTR2B expression was almost entirely restricted to the DSC cluster ( Figure 5B). We further explored the expression dynamics HTR2B during decidualization using a scRNA-Seq time-course dataset (Lucas et al., 2020) and found that its expression is transiently downregulated during the initial inflammatory decidual phase but upregulated upon the emergence of decidual cells and senescence decidual cells after 4 days of differentiation ( Figure  . Additionally, we found that HTR2B was only expressed in human and mouse ESFs, but not in ESFs at TPM ≥2 from other species in a previously generated multispecies ESF RNA-Seq dataset ( Figure 5C). Pseudotime single-cell trajectory analysis of endometrial stromal lineage cells indicates that HTR2B is expressed in most lineage 1 DSCs, which coexpress other genes such as IL15, INSR, and PRDM1 ( Figure 5D); HTR2B is also expressed by a minority of lineage 2 DSCs, ESFs, and PVCs ( Figure 5D). One hundred and ninety-four genes were differentially expressed between HTR2B + and HTR2B − DSCs ( Figure 5E). These genes were enriched in numerous pathways including 'Regulation of Insulin-like Growth Factor (IGF) transport and uptake by Insulin-like Growth Factor Binding Proteins (IGFBPs)', 'Complement and coagulation cascades', 'BMP2-WNT4-FOXO1 Pathway in Human Primary Endometrial Stromal Cell Differentiation', 'IL-18 signaling pathway', and disease ontologies including 'Small-for-dates baby', 'Premature Birth', 'Inflammation', 'Fetal Growth Retardation', 'Pregnancy Complications', 'Hematologic Complications', and 'Spontaneous abortion' ( Figure 5F and Figure 5source data 1). To determine if HTR2B expression was regulated by progesterone, we used previously published RNA-Seq data from human ESFs and ESFs differentiated into DSCs with cAMP/progesterone . HTR2B was highly expressed in ESFs and downregulated during differentiation (decidualization) by cAMP/progesterone into DSCs ( Figure 6A and Figure 6-figure supplement 1). HTR2B has hallmarks of an expressed gene in DSCs, including residing in a region open chromatin assessed by previously published FAIRE-Seq data ( Figure 6B), an H3K4me3 and H3K27ac marked promoter and polymerase II binding, as well as a promoter that makes long-range loops to binding sites for transcription factors that orchestrate decidualization such as the progesterone receptor A isoform (PGR-A), FOXO1, FOSL2, GATA2, and NR2F2 (COUP-TFII) in previously published ChIP-Seq data (see methods) ( Figure 6B). The HTR2B promoter also makes several long-range interactions to transcription factor-bound sites as assessed by H3K27ac HiChIP data generated from a normal hTERTimmortalized endometrial cell line (E6E7hTERT; see methods) ( Figure 6B). Consistent with regulation by these transcription factors, knockdown of PGR, FOXO1, and GATA2 upregulated HTR2B in DSCs ( Figure 6C). HTR2B is also differentially regulated throughout menstrual cycle ( Figure 6D) and Source data 1. Genes that are differentially expressed between HTR2B + and HTR2B − cells, and the pathways/disease ontologies in which they are enriched.     Figure 6F and

Co-option of PDCD1LG2 (PD-L2) in human endometrial cells
Human recruited genes are enriched numerous immune pathway (Figure 2), among these genes are the PD-1 ligand PDCD1LG2 (PD-L2) ( Figure 7A). We found that PDCD1LG2 was expressed by several cell types at the first trimester maternal-fetal interface, including DCs, macrophages, ESFs and DSCs, and multiple trophoblast lineages ( Figure 7B), and is highly expressed in uterine tissues (Figure 7-figure supplement 1). While each of these cell-type populations has individual cells with high-level PDCD1LG2 expression, only 3%-5% of DSCs, 3 % of DCs, 14 % of macrophage, and 66 % of cytotrophoblasts express PDCD1LG2 ( Figure 7C). Consistent with recent recruitment in the human lineage, PDCD1LG2 was highly expressed in human but either moderately or not expressed in ESFs from other species ( Figure 7D; Figure 5-figure supplement 1). The human PDCD1LG2 locus has the hallmarks of an actively expressed gene, such as a promoter marked by H3K27ac, H3K4me3, and H3K4me1, and binding sites for several transcription factors in previously published ChIP-Seq data from DSCs ( Figure 7E). The PDCD1LG2 promoter also makes several long-range interactions to transcription factor-bound sites, including downstream site that is in the region of open chromatin and bound by PGR/GATA/FOXO1 ( Figure 7E). PDCD1LG2 was highly expressed in ESFs and DSCs ( Figure 7F) but downregulated by cAMP/progesterone treatment ( Figure 7G). Knockdown of PGR and FOXO1 up-and downregulated PDCD1LG2 in DSCs, respectively ( Figure 7G). PDCD1LG2 introns also contain several single nucleotide polymorphisms (SNPs) previously associated with gestational duration and number of lifetime pregnancies as assessed by GWAS (Aschebrook-Kilfoy et al., 2015;Sakabe et al., 2020;Zhang et al., 2017), albeit with marginal p values, implicating PDCD1LG2 in regulating gestation length ( Figure 7E).

Co-option of CORIN into human endometrial cells
Among the human recruited genes enriched in disease ontologies related to preeclampsia (Figure 2) is CORIN ( Figure 8A), a serine protease which promotes uterine spiral artery remodeling and trophoblast invasion (Cui et al., 2012;Yan et al., 2000). We found that CORIN was exclusively expressed by a subset of endometrial stromal lineage cells ( Figure 8B and C), dramatically upregulated in DSCs by cAMP/progesterone treatment ( Figure 8D), and highly expressed in uterine tissues (     ( Figure 8E). The CORIN promoter also makes long-range interactions to transcription factor-bound sites as assessed by HiChIP, including an upstream site bound by PGR, FOSL2, GATA2, FOXO1, and NR2F2 in previously published ChIP-Seq data from DSCs ( Figure 8E). Consistent with these observations, knockdown of PGR, NR2F2, and GATA2 downregulated CORIN expression in DSCs ( Figure 8F).

Discussion
Reconstructing the developmental and evolutionary history of anatomical systems is essential for a causally complete explanation for the origins and progression of disease, which has led to the synthesis of evolution and medicine ('evolutionary medicine') (Benton et al., 2021). We used comparative transcriptomics to explore how the functions of the maternal side (endometrium) of the maternalfetal interface evolved, and found that hundreds of genes gained or lost endometrial expression in the human lineage. These recruited genes are enriched in immune functions, signaling processes and genes associated with adverse pregnancy outcomes such as infertility, recurrent spontaneous abortion, preeclampsia, and preterm birth. Among these genes are those that may contribute to a previously unknown maternal-fetal communication system (HTR2B), augment maternal-fetal immunotolerance (PDCD1LG2 also known as PD-L2), and promote vascular remodeling and deep placental invasion (CORIN).

Human-specific remodeling of the endometrial stromal cell transcriptome
The maternal-fetal interface is composed of numerous maternal cell types, all which could have been equally impacted by genes that were recruited into endometrial expression in the human lineage. It is notable, therefore, that the expression of these genes is predominately enriched in endometrial stromal lineage cells, including perivascular mesenchymal stem cells and multiple populations of DSCs. These data suggest that remodeling of the transcriptome and functions of the endometrial stromal cell lineage has played a particularly important role in the evolution of human-specific pregnancy traits. It is also interesting to note that DSCs evolved in the stem lineage of Eutherian mammals (Carter and Mess, 2017;Gellersen et al., 2007;Gellersen and Brosens, 2003;Kin et al., 2016;Kin et al., 2015;Mess and Carter, 2006), coincident with a wave of gene expression recruitments and losses that also dramatically remodeled their transcriptomes (Kin et al., 2015;Lynch et al., 2015). Thus, the endometrial stromal cell lineage has repeatedly been the target of evolutionary changes related to pregnancy, highlighting the importance of DSCs in the origins and divergence of pregnancy traits. These data also suggest that endometrial stromal lineage cells may play a dominant role in the ontogenesis of adverse pregnancy outcomes.

Co-option of serotonin signaling in human endometrial cells
Unexpectedly, human recruited genes are enriched in the serotonin signaling pathway, such as the serotonin receptor HTR2B. Though a role for serotonin in the endometrium has not previously been reported, we found that serotonin treatment effected RAS/MAPK(ERK) and cAMP/PKA signaling pathways, which are essential for decidualization, and that HTR2B is dynamically expressed during menstrual cycle and pregnancy, reaching a low at term. Previous studies have shown that the human placenta is a source of serotonin throughout gestation (Clark et al., 1980;Kliman et al., 2018;Laurent et al., 2017;Ranzil et al., 2019;Rosenfeld, 2020). Remarkably, a body of early literature suggests serotonin might trigger parturition. For example, levels of both serotonin (5-HT) and 5-hydroxyindoleacetic acid (5-HIAA), the main metabolite of serotonin, are highest in amniotic fluid near term and during labor (Jones and Pycock, 1978;Koren et al., 1961;Loose and Paterson, 1966;Tu and Wong, 1976) while placental monoamine oxidase activity (which metabolizes serotonin) is lowest at term (Koren et al., 1965). Furthermore, a single dose of the monoamine oxidase inhibitor paraglyline hydrochloride can induce abortion in humans and other animals (Koren et al., 1966) Consistent with a potential role in regulating gestation length and parturition, use of selective serotonin reuptake inhibitors is associated with preterm birth (Eke et al., 2016;Grzeskowiak et al., 2012;Huybrechts et al., 2014;Ross et al., 2013;Sujan et al., 2017;Yonkers et al., 2012). 5-HIAA also inhibits RAS/MAPK signaling, potentially by competing with serotonin for binding sites on serotonin receptors (Chen et al., 2011;Klein et al., 2018;Schmid et al., 2015). Collectively, these data suggest a mechanistic connection between serotonin/5-HTAA, and the establishment, maintenance, and cessation of pregnancy.

Co-option of PDCD1LG2 (PD-L2) into human endometrial cells
Among the genes with immune regulatory roles that evolved endometrial expression in the human lineage is the programmed cell death protein 1 (PD-1) ligand PDCD1LG2. PD-1, a member of the immunoglobulin superfamily expressed on T cells and pro-B cells, regulates a critical immune checkpoint that plays an essential role in downregulating immune responses and promoting self-tolerance by suppressing T-cell inflammatory activity (Patsoukis et al., 2020). PD-1 has two ligands, CD274 (PD-L1) and PDCD1LG2 (PD-L2), which upon binding PD-1 promote apoptosis in antigen-specific T cells and inhibit apoptosis in anti-inflammatory Tregs (Patsoukis et al., 2020). Unlike CD274, which is constitutively expressed at low levels in numerous cell types and induced by IFN-gamma, PDCD1LG2 expression is generally restricted to professional antigen-presenting cells (APCs) such as DCs and macrophages and has a fourfold stronger affinity for PD-1 than does CD274 (Ghiotto et al., 2010;Latchman et al., 2001;Sharpe et al., 2007;Sharpe and Pauken, 2018) Remarkably, this higher affinity emerged in the Eutherian stem lineage (Philips et al., 2020). These data suggest that a subpopulation of human DSCs have co-opted some of the immune regulatory functions of professional APCs, which may have been significantly augmented in the human lineage. While more mechanistic studies will help define the precise role of decidual cells in the establishment and maintenance of maternal-fetal immunotolerance, a role for decidual PDCD1LG2 in pregnancy is strongly suggested by its association with variants linked to gestational length and number of lifetime pregnancies (parity) (Aschebrook-Kilfoy et al., 2015;Sakabe et al., 2020;Zhang et al., 2017).

Co-option of CORIN into human endometrial cells
Placental invasiveness varies dramatically in Eutherians, but the cellular and molecular mechanisms responsible for this variation are ill defined. One of the genes that may play a role in the evolution of deeply invasive hemochorial placentation is the serine protease CORIN, which converts pro-atrial natriuretic peptide (pro-ANP) to biologically active ANP (Yan et al., 2000). CORIN-mediated ANP production in the uterus during pregnancy has been shown to promote spiral artery remodeling and trophoblast invasion (Cui et al., 2012). These data implicate co-option of CORIN into endometrial expression may have contributed to the evolution of particularly deep trophoblast invasion and extensive spiral artery remodeling in humans and other great apes (Carter et al., 2015;Pijnenborg et al., 2011a;Pijnenborg et al., 2011b;Soares et al., 2018). CORIN expression is also significantly lower in patients with preeclampsia than in normal pregnancies (Cui et al., 2012), suggesting that the co-option of CORIN into human endometrium may predispose humans to preeclampsia. Additional evolutionary and molecular studies will be required to establish a mechanistic connection between the co-option of CORIN into the endometrium, the evolution of hemochorial placentation, and the origins of preeclampsia in the human lineage.

Caveats and limitations
A limitation of this study is our inability to determine with precise phylogenetic resolution the lineages in which some gene expression changes occurred. For example, we lack pregnant endometrial samples from Hominoids other than humans (chimpanzee/bonobo, gorilla, orangutan, and gibbon/ siamang), thus we are unable to identify truly human-specific gene expression changes. Similarly, we lack endometrial gene expression data from multiple human populations exposed to differing environmental stresses, and therefore are unable to determine the range of physiologically 'normal' gene expression or the reaction norms of individual and collective gene expression levels. Our functional genomic and experimental studies are also restricted to an in vitro cell culture system, which makes it difficult to assess the in vivo impact of gene expression changes on the biology of pregnancy. These limitations are not unique to our study and impact virtually all investigations of Hominoid development and disease, particularly the ones of human-specific traits. Endometrial organoids and iPSCderived ESFs, however, are promising systems in which to study the development of these traits and disease susceptibility that circumvents the limitations of studying human biology (Abbas et al., 2020;Boretto et al., 2017;Marinić et al., 2020;Rawlings et al., 2021;Turco et al., 2017). Our gene expression dataset also represents only a snapshot in time of gestation, rather than a comprehensive time course of endometrial gene expression throughout gestation. Interestingly however, the expression changes we identified from these early time points are enriched in disease ontology terms related to adverse pregnancy outcomes that span the length of gestation including infertility, recurrent spontaneous abortion, preeclampsia, and preterm birth. These findings suggest that atypical gene expression patterns and physiological changes at the earliest stages, perhaps even processes occurring in the endometrium before pregnancy (e.g., decidualization of ESFs into DSCs), may predispose to multiple adverse outcomes, including those at the latter stages like preterm birth (birth before 37 weeks). An important focus of future studies should be collecting endometrial samples across species and from multiple stages of pregnancy, particularly close to term, when the mechanisms that maintain gestation cease and those that initiate parturition are likely to be activated.

Conclusions
We found that hundreds of genes gained or lost endometrial expression in humans, including genes that may contribute to a previously unknown maternal-fetal communication system (HTR2B), enhanced

Endometrial gene expression profiling
Anatomical terms referring to the glandular portion of the female reproductive tract (FRT) specialized for maternal-fetal interactions or shell formation are not standardized. Therefore, we searched the NCBI BioSample, Sequence Read Archive (SRA), and Gene Expression Omnibus (GEO) databases using the search terms 'uterus', 'endometrium', 'decidua', 'oviduct', and 'shell gland' followed by manual curation to identify those datasets that included the region of the FRT specialized for maternal-fetal interaction or shell formation. Datasets that did not indicate whether samples were from pregnant or gravid females were excluded, as were those composed of multiple tissue types. For all RNA-Seq analyses, we used Kallisto (Bray et al., 2016) version 0.42.4 to pseudoalign the raw RNA-Seq reads to reference transcriptomes (see Figure 1-source data 1 for accession numbers and reference genome assemblies) and to generate transcript abundance estimates. We used default mechanisms for maternal-fetal immunotolerance (PDCD1LG2 also known as PD-L2), and deep placental invasion (CORIN). These results demonstrate that gene expression changes at the maternalfetal interface likely underlie human-specific pregnancy traits and adverse pregnancy outcomes. Our work also illustrates the importance of evolutionary studies for investigating human-specific traits and diseases. This 'evolutionary forward genomics' approach complements traditional forward and reverse genetics in model organisms, which may not be relevant in humans, as well as commonly used methods for characterizing the genetic architecture of disease, such as quantitative trait mapping and GWASs. Specifically, our data demonstrate the importance of evolutionary medicine for a mechanistic understanding of endometrial (dys)function, and suggest that similar studies of other tissue and organ systems will help identify genes underlying normal and pathological anatomy and physiology. We anticipate that our results will further the synthesis of evolution and medicine and may contribute to the development of interventions for adverse pregnancy outcomes such as preterm birth.

Materials and methods
parameters bias correction, and 100 bootstrap replicates. Kallisto outputs consist of transcript abundance estimates in TPM, which were used to determine gene expression. To ensure that human decidua samples were free from trophoblast contamination, we compared the expression of placental enriched genes in RNA-Seq data from human placenta, a human ESF cell line, a human decidual stromal (DSC) cell line, and human first trimester decidua. These results suggest that there is likely no trophoblast contamination of human first trimester decidua samples (Table 1), thus inferences of gene expression gains in the human lineage are unlikely to be the result of trophoblast contamination. Next, we compared two different gene expression metrics to reconstruct the evolutionary history of endometrial gene expression: (1) TPM, a quantitative measure of gene expression that reflects the relative molar ratio of each transcript in the transcriptome; and (2) binary encoding, a discrete categorization of gene expression that classifies genes as expressed (state = 1) or not expressed (state = 0). For binary encoding we transformed transcript abundance estimates into discrete character states, such that genes with TPM ≥2.0 were coded as expressed (state = 1), genes with TPM <2.0 were coded as not expressed (state = 0), and genes without data in specific species coded as missing (state = ?); see Box 1 for a detailed justification of the TPM ≥2 cutoff. The TPM coded dataset grouped species randomly (Figure 1-figure supplement 1A), whereas the binary encoded endometrial gene expression dataset generally grouped species by phylogenetic relatedness (Figure 1-figure supplement  1B), suggesting greater signal to noise ratio than raw transcript abundance estimates. Therefore, we used the binary encoded endometrial transcriptome dataset to reconstruct ancestral gene expression states and trace the evolution of endometrial gene expression changes across vertebrate phylogeny ( Figure 1A). Orthology assessment was inferred using Ensembl Compara.

Ancestral transcriptome reconstruction
Ancestral states for each gene were inferred with the empirical Bayesian method implemented in IQ-TREE 2 (Minh et al., 2020;Nguyen et al., 2015) using the species phylogeny shown in Figure 1A and the best-fitting model of character evolution determined by ModelFinder (Kalyaanamoorthy et al., 2017). The best-fitting model was inferred to be the General Time Reversible model for binary data (GTR2), with character state frequencies optimized by maximum likelihood (FO), and a FreeRate model of among site rate heterogeneity with four categories (R4) (Soubrier et al., 2012). We used ancestral transcriptome reconstructions to trace the evolution of gene expression gains (0 → 1) and losses (1 → 0) from the last common ancestor of mammals through to the Hominoid stem-lineage limiting our inferences to reconstructions with BPPs ≥0.80 ( Figure 1A and Figure 1-source data 2). Ancestral reconstructions with BPP ≥0.80 were excluded from over representation analyses.

Data exploration and MDS
We used classical MDS to explore the structure of extant and ancestral transcriptomes. MDS is a multivariate data analysis method that can be used to visualize the similarity/dissimilarity between samples by plotting data points (in this case transcriptomes) onto two-dimensional plots. MDS returns an optimal solution that represents the data in a two-dimensional space, with the number of dimensions (k) specified a priori. Classical MDS preserves the original distance metric, between data points, as well as possible. MDS was performed using the veganR package (Oksanen et al., 2019) with four reduced dimensions. Transcriptomes were grouped using K-means clustering with K = 2-6, K = 5 optimized the number of distinct clusters and cluster memberships (i.e., correctly grouping species by phylogenetic relationship, parity mode, and placenta type).

Reanalyses of Vento-Tormo et al., 2018 endometrial scRNA-Seq data
Maternal-fetus interface 10× Genomics scRNA-Seq data were retrieved from the E-MTAB-6701 entry as a processed data matrix (Vento-Tormo et al., 2018). The RNA counts and major cell-type annotations were used as provided by the original publications. Seurat (v3.1.1) (Butler et al., 2018), implemented in R (v3.6.0), was used for filtering, normalization, and cell types clustering. The subclusters of cell types were annotated based on the known transcriptional markers from the literature survey. Briefly, we performed the following data processing steps: (1) cells were filtered based on the criteria that individual cells must be expressing at least 1000 and not more than 5000 genes with a count ≥1; (2) cells were filtered out if more than 5 % of counts mapping to mitochondrial genes; (3) data normalization was performed by dividing uniquely mapping read counts (defined by Seurat as unique molecular identified [UMI]) for each gene by the total number of counts in each cell and multiplying by 10,000. These normalized values were then log-transformed. Cell types were clustered by using the top 2000 variable genes expressed across all samples. Clustering was performed using the 'FindClusters' function with essentially default parameters, except resolution was set to 0.1. The first 20 PCA dimensions were used in the construction of the shared-nearest neighbor (SNN) graph and the generation of twodimensional embeddings for data visualization using UMAP. Major cell types were assigned based on the original publication samples' annotations, and cell subtypes within major cell types were annotated using the subcluster markers obtained from the above parameters. We then chose the decidual and PV cells to perform the single-cell trajectory, pseudotime analysis, and cell ordering along an artificial temporal continuum analysis using Monocle2 (Qiu et al., 2017). The top 500 differentially expressed genes were used to distinguish between the subclusters of decidua and PV cell populations on pseudotime trajectory. The transcriptome from every single cell represents a pseudo-time point along an artificial time vector that denotes decidual and PV lineages' progression, respectively. To compare the differentially expressed genes between HTR2B-positive and HTR2B-negative cells, we first divided the decidual and PV datasets into those groups of cells that either express HTR2B with a count ≥1 and those with zero counts. We then performed differentially expressed genes analysis between the mentioned two groups of cells using the bimodal test for significance.
To calculate the enrichment score of human-gain genes in each cell type, we first transformed the data into a pseudobulk expression matrix by averaging all genes' expression in each cell type. We then calculated the fraction of human-gained genes expressed (Observed) and the proportion of the rest of the genes expressed in each cell type (Expected). The enrichment ratio shown on the plot is the ratio of Observed and Expected values for each cell type. The p value was calculated using a two-way Fisher exact test followed by Bonferroni correction.

Reanalyses of HTR2B, PDCD1LG2, and CORIN expression across multiple endometrial scRNA-Seq datasets
Transcriptomic dynamics of human endometrium in vivo. Data mined from publicly available database repr oduc tive cell atlas. org (Garcia-Alonso et al., 2021). Data show a cellular map of the human endometrium from combinatorial transcriptomics (scRNA-Seq and single-nuclei RNA sequencing [snRNA-Seq]) alongside spatial transcriptomics methods (10× Genomics Visium slides and highresolution microscopy) representing 98,568 cells from fifteen individuals grouped into five main cellular types. No reuse allowed without permission.
Single-cell analysis of peri-implantation endometrium. Six LH-timed endometrial biopsies were processed for Droplet generation and single-cell sequencing (Drop-Seq) as described in Lucas et al., 2020. Anonymized endometrial biopsies were obtained from women aged between 31 and 42 years with regular cycles, body mass index between 23 and 32 kg/m 2 , and the absence of uterine pathology on transvaginal ultrasound examination. t-Distributed stochastic neighbour embedding (t-SNE) analysis assigned 2831 cells to four clusters, designated based on canonical marker genes as ECs (n = 141), epithelial cells (EpC; n = 395), immune cells (IC; n = 352), and ESFs (n = 1943). Data are available in the GEO repository GSE127918.

Deconvolution of in vitro cell types from endometrium
Three independent midluteal biopsies were used to isolate, culture, and sequence different endometrial cell types. Following enzymatic digestion, EpC were separated from the stromal cell fraction as described (Barros et al., 2016). PVCs and ESFs were then subjected to magnetic activated cell sorting (MACS) using W5C5 antibody (Masuda et al., 2012). PVC and ESFs were maintained in standard cultures as well as subjected to colony-forming unit assays. Total RNA was then extracted from the resulted clones, designated endometrial mesenchymal stem cells (eMSCs) and transit amplifying cells (TA), respectively. The standard PVC and ESF cultures were propagated until 90 % confluence and then subjected to total RNA extractions. Primary EpC were subjected to gland organoid formation (Turco et al., 2017). uNK cells were also isolated although not cultured. After overnight incubation, the supernatant of the stromal cell fraction was subjected to MACS to isolate uNK cells using a PE-conjugated antihuman CD56 monoclonal antibody. Libraries were prepared using TruSeq RNA Library preparation kit V2 and sequenced on HiSeq 4000 with 75 bp paired-end reads. Data are presented as TPM (Diniz-da-Costa et al., 2021 [unpublished thesis] Endometrium through the menstrual cycle Data mined from the publicly available database GDS2052 (Talbi et al., 2006). Data are presented as average gene counts.

Single-cell analysis of the decidual pathway in vitro
Primary ESFs were decidualized with a progestin (medroxyprogesterone acetate, MPA) and a cyclic adenosine monophosphate analog (8-bromo-cAMP, cAMP) for 8 days. Cells were recovered every 48 hr and subjected to single-cell analysis using nanoliter droplet barcoding and high-throughput RNA sequencing. Approximately 800 cells were sequenced per time point, yielding on average 1282 genes per cell. After computational quality control 4580 cells were assigned to 7 transcriptional cell states (6 presented) using SNN and t-SNE methods and presented as transcriptional states. Data are available in the GEO repository GSE127918. Data are presented as TPM (Lucas et al., 2020). We note that Vento-Tormo et al. identified five populations of cells in the endometrial stromal lineage, including two perivascular populations (likely reflecting the mesenchymal stem cell-like progenitor of ESFs and DSCs) and three cell types they call 'decidual stromal cells' and label 'S1-3'. However, based on the gene expression patterns of 'dS1-3' (shown in Vento-Tormo et al. Figure 3a), only 'dS3' are decidualized, as indicated by expression of classical markers of decidualization such and PRL (Tabanelli et al., 1992) and IGFBP1/2/6 (Tabanelli et al., 1992;Kim et al., 1999). In stark contrast, 'dS1' do not express decidualization markers but highly express markers of ESFs such as TAGLN and ID2, as well as markers of proliferating ESFs including ACTA2 (Kim et al., 1999). 'dS2' also express ESFs markers (TAGLN, ID2, ACTA2), but additionally LEFTY2 and IGFBP1/2/6, consistent with ESFs that have initiated the process of decidualization. These data indicate that the 'dS1' and 'dS2' populations are both ESFs, but 'dS2' are ESFs that have initiated decidualization (because they express IGFBPs but not PRL), and that 'dS3' are DSCs. Vento-Tormo et al. show that the differences in gene expression between 'dS1-3' are related to their topography in the endometrium, but degree of decidualization ('dS1'/ESF1 < 'dS2'/ESF2 < 'dS3'/DSC) is also linked to differential gene expression.

Endometrial stromal lineage cell nomenclature
Consistent with this, other scRNA-Seq studies have identified two ESF populations and one DSC population in the first trimester decidua, and used pseudotime analyses to show that they represent different states of differentiation from ESFs to mature DSCs (Suryawanshi et al., 2018). Therefore, we prefer to use the perivascular/ESF/DSC nomenclature because it more accurately reflects the biology and gene expression profile of these cell types than the 'dS1-3' naming convention. We also note that while it is generally thought that ESFs are absent from the pregnant uterus, ESFs retain a presence in the endometrium from the first trimester until term (Richards et al., 1995;Suryawanshi et al., 2018;Muñoz-Fernández et al., 2019;Sakabe et al., 2020).

Overrepresentation analyses
We used WebGestalt v. 2019 (Liao et al., 2019) to identify enriched ontology terms using overrepresentation analysis (ORA). We used ORA to identify enriched terms for three pathway databases (KEGG, Reactome, and Wikipathway), three disease databases (Disgenet, OMIM, and GLAD4U), and a custom database of genes implicated in preterm birth by GWAS. The preterm birth gene set was assembled from the NHGRI-EBI Catalog of published GWASs (GWAS Catalog), including genes implicated in GWAS with either the ontology terms 'Preterm Birth' (EFO_0003917) or 'Spontaneous Preterm Birth' (EFO_0006917), as well as two recent preterm birth and birth weight GWASs (Warrington et al., 2019;Sakabe et al., 2020) using a genome-wide significant p value of 9 × 10 -6 . The custom gmt file used to test for enrichment of preterm birth associated genes is included as a supplementary data file to (Figure 2, Figure 2-source data 1).
We also used previously published RNA-Seq and microarray gene expression data generated from human ESFs and DSCs that were downloaded from NCBI SRA and processed remotely using Galaxy platform (https:// usegalaxy. org/; Version 20.01) for RNA-Seq data and GEO2R for microarray data. RNA-Seq datasets were transferred from SRA to Galaxy using the Download and Extract Reads in FASTA/Q format from NCBI SRA tool (version 2.10.4+ galaxy1). We used HISAT2 (version 2.1.0+ galaxy5) to align reads to the Human hg38 reference genome using single-or paired-end options depending on the dataset and unstranded reads, and report alignments tailored for transcript assemblers including StringTie. Transcripts were assembled and quantified using StringTie (v1.3.6) (Pertea et al., 2016;Pertea et al., 2015), with reference file to guide assembly and the 'reference transcripts only' option, and output count files for differential expression with DESeq2/edgeR/limma-voom. Differentially expressed genes were identified using DESeq2 (Love et al., 2014) (version 2.11.40.6+ galaxy1). The reference file for StringTie guided assembly was wgEncodeGencodeBasicV33. GEO2R performs comparisons on original submitter-supplied processed data tables using the GEOquery and limma R packages from the Bioconductor project (https:// bioconductor. org/). These datasets included gene expression profiles of primary human ESFs treated for 48 hr with control nontargeting, PGRtargeting (GSE94036), FOXO1-targeting (GSE94036), or NR2F2 (COUP-TFII)-targeting (GSE47052) siRNA prior to decidualization stimulus for 72 hr; transfection with GATA2-targeting siRNA was followed immediately by decidualization stimulus (GSE108407). Probes were 206638_at (HTR2B), 220049_s_at (PDCD1LG2), and 220356_at (CORIN) for GSE4888 (endometrial gene expression throughout menstrual cycle) and for GSE5999 (gene expression in basal plate throughout gestation). Multispecies RNA-Seq analysis of ESFs and DSCs is from GSE67659.

Immunofluorescent staining for endometrial HTR2B
Endometrial biopsies were fixed overnight in 10 % neutral buffered formalin at 4 °C and wax embedded in Surgipath Formula 'R' paraffin using the Shandon Excelsior ES Tissue processor (Thermo Fisher). Tissues were sliced into 3 μm sections on a microtome and adhered to coverslips by overnight incubation at 60 °C. Deparaffinization and rehydration were performed through xylene, 100 % isopropanol, 70 % isopropanol, and distilled water incubations. Following antigen retrieval, slides were washed, blocked, and incubated in primary HTR2B antibody (1:200; Fisher Scientific) overnight at 4 °C. After washing three times, slides were incubated with Alexa Fluor 594 (1:1000; Fisher Scientific) for 2 hr, washed and mounted in ProLong Gold → Antifade Reagent with DAPI (Cell Signaling Technology). Slides were visualized using the EVOS Auto system, with imaging parameters maintained throughout image acquisition.
A total of 10 4 ESFs were plated per well of a 96-well plate, 18 hr later cells were transfected in Opti-MEM (31985070, Thermo Fisher Scientific) with 100 ng of luciferase reporter plasmid, 10 ng Renilla control plasmid, 0.25 μl of Lipofectamine LTX (15338100, Thermo Fisher Scientific) and 0.1 μl Plus Reagent as per the manufecturer's protocol per well; final volume per well was 100 μl. Luciferase reporter plasmids were synthesized (GenScript) by cloning the response elements from the pGL4. ESFs were incubated in the transfection mixture for 6 hr. Then, ESFs were washed with warm PBS and incubated in the maintenance medium overnight. The next day, the medium in half of the wells was exchanged for the differentiation medium consisting of DMEM with Phenol Red and GlutaMAX (10566-024, Thermo Fisher Scientific), supplemented with 2 % fetal bovine serum (26140-079, Thermo Fisher Scientific), 1 % sodium pyruvate (11360070, Thermo Fisher Scientific), 1 μM medroxyprogesterone 17-acetate (MPA; M1629, Sigma Aldrich), and 0.5 mM 8-Bromoadenosine 3′,5′-cyclic monophosphate (8-Br-cAMP; B5386, Sigma Aldrich). After 48 hr, serotonin (5-HT; H9523, Sigma Aldrich) was added to the wells with both maintenance and differentiation medium (for each in triplicates) in the following concentrations: 50 μM, 200 μM, and 1 mM; vehicle control (0 μM) was water. After 6 hr of incubation, we used a Dual Luciferase Reporter Assay (Promega) to quantify luciferase and Renilla luminescence following the manufacturer's Dual Luciferase Reporter Assay protocol.

Cell lines
Human hTERT-immortalized endometrial stromal fibroblasts were purchased from ATCC (T-HESC; CRL-4003, ATCC). Their identity has been authenticated by ATCC, and was determined to be mycoplasma free.