Tracing the cis-regulatory changes underlying the endometrial control of placental invasion

The ELI hypothesis states that bovine (and other epitheliochorial) stromal fibroblasts have evolved to be more resistant to trophoblast invasion and that endometrial invasibility also is correlated with susceptibility to cancer malignancy between species (2) (Fig. 1A). In this study, we use ESFs from nine species with invasive and noninvasive placentation (Fig. 1A). We use a linear regression model to estimate the extent to which cis-regulatory TFBS numbers could statistically explain gene expression differences across species (Fig. 1B). Using gene expression quantified as transcripts per million (TPM), of 8,639 genes with 1–1 orthology among our species sample, we first calculated their square root TPM values to decrease the influence of extremely high expression values. We then subtracted the gene-wise mean sqrt(TPM) across all species from the square root of the species specific sqrt(TPM). This scaled expression measure was fit to a linear model with the number of TFBSs in promoter regions as the predictor variable.

Regressing the TFBS number for 572 TFs for which consensus binding site motifs are available in the JASPAR database (17), we identified 187 TFs at a false discovery rate (FDR) less than 0.05, which could explain a significant portion of the interspecies differences in gene expression in endometrial stromal cells (Fig. 1C). A large fraction of TFBSs statistically influence interspecies gene expression variation, indicating that the gene expression profile of endometrial stromal cells has evolved, in part, through changes in cis-regulatory elements (Fig. 1C).

We found that β, the regression coefficient describing the extent to which TFBS abundance could explain variation in gene expression, was positive for most TFs (Fig. 1D). A positive β indicates that the gain of TFBS is correlated with increased expression of the target gene, suggesting that most TFs acted as positive regulators of gene expression (i.e., activating TFs) in ESFs. The range of β values was similar for positive and negative regression coefficients (Fig. 1D). We found that for many TFs, |β| was above 0.025 for the whole transcriptome, with high significance (FDR < 10⁻²⁰) (Fig. 1E and Dataset S1). Notably, averaged β over all genes did not reflect the large effect a single TF may have on the expression of individual genes (Fig. 1 E and F). For example, the regression coefficient |β| for the TF MGA could be 20 times higher for some genes than the average β over all genes, demonstrating specificity for particular target genes (Fig. 1F). The interspecies distribution of the number of binding sites for the same TF for all genes showed that enrichment for specific TFBS was high for a small number of genes among hemochorial species (Fig. 1G). As an example, Fig. 1H shows the expression of the target gene LRCH3, as a function of the number of binding sites for MGA.

Line_NoWe then identified the top TFs with the highest R² values (Fig. 1I). These TFs consisted of various nuclear receptors including NR2F6, a steroid receptor which is also a repressor of IL17, the Estrogen Receptor ESR1, as well as SREBF1, the sterol regulatory element binding TF 1 (18). NR2F6 coregulates and inhibits thyroid receptor (19) and oxytocin receptor (OXTR) (20, 21), while ESR1 directly regulates OXTR expression (22) and is considered essential for uterine receptivity and endometrial decidualization (23, 24). In addition, among the top candidates an important factor with a negative β value was FOXA2, the Forkhead Box A2, an essential TF for fertility and uterine function (25, 26). In addition, other top TFs with high β values were Max Dimerization Protein MGA, a negative regulator of MYC pathway and implicated in endometrial cancers (27, 28), and PAX5 is involved in gene regulation in many tissues. Overall, we found that differences in TFBS numbers could explain variation in gene expression to a modest extent overall, while there are also many TFs explaining a large extent of interspecies variation for specific genes. Next, we sought to identify stromal genes whose expression correlates with the invasive phenotype.

To explore the evolutionary history of genes that have changed expression along with the advent of noninvasive epitheliochorial placentation, we expanded the number of eutherian species in our transcriptomic analysis (2). We isolated and cultured the stromal fibroblasts from endometrial tissues of rabbit, guinea pig, rat, cat, horse, humans, sheep, dog, and cow. This selection of species comes from two eutherian clades, the Euarchontoglires (primates, rodents, and others) and Laurasiatheria (ungulates, carnivores, and others), together forming the clade of Boroeutheria (9). This balanced taxon sample contained species with invasive placentation (human, rabbit, rat, and guinea pig) and less-invasive endotheliochorial or epitheliochorial placentation (cat, dog, horse, sheep, and cow) (Fig. 2A). Our aim was to uncover genes that render hemochorial species to be more invasible (ELI^up genes) or endow the stroma of epitheliochorial species with the resistance to trophoblast invasion (ELI^dn genes) (Fig. 2B). We created a rough invasibility index based on the anatomical presentation of the fetal–maternal interface in these species (9, 29). Specifically, species with epitheliochorial placentation were given a rank of 0, i.e., no invasion (cow, sheep, and horse); carnivores with endotheliochorial placentation were ranked 1 (cat and dog); and rodents with hemochorial placentation ranked as 2 (rabbit, guinea pig and rats), while humans, which along with primates have a particularly invasive hemochorial placentation, were ranked 3 (Fig. 2A). In order to identify genes that correlate with the placental invasibility phenotype, we calculated the correlation between gene expression and the placental invasion phenotype with both a naïve linear model and a phylogenetically informed linear model (Methods) (Fig. 2C and Dataset S2). The phylogenetic linear model corrects for correlations that arise due to shared descent (30). As expected, the phylogenetic linear model showed less significant P values than the naïve model since it allows for fewer degrees of freedom. Nevertheless, the same genes show high correlations in either model showing that confounding effect of phylogenetic relatedness is minimal (Fig. 2C). Therefore, either model can be used to prioritize candidate genes correlated with invasibility.

Fig. 2D shows an example of a gene with higher expression correlated with increased placental invasiveness (defined as ELI^up gene), while Fig. 2E shows an example of a gene whose expression is anticorrelated to invasibility (defined as ELI^dn gene). The ELI^up genes are purportedly proinvasible as their stromal expression increases along with increase in placental invasion, while ELI^dn genes are purportedly antiinvasible as they increase in expression with reduced stromal invasion. Fig. 2 F and G show the top 10 genes ranked in order of their positive and negative influence on invasibility, respectively. We selected a balanced set of 200 statistically most significant genes (all had FDR < 0.05), with 100 of them from the ELI^up set and 100 from ELI^dn. This gene set, termed “ELI gene set,” is a rank ordered set of genes which show consistent differences across species from the least invasive to the most invasive placentation. This gene set was used to identify candidate TFs that regulate stromal invasibility.

Modeling specifically the ELI-specific set of 200 genes, we found that TFBS could statistically explain a much higher fraction of gene expression differences than for all 1–1 orthologous genes (Fig. 3 A and B). For many TFs, both the absolute β values and the fraction of variance explained (Figs. 1I and 3C) were substantially higher than for the whole gene set, with most TFs exhibiting a positive β value. Many of these TFs were related to the myosin enhancer family of TFs, with all major MEF TFs, e.g., MEF2A-D, showing significant effects. In addition, we found that GATA TFs were highly enriched among the TFs showing significant effects. These included GATA1 and GATA2, which are classically referred to as hemotopoeitic TFs (31), and GATA5 and GATA6 (32), which are classically referred to as cardiac specific TFs. Crucially, GATA1 and GATA2 are key TFs up-regulated in cancer associated fibroblasts and induce increased cancer dissemination (33). All four GATA factors are also important in trophoblast development (34), and GATA2 acts as a key TF regulating decidualization (16, 35). These data highlight genomically encoded transcriptional regulation of decidual genes in species with increased endometrium–trophoblast interaction. Notably, we found relatively few TFs with significant negative β values, but they included key developmental TFs like PAX6, SOX10, and SIX2.

TFs displaying the highest statistical significance showed R² values 50 times higher than for similar calculations aggregated over all genes (Figs. 1I and 3 B and C). Individual genes likely to be involved in regulating placental invasibility (such as UACA) can be seen to be affected by differences in the number of binding sites of particular TFs (such as GATA2) (Fig. 3D). Further, as for the whole transcriptome, an aggregated R² value for a TF over all ELI genes masked the substantial effect of TFBS in explaining the variance in gene expression for individual target genes (Fig. 3E). As an example, binding sites for GATA2 could explain more than 50% of the interspecies variance for many ELI genes. Interestingly, GATA2 is acting as an enhancer for some genes and a repressor of gene expression for others (Fig. 3E). For many ELI genes, a substantial extent of the interspecies transcriptional difference is explained by changes in the number of cis-regulatory binding sites (Fig. 3 F and G and Dataset S3). We then sought to experimentally test whether these TFs could affect the downstream gene expression. As predicted, qPCR confirmed that gene KO of TF GATA2 by CRISPR/Cas9 gRNAs significantly reduced the expression of Sash1 and Setdb1, both implicated in preeclampsia (Fig. 3 H and I).

Previously, we experimentally demonstrated that human stromal fibroblasts could be rendered more resistant to trophoblast invasion by modifying gene expression to match that of bovine stromal fibroblasts (2). Here we tested whether interspecies variation in these experimentally validated ELI genes could be explained by cis-regulatory evolutionary changes (Fig. 3J). Indeed we found that among 10 experimentally validated ELI genes (ACVR1, BMP4, CD44, DDR2, LPIN1, MGAT5, NCOA3, STARD7, TGFB1, and WNT5B), many TFBSs correlated with their interspecies expression differences. These included homeobox TFs like PITX1 and GSC; sterol binding factors like SREBF1, SREBF2, and NR4A1; and nuclear signal transducers like NFKB2 and STAT6. Notably, the experimentally validated genes were positively correlated with the invasibility phenotype (i.e., the validated genes are ELI^up genes), and the regulating TFs also had positive β values.

Overall, these data strongly suggest that cis-regulatory genomic sequences explain variation in gene expression associated with increased stromal invasibility or stromal resistance to placental invasion.

Phylogenetic analysis has settled the debate over the evolution of noninvasive epitheliochorial placentation in eutherians, which is now considered to have evolved after the ancestral invasive hemochorial placentation (10, 14, 15). However, anatomical reports suggest that among primates, placentation has continued to evolve into a very aggressively invading form (36, 37). We asked whether cis-regulatory genomic changes in the stroma could specifically explain gain or loss of the phenotype of placental invasibility. Toward this objective, we tested if the ELI genes that correlate with increased (or decreased) placental invasion are differently influenced by TFBS.

We used the top 100 genes which were the most positively correlated (ELI^up genes) and the top 100 genes most anticorrelated (ELI^dn genes) with the placental phenotype (FDR < 0.05 for all selected genes). Separate treatment of either gene set in our model could be influenced by systematic interspecies differences in overall TFBS frequencies. To avoid this artifact, we normalized the frequencies of each TFBS across the species by subtracting the TFBS abundance in the genome from the TFBS number at the promoter of each ELI gene. Overall, we found that a higher number of TFs showed enhancing effects on the expression of ELI^dn genes, while relatively more TFs showed repressive effects on expression of ELI^up genes (Fig. 4A). According to our model, this suggests that the majority of TFs positively regulate stromal resistance to invasion.

For each TF, we analyzed the effect size of TFBS number changes for either gene set (ELI^up and ELI^dn) (Fig. 4B). We found that while few TFs act as either enhancers or repressors for both ELI^up and ELI^dn genes (quadrants I and III), there were also many TFs which acted as antagonistic regulators for ELI^up and ELI^dn genes (quadrants II and IV). The tables in Fig. 4B show examples of TFs which regulate a given ELI^up and ELI^dn gene either as an activator or a repressor. Indeed, for many genes, a large fraction of their interspecies variation is explained by a single TF. The statistical significance of the regulation by each TF for individual ELI genes and for the two sets of 100 ELI^up and ELI^dn genes are listed in Dataset S3.

Quadrant I lists TFs activating both ELI^up and ELI^dn genes and includes many GATA factors and MEF2 factors. All GATA factors and MEF2 factors are key TFs regulating human and mouse ESF decidualization. We speculate that these TFs may regulate resistive genes in ESFs of epitheliochorial placentation and also genes which are present in the more resistive decidualized fibroblasts in hemochorial species. Other TFs also suggest a link to decidualization and preparation for pregnancy including estrogen responsive factor, HOXC10 (38), and early growth response gene 1 (EGR1) (39). EGR1 is induced by placenta growth factor and is an inducer of hypoxia inducible factor 1 (HIF1α) necessary for early implantation (40). Quadrant III shows TFs which acted as repressors for both ELI^up and ELI^dn genes.

In contrast, quadrant II shows TFs which act as repressors for ELI^up genes and activators for ELI^dn genes and thus are positive regulators of resistance to invasion. Since ELI^dn genes are more highly expressed in species with epitheliochorial placentation, these TFs may act as important regulators promoting stromal resistance, by both activating resistive genes and downregulating genes which promote invasion. These TFs include POU2F2/3, SRY, Nkx3.2, and TBX19. Quadrant IV lists TFs, including SIX3 and E2F2, which act as antagonistic repressors of proinvasibility ELI^up and activators of antiinvasibility ELI^dn genes and thus could be positive regulators of invasibility. The presence of antagonistic transcriptional regulators suggests that the different phenotypic outcomes are achieved partially by shared mechanisms. These proinvasibility TFs present attractive targets to inhibit the invasive processes by gene knockdowns. This would, for instance, be desirable in tumor-associated fibroblasts as the KD of such a gene is expected to increase the resistance of these cells to invasive cancer cells. Furthermore, the KD of any one of these TF genes will affect the expression a number of effector genes that influence stromal invasibility.

Our model has identified candidate stromally expressed TFs regulating genes correlated with placental invasion among mammals. We therefore asked if modulating the expression of these TFs could influence stromal invasibility, the key phenotype posited to have evolved according to the ELI framework. Unlike collective migration of cells, stromal invasion is composed of complex interactions between the stromal cells and the invading cells, as well as with the extracellular matrix (1). Invasion into a stromal monolayer is also a kinetically slow process requiring many days of observation to quantitatively resolve small phenotypic differences. We have therefore previously described a quantitative stromal invasion platform composed of juxtaposed populations of invasive cells and stromal fibroblasts on a nanofabricated substrate mimicking the topography of collagen fibers in connective tissues (2, 41). The anisotropic arrangement of the submicron features not only provides a physiological context to observe stromal invasion but also aligns the actomyosin machinery within the cells in one direction, thereby increasing the speed of invasion and sensitivity of observation (Fig. 5A).

We choose two TFs to investigate: GATA2 and TFDP1, a dimerization partner of E2F, which were predicted to significantly enhance expression of ELI^up genes in ESFs (P values for ELI^up genes were 7.987 × 10⁻¹¹ and 1.086 × 10⁻³). We hypothesized that since human ESFs (hESFs) have an overall higher constitutive expression of ELI^up genes compared to their expression in other species, knockdown of these TFs should increase stromal resistance to invasion. Using CRISPR/Cas9-based gene KO, we targeted these TFs in hESFs, as well as in BJ skin fibroblasts. Unperturbed trophoblasts (HTR8), breast epithelial cells (MCF10A1), or malignant melanoma cells (A375), all of human origin, were patterned along with either the endometrial fibroblasts or skin fibroblasts (BJ) (Fig. 5B). Morphological changes in the interface were observed using live cell fluorescence and phase-contrast microscopy for 24 h to identify the extent of aerial invasion into the stroma. It is possible for the invasive fronts to deeply invade the stromal monolayer without a homogeneous aerial movement. Invasive fronts involve specification of a leader cell, which forges its way through the stromal monolayer, followed by other cells maintaining cell–cell connection through adherens junctions, and are commonly described in cancer dissemination and other invasive processes (42, 43) (Fig. 5 C and D).

Using the stromal invasion assay, described earlier (2), we performed experiments to measure the effect of GATA2 and TFDP1 gene silencing in hESFs and human skin fibroblasts BJ cells and its effect on the invasion of trophoblasts, breast epithelial cells, and melanoma cells for 24 h (Fig. 5E and Dataset S3).

We found that both GATA2 and TFDP1 gene KO in the hESFs resulted in significant reduction in the invasion by HTR8 trophoblasts (Fig. 5F). Furthermore, the morphological nature of invasion by HTR8 cells was also influenced by the gene KO. To focus on the depth of stromal invasion by the invasive forks (Fig. 5 C and E), we calculated the percentage of the vertical segment occupied by the invading monolayer at every location in the direction of invasion, showing stark differences in the profile along the depth of invasion by TFDP1 silencing (Fig. 5G). We also identified the tip for each invading fork and calculated the depth forged by each fork from its initial location (Fig. 5H), finding that TFDP1 silencing in hESFs significantly reduced the depth of invasion by HTR8 cells (Fig. 5H).

We also tested the effect of modulating stromal gene expression on invading breast epithelial MCF10A1 cells. Again, silencing of both TFs in BJ fibroblasts showed a significant effect in enhancing stromal resistance to MCF10A1 invasion (Fig. 5I). Notably, TFDP1 silencing in BJ also significantly reduced not only the overall aerial invasion of MCF10A1 but also the depth forged by the invading forks (Fig. 5J). Finally, we also showed that silencing both TFs in BJ cells significantly enhanced their resistance to invasion by malignant melanoma cells, A375 (Fig. 5K).

We previously have proposed that stromal fibroblasts in other tissues may have coevolved with endometrial fibroblasts, resulting in a more resistive stroma in species with epitheliochorial placentation (2). These experiments demonstrate that the TFs predicted by our comparative placentation analysis to regulate endometrial vulnerability to placental invasion also regulate the invasibility of skin fibroblasts by cancer cells.

hESFs were immortalized by the Charles Lockwood laboratory and obtained from the Gil Mor group. Human SFs (BJ) were purchased from American Type Culture Collection (ATCC) (CRL-2522). ESFs from cow (Bos taurus), dog (Canis lupis), guinea pig (Cavia porcellus), horse (Equus caballus), cat (Felis catus), rabbit (Oryctolagus cuniculus), sheep (Ovis aries), and rat (Rattus norvegicus) were obtained by enzymatic digestion of uterine tissues and attachment-based positive selection.

RNA was isolated using RNeasy micro kit (QIAGEN) and resuspended in 15 μL of water. RNA integrity number of over 8 was observed in the samples measured by Agilent Bioanalyzer 2100. mRNA libraries prepared using TruSeq RNA Library (Illumina) were sequenced on Illumina HiSeq to generate 30 to 40 million reads per sample (single-end 75 base pair reads), and a high Q score was observed (Q > 30) for the sequenced data. A transcript-based quantification approach was used to analyze RNA sequencing (RNA-seq) data using the program kallisto (55). The Ensembl release 99 (56) gene annotation model was used to align raw reads. In order to facilitate comparative gene expression across species, a one-to-one ortholog dataset consisting of 8,639 species was formulated (composed of human, rat, rabbit, horse, guinea pig, cat, sheep, cow, horse, dog, and opossum) such that the sum of TPMs across these 8,639 genes for each species totals to 10⁶.

We obtained the gene expression $e_{g,s}$ of gene $g$ in species $s$ in TPM. For each species, we also rated the degree of placental invasibility $r_s$, with a value from 0 to 3. In order to find the stromal genes likely involved in modulating the degree of invasibility, we fit the linear modelwhere ${ϵ}_{g,s}$ is the residual. This is fitted independently for each gene, and positive ${\gamma }_g$ represent genes that are up-regulated in species with more invasive placenta, while negative ${\gamma }_g$ indicates the opposite. The statistical significance is calculated for testing against the null hypothesis of ${\gamma }_g=0$. The P values for all genes are corrected for multiple testing using the false discovery rate method. In addition, this linear model was also fitted using a phylogenetic linear model, with an assumption of Brownian motion processes producing both the traits (gene expression and species invasibility). The phylogenetic linear model was fitted using the phylolm package (30).

For the cell patterning, a PDMS (polydimethylsiloxane) mold, which was fabricated by a stereolithographic plastic mold, was used. To fabricate a PDMS stencil, the monomer and cross-linker (with a ratio of 10:1) of PDMS were mixed and cured at 80 °C for 4 h and then cast in the predesigned mold. The stencil was placed on the substrate after washing with isopropanol and drying with N₂ steam. The device was kept in a vacuum to remove air bubbles under the stencil. The Invasive cells were labeled by cell tracker green and were seeded at a density of 5 × 10⁵ cells and attached to the surface overnight. The stencil was removed carefully using blunt-end tweezers. The unlabeled stromal cells were seeded at a density of 5 × 10⁵ to fill and attach that area covered by the stencil before. The unattached cells were washed off after 5 h of incubation. The plate was mounted on the live microscopy stage.