Research ArticleHEALTH AND MEDICINE

A single-cell survey of the human first-trimester placenta and decidua

See allHide authors and affiliations

Science Advances  31 Oct 2018:
Vol. 4, no. 10, eaau4788
DOI: 10.1126/sciadv.aau4788

Abstract

The placenta and decidua interact dynamically to enable embryonic and fetal development. Here, we report single-cell RNA sequencing of 14,341 and 6754 cells from first-trimester human placental villous and decidual tissues, respectively. Bioinformatic analysis identified major cell types, many known and some subtypes previously unknown in placental villi and decidual context. Further detailed analysis revealed proliferating subpopulations, enrichment of cell type–specific transcription factors, and putative intercellular communication in the fetomaternal microenvironment. This study provides a blueprint to further the understanding of the roles of these cells in the placenta and decidua for maintenance of early gestation as well as pathogenesis in pregnancy-related disorders.

INTRODUCTION

The first-trimester human placenta and maternal decidua interact dynamically in a highly regulated manner to enable establishment of pregnancy; provide physical support and immunologic tolerance; facilitate maternal-fetal transfer of nutrients, waste, and gas exchange; and produce hormones and other physiologically active factors (1). Abnormalities in this nascent period of placental development can manifest in miscarriage and complications later in pregnancy including preeclampsia, fetal growth restriction, and placenta accreta, which are leading causes of maternal and neonatal death (2).

Structurally, the placenta is a complex and heterogeneous organ derived from trophectoderm and extraembryonic mesoderm. It consists of villous projections comprising an outer layer of multinucleated syncytiotrophoblasts (SCTs) overlying villous cytotrophoblasts (VCTs) and an inner core consisting of fibroblasts (FBs), fetal macrophages [Hofbauer cells (HCs)], and capillary networks (3, 4). Placental villi anchor to the maternal decidua or terminate as floating villi within intervillous spaces that are bathed in maternal blood. Extravillous trophoblasts (EVTs) migrate away from the villi to invade the maternal decidua, and a subset of EVTs are involved in the remodeling of maternal spiral arteries and uterine glands. Despite the pivotal roles of the first-trimester human placenta and decidua, we do not currently have high-resolution, single-cell transcriptomic profiles of all subtypes of cells present, providing insights into molecular characteristics, the cellular origin of critical secreted factors, and the cellular origin of placenta-derived extracellular RNAs.

Here, we performed single-cell RNA sequencing (scRNA-seq) analysis of human first-trimester placental and decidual cells using 10x Genomics (referred to as 10x) and Drop-seq platforms. We identified gene signatures, crucial transcription factors (TFs), and abundantly expressed receptors and ligands for each cell type of villi and decidua. We distinguished resting and highly proliferative subtype in villi and decidua and uncovered previously unknown subtypes of the FB-like cells. Collectively, we charted a comprehensive map of cell types in human first-trimester placenta and decidua.

RESULTS

Strategy for scRNA-seq and data integration

We performed scRNA-seq in freshly collected first-trimester (6 to 11 weeks of gestation) villi (n = 8) and decidua (n = 6) samples using a custom-built Drop-seq (5) or commercial 10x platform (6) (fig. S1A and table S1). We obtained 14,341 and 6754 high-quality scRNA-seq profiles from villi and decidua samples, respectively. To evaluate the impact of tissue dissociation on transcriptome profiles, we averaged the gene expression of all villi scRNA-seq datasets, compared it with villi bulk RNA-seq, and found that they correlated well (r = 0.86, Pearson correlation; fig. S1B). Some of the genes that were elevated in scRNA-seq data were JUN, FOS, and HSPA1A, all of which were recently shown to be associated with enzymatic digestion of tissue required for preparation of the single-cell suspension (7). The genes known to be up-regulated in hypoxic environments, such as HIF1A and HIF2A, were expressed at similar levels in bulk and scRNA-seq datasets, indicating minimal or no hypoxia-induced stress–related effects on transcriptome profiles. The high expression level of hemoglobin subunit genes, HBZ and HBE1, in the scRNA-seq dataset likely represents differences in number of erythroblasts (EBs) between the two datasets. Although the 10x method detected more transcripts [unique molecular identifiers (UMIs)] and genes compared to Drop-seq (fig. S1, D and E), we observed good correlation (r = 0.89, Pearson correlation) between the 10x and Drop-seq expression data (fig. S2B). We collectively analyzed datasets from these platforms after cross-platform data integration using recently described Seurat V2.0 method (fig. S1C) (8). Differential gene expression analysis (Wilcoxon rank sum test) revealed nearly identical segregation of gene expression markers between the 10x and Drop-seq datasets (fig. S2, A and C).

Cell types of villi

We applied a graph-based clustering method and identified nine major clusters of cells as visualized by t-distributed stochastic neighbor embedding (t-SNE; Fig. 1A) (8). Differential gene expression analysis guided by established markers revealed these cell clusters as trophoblasts, VCTs, SCTs, EVTs, three FB cell clusters (FB1 to FB3), vascular endothelial cells (VECs), EBs, and HCs. VCTs, SCTs, and EVTs constituted 41% of cells dissociated from villi tissue, with VCT (90%) being the most abundant trophoblast, followed by SCT (7.7%) and EVT (2.3%; Fig. 1D and fig. S3A). FBs constituted the second most abundant cell type, representing 35% of cells in villi, followed by VECs (9.2%), EBs (8.4%), and HCs (6.8%).

Fig. 1 Single-cell expression atlas of first-trimester villi.

(A) Cell type assignment following t-SNE–based visualization of expression differences for 14,341 single cells using established lineage markers. (B) Dot plots showing expression of known lineage markers and coexpressed lineage-specific genes. (C) Feature plots showing proliferating subpopulations of VCT, EB, FB1, and FB2 cells based on coexpression of proliferation marker MKI67, TOP2A, TK1, and PCNA. (D) Box plots highlighting the contribution of individual samples (n = 8) to each cell cluster. (E) TF enrichment analysis showing the most abundant (maximum of 10) and specific of TFs of major cell groups and individual cell types. (F) Immunofluorescence staining for FB2-specific REN (green) and pan-FB marker VIM (red). Scale bars, 25 μm.

Trophoblasts. Trophoblasts share expression of KRT7 and PERP across all subtypes and can be further subclassified into VCT, SCT, and EVT by sublineage markers such as PARP1, ERVFRD-1, and HLA-G, respectively (Fig. 1B) (9). SCTs express retroelement-like genes including ERVW-1, which is important for syncytium formation. CGA, which codes for the common α subunit of the hormones hCG (human chorionic gonadotropin), LH (luteinizing hormone), FSH (follicle-stimulating hormone), and TSH (thyroid-stimulating hormone), is primarily expressed by SCTs and EVTs. The aromatase CYP19A, which converts androgens to estrogens, is specifically expressed by SCTs, whereas HSD3B1, which is essential for the production of progesterone for pregnancy maintenance, is mainly expressed by EVTs and SCTs. SCTs also express the transforming growth factor–β superfamily member GDF15, which is elevated in maternal plasma during pregnancy and dysregulated in preeclampsia and diabetic pregnancies (10). We observed a series of nutritional, steroid, and hormone solute carrier (SLC) family proteins in all placental cell types. SLC27A2 was specifically expressed by VCTs, SLC22A11 was expressed by both SCTs and VCTs, SLC13A4 and SLC52A1 were highly expressed by SCTs but showed negligible expression in other placental cell types, whereas SLCO4A1 was predominantly expressed by EVTs. We also observed that LCMT1-AS2, an antisense noncoding RNA to LCMT1, and INSL4, an insulin family hormone, were abundantly and specifically expressed by SCTs.

We characterized EVTs by the expression of HLA-G, DIO2, and LAIR2 (Fig. 1B) (11), but we also found novel EVT markers, including the serine peptidase HTRA4, previously reported to be up-regulated in the early onset but not late onset of preeclampsia along with increased levels in maternal serum (12). The extracellular matrix (ECM) glycoprotein MFAP5, a component of microfibrils shown to induce cell motility and invasion of cancer cells (13), was also specifically expressed by EVTs, consistent with a role in EVT invasion and mobility. Last, we identified EVTs and VECs as the primary sources of prostaglandin catabolizing enzyme HPGD. Prostaglandins play crucial roles in the initiation and progression of labor in most mammalian species (14). We performed differential gene expression analysis of EVT and VCT cell types (Wilcoxon rank sum test) and compared the top differentially expressed genes (>1.5 log2 fold change) with the list provided by Apps et al. (15). We observed that established EVT markers such as HLA-G, LAIR2, and DIO and VCT markers such as PAGE4, PEG10, and PARP1 were identified in both the studies (fig. S9). In addition, there were several other genes that were exclusively identified by both our and Apps et al.’s study.

FB-like cells. Three FB subtypes, referred to as FB1, FB2, and FB3, were identified by gene expression. FB1 formed the most abundant FB cluster (58.2%), followed by FB2 (27.4%) and FB3 (14.4%; Fig. 1D and fig. S3A); FB1 and FB2 also comprised proliferating cells (Fig. 1C). FBs share expression of ECM genes such as COL1A1, COL1A2, and COL3A1 (Fig. 1B) and smooth-muscle actin ACTA2. Furthermore, first-trimester FBs are the primary source of DLK1 (Fig. 1B), an imprinted gene encoding an endocrine signaling molecule present at high concentration in maternal circulation during late pregnancy, and its level is strongly associated with fetal growth in mouse and humans (16). Last, the placenta-specific EGFL6 gene, known to promote endothelial cell migration and angiogenesis, was also FB-specifically expressed. FB1 and FB3 showed a characteristic resemblance to myofibroblasts by expressing DES, MFAP4, OGN, and S100A4 genes. FB3 additionally expressed proinflammatory genes such as IL6, PTGDS, CFD, CXCL2, and BDKRB1. FB2s are characterized by the unique expression of REN and AGTR1, genes implicated in the regulation of blood pressure, sodium, and fluid homeostasis, as well as IGFBP7 and AREG. Immunostaining identified the REN-expressing FB2s in the stroma of villi (Fig. 1F and fig. S4).

Other cell types in villi. VECs form the vasculature of the placenta and were identified by the expression of canonical markers including CDH5 and PECAM1 (Fig. 1B). EBs specifically expressed hemoglobin subunit genes such as HBE1, HBG1, HBA1, HBM, HBG2, HBA2, and HBZ. EBs also expressed sialic acid–rich erythrocyte membrane proteins glycophorins (GYPA, GYPB, and GYPC). HCs abundantly and specifically expressed established macrophage markers (CD163, CD14, CSF1R, and CD68), the inflammatory cytokines CCL4 and CXCL8, and allograft inflammatory factor 1 (AIF1). The maternal immunoglobulin G transporter FCGRT was primarily expressed in HCs, and to some extent by FBs and VECs, contrary to an earlier report indicating its expression in SCTs (17). More detailed analysis of each cluster of cells furthermore identified proliferating and nonproliferating subpopulations within VCT, EB, FB1, and FB2 cells based on coexpression of the proliferation markers MKI67, TOP2A, TK1, and PCNA (Fig. 1C).

Cell types of the decidua

Graph-based clustering analysis identified 11 distinct cell clusters characterized by the expression of lineage markers specific for decidualized stromal cells (DSCs), two distinct decidual FB populations (FB1 and FB2), smooth muscle cells (SMCs), endometrial epithelial cells (EECs), two populations of natural killer cells (NK1 and NK2), antigen-presenting cells (APCs), T cells (TCs), lymphatic endothelial cells (LECs), and VECs (Fig. 2A). Decidual tissue was composed of about 48.7% of cells expressing high levels of ECM genes such as COL1A1, COL1A2, COL3A1, and DCN including DSC, FB1, FB2, and SMCs (Fig. 2, B and D, and fig. S3B). Leukocytes (NK1, NK2, APCs, and TCs) were accounted for 40.4% of cells in the decidua, and the remaining cells were 4.9% EEC, 2.1% LEC, and 3.9% VEC.

Fig. 2 Single-cell expression atlas of first-trimester decidua samples.

(A) Cell type assignment using established lineage markers following t-SNE–based visualization of 6754 single cells. (B) Expression data dot plots of known lineage markers and coexpressed lineage-specific genes. (C) Left: t-SNE visualization of APCs. Right: Differentially expressed genes in APC subpopulations. Wilcoxon rank sum test. (D) Box plots represent the contribution of samples (n = 6) to each cell cluster. (E) TF enrichment analysis showing the most abundant (maximum of 10) and specific expression of TFs by major cell groups and individual cell types.

ECM-expressing cells. DSCs represented 27% of ECM-expressing cells and are defined by the expression of not only PRL and IGFBP1, established markers of decidualization (Fig. 2B and fig. S3B) (18), but also APOA1, CHI3L2, SERPINA3, IL1B, and PROK1. Decidual FB1 (44.5%) and FB2 (21.5%) cells are differentiated by the unique expression of ARC, GEM, BDKRB1, and PLIN2 in FB2. Vascular SMCs (7%) are differentiated on the basis of their expression of MYH11, ACTA2, and RGS5 in addition to PI15 and NDUFA4L2.

To address stromal cell differentiation, we subjected 1524 stromal cells of decidua (DSC, FB1, and FB2) profiled by 10x to diffusion map analysis (destiny package) and placed these cells in pseudotemporal order (19). On the basis of the expression of 1751 variable genes according to graph-based clustering (Seurat), we observed two trajectories originating from the FB1 population toward DSCs and FB2s with a gradual increase in expression of PRL and GEM, respectively (Fig. 3A). Further analysis of gene expression on each branch revealed graded expression for several genes. Most notably, IGF1 was specific to the FB2 branch, whereas insulin-like growth factor–binding proteins (IGFBPs) were expressed in a cell lineage–dependent manner along pseudotemporal trajectories (IGFBP5 on FB2, IGFBP6 on FB1, and IGFBP1 on DSC), indicating distinct regulation of IGF1 function in the decidual microenvironment (Fig. 3, A and C). Decidualization is characterized by extensive tissue remodeling to facilitate implantation of the blastocyst and to regulate EVT invasion through regulation of ECM genes (20). We found that although most collagen genes showed abundant expression on FB1 cells, COL4A5, COL14A1, and COL4A6 were uniquely expressed on the DSC trajectory (Fig. 3C and fig. S8). Polyploidization is the hallmark of mature DSCs, which is also characterized by an increased number of genes expressed (21). In concordance, we observed that DSCs expressed a larger number of genes and transcripts (P < 0.05, unpaired t test) compared to the FB1 and FB2 cells (Fig. 3B).

Fig. 3 Pseudotemporal ordering of stromal cells of the decidua.

(A) Diffusion maps of 1524 stromal cells from decidua showing expression pattern of genes over pseudotime trajectory. (B) Bar graphs representing number of UMIs and genes for FB1, FB2, and DSC cell types (n = 5, means ± SEM; *P < 0.05, unpaired t test). (C) Expression pattern over pseudotime for the top 50 most variable genes sorted along the DC1 coordinate. Box plots showing average expression level of genes.

Endometrial epithelial cells. EECs are characterized by the abundant expression of PAEP and epithelial markers such as KRT8 and CDH1 (Fig. 2B) (22), which play a crucial role during embryonic implantation (20). We observed that IGFBP1, which is expressed by DSCs, was also expressed by EECs and that EECs specifically express the monoamine transporter SLC18A2, along with GPX3, PCCA, SLPI, DEFB1, CP, RIMKLB, and LCN2. We also identified three long noncoding RNAs (LINC01320, LINC01502, and LINC01541) exclusively expressed by EECs.

Leukocyte subtypes of the decidua. NK cells are prominent in the early decidua but vanish from the endometrium before menses in the absence of pregnancy (23). NK cells expressed the canonical markers NKG7 and NCAM1 (Fig. 2B), the cytotoxic granzyme GZMA and GZMB, perforin PRF1, and the inflammatory chemokines XCL1 and XCL2. We differentiated a resting (NK1, 58% of leukocytes) and a proliferating (NK2, 8.9% of leukocytes) subpopulation based on their expression of MKI67 and TOP2A (Fig. 2D and fig. S3B). Gene ontology enrichment analysis comparing NK2 to NK1 gene expression revealed that the top enrichment processes were cell cycle process, cellular component organization, and cell proliferation, further confirming the activated and proliferative status of NK2 cells (Fig. 4C).

Fig. 4 Differential gene expression of villi and decidual cell type averaged profiles.

(A) Annotation of fetal and maternal origin for each cell type based on the expression of DDX3Y, EIF1AY, and XIST based on male fetus data (villi, n = 7; decidua, n = 4). (B) Unsupervised clustering analysis of the top 60 genes expressed in each ECM-expressing cell types of villi and decidua (except SMCs). (C) Left: Heatmap reporting average expression level [log2(TPM + 1)] of the top 10 differentially expressed genes in decidual NK1 and NK2 (top) and NK lineage marker genes (bottom). Right: Gene ontology (GO) terms enriched for genes in each of these panels. (D) Scatter plots depicting villi- versus decidual-specific gene expression [log2(TPM + 1)] based on log2 fold change (FC). Left, VECs; right, MACs/HCs.

TCs occurred at 14.2% frequency among all leukocytes and were annotated on the basis of expression of CD3D and CD69 (Fig. 2B). TCs also showed elevated levels of expression of genes involved in cell-cell signaling such as LTB, CXCR4, CCL5, and IL7R. Expression of IFNG, GZMA, and GZMK established these cells as cytotoxic TCs.

APCs constituted 19.4% of decidua leukocytes according to expression of CD74, human leukocyte antigen (HLA) class II α chain (HLA-DRA and HLA-DPA1) and β chain (HLA-DRB1 and HLA-DPB1) genes (Fig. 2B). Subclustering analysis of APCs identified three subtypes (Fig. 2C), further classified as maternal macrophages (MACs) and dendritic cells (DC1 and DC2). MACs specifically expressed MS4A4A, STAB1, SEPP1, and MS4A7, characteristic for “alternatively activated” M2 macrophages. Dendritic cell types DC1 and DC2 were of the CLEC9A+ and CD1C+ subtypes, respectively, as described in a recent scRNA-seq study of human blood (24).

Lymphatic and vascular endothelial cells. LECs, identified by the expression of LYVE1, and VECs, identified by the expression of CDH5 (Fig. 2B), represented 2.1 and 3.9% of all decidual cells, respectively (Fig. 2D and fig. S3B). Expression of STMN1, CD9, FABP5, ADM, and C10orf10 was characteristic for LECs, while VECs expressed ACKR1, PCDH17, FAM167B, AQP1, HP, and IFI27.

TFs of cell types in villi and decidua

We used TFs with expression greater than 0.75 TPM (transcripts per million) and expressed by more than 10% cells for enrichment analysis. Within villi, EBs, FB3, and VCTs were distinguished by the expression of the Krüppel-like family (KLF) of TFs KLF1, a known regulator of erythropoiesis (25), KLF4, and KLF5, respectively (Fig. 1E). NR4A2, involved in the macrophage-mediated inflammation process, was expressed by HCs. In addition, several other TFs were enriched in a cell type–specific manner.

Within the decidua, ELF3, PAX8, and EHF, implicated in epithelial cell survival and function, were abundant and specific in EECs (Fig. 2E). Immune cells, such as NK cells, TCs, MACs, and DC1 cells, preferentially expressed known immune cell–specific TFs EOMES, BAFT, KLF2, and ZNF366, respectively.

Differential gene expression analysis of villi and decidua cell types

To differentiate embryonic/fetal and maternal cellular origin, we studied pregnancies with male fetuses and took advantage of sex-specific expressed Y chromosome–encoded genes (DDX3Y and EIF1AY) and female-specific XIST, which is involved in X chromosome inactivation. Using male embryo datasets, we identified that all cell types identified in the villi sample expressed male-specific DDX3Y and EIF1AY and virtually no XIST, all indicative of their embryonic/fetal origin, whereas all cell types identified in the decidua samples showed an opposite expression pattern corresponding to a maternal origin (Fig. 4A).

On the average expression profiles of all FB and DSC subtypes, we performed unsupervised clustering analysis of the union of the top 60 expressed (Fig. 4B). Villi FB1, FB2, and FB3 expressed DLK1, HAPLN1, GPC3, EGFL6, and S100A10 among other genes. Decidua FB1, FB2, and DSC expressed APOD and SPARCL1, and DSC showed very high expression of PRL and IGFBP1. On the basis of these observations, we indicate that FB-like cells adopt unique fates based on their tissue environment.

Similarly, we compared the VECs of decidua and villi (Fig. 4D, left). Most notably, EGFL7, known to be involved in vascular tube formation (26), was abundantly expressed by villi VECs, along with focal adhesion structure and signaling genes CAV1, MARCKS, and CTNNB1, and tissue factor pathway inhibitors TFPI and TFPI2. In contrast, IGFBP7, a known inhibitor of angiogenesis of VECs (27), was expressed by decidual VECs. Furthermore, interferon-stimulated genes such as IFI27 and IFITM1 were abundantly expressed in decidual VECs, as was ACKR1, a nonspecific receptor for several chemokines, and A2M, an inhibitor of inflammatory cytokines.

The comparison of expression patterns of MACs with their embryonic/fetal counterpart HCs as a subpopulation also revealed notable differences. Embryonic/fetal HCs expressed higher levels of SPP1, PLIN2, HMOX1, CD36, and LYVE1 (Fig. 4D, right), whereas MACs showed higher expression of HLA genes (HLA-DRA, HLA-DPA1, HLA-DRB1, HLA-DPB1, HLA-DQA1, and HLA-DMA) and invariant chain CD74, indicating strong antigen-presenting capacity of maternal MACs.

Potential ligand-receptor interactions

To investigate intercellular communications in tissues, we visualized average expression levels of most abundant ligands and their cognate receptors for each cell type’s ligand-receptor pair (Fig. 5). In villi, membrane-anchored ligand DLK1 showed the highest expression (>75 TPM), restricted mostly to stromal cells, whereas the receptor, NOTCH3, was expressed by FB2s, consequently implicating this interaction in the maintenance of the stroma (Fig. 5A) via paracrine signaling. High levels of plancental growth factor (PGF) ligand and soluble fms-like tyrosine kinase 1 (FLT1) receptor in maternal serum are predictive of preeclampsia (28). PGF-FLT1 interactions are abundant in villi initiated by PGF expression from SCTs and EVTs, and FLT1 expression in EVTs, VECs, FB2s, and FB3s. IL6, implicated in the regulation of acute inflammation, was expressed in FB3s, suggestive of interaction with its receptor, IL6ST, expressed in EVTs, and stromal cells (FB1, FB2, and FB3). CCR1, a regulator of trafficking of macrophages (29), was expressed in HCs, and the corresponding ligand, CCL2, was abundantly expressed in VECs and FB2s, suggesting its involvement in chemoattracting HCs in villi. In addition, we observed high expression of CCL2, CCL4, and CCL5 by HCs, suggesting a CCR1-dependent positive feedback loop for HCs infiltration.

Fig. 5 Ligand-receptor interaction analysis.

Putative signaling within villi (A), decidua (B), and between villi and decidua cell types (C) with the size of the arrow stem proportionate to expression levels of the ligands. All the arrows are pointing to the receptors. Only ligands with >5 TPM and receptors with >1.5 TPM average expression were used for interaction analysis.

In decidua, CCL21 was the most abundant cytokine expressed by LECs and is likely involved in chemotaxis of NK1 and NK2 cells into decidua binding to their CXCR3 receptor (Fig. 5B). In decidua, the chemokines that potentially undergo transcytosis via ACKR1 (30) were CCL2 in MACs and DC2, CCL5 in TCs and NK2 cells, and CXCL8 in MACs, DC2, and FB2 cells. PRL-PRLR communication was restricted to stromal cells. IFNG binds to heterodimeric IFNGR1 and IFNGR2 receptor to regulate several biological processes including inflammation. IFNG expression was restricted to TCs, whereas IFNGR1 and IFNGR2 were expressed on MACs and DC2, indicating immune cell cross-talk in the decidual microenvironment.

In vitro, trophoblasts show concentration-dependent migration in response to CCL21 (31), which we identified to be uniquely expressed by decidual LECs (Fig. 5C). Its receptor CCR7 was expressed by SCTs and VCTs, indicating a likely role of CCL21-CCR7 interactions in trophoblast migration in decidua.

DISCUSSION

The placenta and maternal decidua undergo remarkable transformations to enable establishment and support of pregnancy. Our high-resolution single-cell survey of the first-trimester human placenta and maternal decidua established a cellular transcriptomic atlas from thousands of cells of each tissue type and integrated the data across two recent scRNA-seq platforms. Our study complements previous attempts of scRNA profiling that did not evaluate the decidua and the focus was limited to second- and third-trimester placenta when most dynamic changes have already been completed (32, 33). The overall proportion of trophoblasts in the first (41%) and third trimesters (36%), however, is comparable. Across trimesters, VCTs remained the predominant trophoblast cell type, followed by SCTs and EVTs, except when the tissue sample was retrieved from the deeper fetomaternal interface, resulting in an increase in the EVT proportion in parallel to an increase in decidual cell types. Both first- and third-trimester VCTs showed a subpopulation of proliferating cells expressing MKI67 and TOP2A, likely reflecting the continuous process of VCT to EVT differentiation throughout the pregnancy. The third-trimester EVTs highly expressed the protease MMP11 (28), specific for substrates such as IGFBP1 and COL6A3, both of which are highly expressed by decidual stromal cells (Fig. 3C and fig. S8). Our results showed that first-trimester EVTs expressed MMP12, known for its role in dampening the inflammation process, rather than MMP11, which degrades collagen (data file S1) (34). These results indicate a possible role shift of EVTs from being anti-inflammatory in the early stage to highly invasive in the later stage of pregnancy. The addition of the decidua allowed us also to examine the cellular composition and dynamic interaction with the placenta. Although exclusively in silico, by restricting analysis to the only highly expressed receptor-ligand interactions, these predictions remain biologically important. We also characterized the relative proportion of each cell type in villi and decidual tissue. However, the estimate is likely to be influenced by factors including sensitivity to dissociation methods and cell size. The proportion of large, multinucleated SCTs is likely underestimated because of their inefficient encapsulation during the droplet generation process. Despite these challenges, we captured 449 high-quality SCTs identified by established and novel marker gene expression.

In summary, our analysis yielded transcriptome definitions of 20 major cell populations (9 from villi and 11 from decidua), insights into the complex and layered regulatory code for maintaining differentiated cells, and a map of the interactome between the most abundantly expressed ligands and receptors within and between villi and decidua cell types. Our analysis also revealed an unexpected cellular subtype complexity within each tissue type, including the presence of the distinct FB populations within the placenta and the decidua. This survey will guide future molecular elucidation of the dynamic changes within and between the placenta and decidua and enhances our understanding of how dysregulation of these processes including genetic changes may lead to miscarriage and pregnancy complications.

MATERIALS AND METHODS

Tissue collection and dissociation

First-trimester elective termination tissue samples were obtained either via manual vacuum aspiration or dilation and curettage and temporarily floated in cold water. Gestational sac, chorionic villi, decidua, and blood clots were identified microscopically. The chorionic villi and decidua were separated using forceps and scissors and then transported to the laboratory in phosphate-buffered saline (PBS) maintained at 37°C. For cell dissociation, 200 μl of freshly prepared dissociation solution composed of Liberase TL (2 mg/ml; Sigma-Aldrich) and 1800 μl of PBS were added to 0.1 to 0.3 g of tissue samples and incubated at 37°C for 15 min. The tissue was then manually disaggregated using a 2-ml syringe carrying a 16-gauge needle and pulling up and down the solution gently for 10 times, followed by two further rounds of needle aspiration in 15-min intervals now using an 18-gauge needle. Thereafter, 2 ml of fetal bovine serum (FBS) was added, and cells were passed through a Falcon 40-μm cell strainer (#352340, Corning) and collected by centrifugation at 300g for 5 min. Cells were washed once and resuspended in 100 μl of freshly prepared PBS–bovine serum albumin (BSA) (1× PBS and 0.04% BSA). Cell viability was assessed using Trypan blue (#1450013, Bio-Rad) exclusion method on a TC20 automated cell counter (Bio-Rad). Placenta_23 villi were processed on both Drop-seq and 10x platforms as the sample yielded enough cells. For Placenta_22 decidua, we processed two different tissue samples (ID: P2D_DS and P5D_DS) from the same individual.

Single-cell capture, library preparation, and sequencing

scRNA-seq libraries were generated using a Drop-seq protocol or the Chromium Single Cell 3’ Reagent Kit (10x Genomics). Drop-seq was performed as described earlier (5). For Drop-seq, single-cell suspensions were prepared in 1× PBS-BSA at a concentration of 100,000 cells/ml. Barcoded beads (ChemGenes) were resuspended in lysis buffer at concentration of 120,000 beads/ml. Syringes loaded with suspension of beads, suspension of cells, and oil were connected to a custom-built microfluidics chip (FlowJEM), and monodispersed droplets were generated. Droplets were lysed, complementary DNA (cDNA) was generated and amplified, and quality was assessed using the TapeStation (Agilent 2200). Thereafter, the cDNA was tagmented and amplified using the Nextera XT DNA Sample Prep Kit (Illumina) using primers described previously (Drop-seq Laboratory Protocol, version 3.1; http://mccarrolllab.com/dropseq/). The library was purified and sequenced using custom Read1 primer on the Illumina NextSeq 500 platform (using High Output v2 kit, Illumina) as follows: 20 base pairs (bp) (Read1) and 60 bp (Read2). For 10x scRNA-seq, the procedure was performed according to the manufacturer’s instruction using Chromium Single Cell 3’ Reagents Kits V2 (10x Genomics). The library was sequenced on Illumina HiSeq 2500 platform as follows: 26 bp (Read1) and 98 bp (Read2).

Bulk RNA-seq and analysis

Frozen tissue (350 mg) of postpartum first-term placenta in 2-ml microcentrifuge tubes was lysed by milling (MM301, RETSCH Mill) for 5 min at 27 cycles/s using 3-mm Tungsten Carbide Beads (#69997, QIAGEN) in 359 μl of ED2S buffer [36 mM citric acid, 44 mM NaOH, 0.4% (w/v) sacosyl, 1.5 M guanidine isothiocyanate, 183 mM NaCl, and 48 mM 2-mercaptoethanol in phenol]. After lysis, the suspension was centrifuged at 13,000g, 1 ml of the supernatant was transferred into a new tube, and 298 μl of PBSP [125 μl of 30% (w/v) sodium dodecylsulfate, 66 mM tris-HCl (pH 7.5), 19.8 mM EDTA, 265 mM 2-mercaptoethanol (pH 8.0), and 175 μl of PBS] was added, followed by mixing, incubation at 60°C for 90 s, and cooling on ice for 90 s. Phase separation was induced by adding 126 μl of chloroform, vortexing, and spinning at 21,000g for 5 min at room temperature. Aqueous supernatant (650 μl) was transferred into a new tube, and 1300 μl of VB2G [98.2% isopropanol, 7.2 mM MgCl2, 2.4 mM CaCl2, 1 M guanidinium thiocyanate (GITC), and 5.0 mM tris(2-carboxyethyl)phosphine (TCEP)] was added, followed by sample loading to QIAGEN RNAEasy MiniElute spin columns using a vacuum manifold. After loading, columns were washed twice with 970 μl of buffer EWL (18 mM NaCl, 2.7 mM MgCl2, 0.9 mM CaCl2, 0.5% Triton X-100, 360 mM GITC, and 5 mM TCEP), once with 970 μl of 100% (v/v) ethanol, and twice with 500 μl of 80% (v/v) ethanol. Membranes were dried by centrifugation at 17,000g for 5 min, and samples were eluted in 20 μl of 10 mM tris-HCl (pH 7.4) by centrifugation under the same conditions. Total RNA (100 ng) was used to generate RNA-seq libraries using Illumina TruSeq stranded mRNA LT kit. Libraries were denatured and sequenced on Illumina HiSeq 2500 sequencer to generate 101-bp paired-end reads.

Bulk RNA-seq analysis was performed using standard RNA-seq pipeline. Briefly, after quality control, reads were aligned to human (hg38) genome reference using STAR aligner (STAR_2.5.1a), allowing no more than three mismatches. Counts for genes were generated using the function featureCounts in the R/Bioconductor package “Rsubread,” followed by TPM calculation.

scRNA-seq data processing

The Drop-seq and 10x sequencing data were processed using the standard pipeline (Drop-seq core computational protocol V1.2, http://mccarrolllab.com/dropseq/) with minor modifications. For Drop-seq Read1, bases 1 to 12 were tagged with cell-barcode “XC,” and bases 13 to 20 were tagged with UMI “XM.” For 10x Read1, bases 1 to 16 were tagged with cell barcode XC, and bases 17 to 26 were tagged with UMI XM. Read2 was trimmed at the 5′ end to remove any adaptor sequence, and the 3′ end was trimmed to remove poly(A) sequences of length six or more and then aligned to human (hg38) genome reference using STAR aligner (STAR_2.5.1a), allowing no more than three mismatches. Gene expression matrix was then generated using the “MIN_BC_READ_THRESHOLD=2” option to retain UMIs with two or more reads.

Canonical correlation analysis, Seurat alignment, and t-SNE visualization

Gene expression analysis and cell type identification were performed independently for villi and decidua samples using Seurat V2.0 (8). The gene expression matrix for each sample was generated, and ubiquitously expressed ribosomal protein–coding (RPS and RPL) and MALAT1 noncoding RNA genes were removed. Seurat objects were created for individual samples. Only those genes that were expressed in more than three cells and cells that expressed more than 100 genes were retained. All 10x Seurat objects for individual samples were merged into one 10x combined object (MergeSeurat), followed by scaling data (ScaleData) and finding variable genes (FindVariableGenes). All Drop-seq Seurat objects for individual samples were processed through similar steps as described above to generate a single Drop-seq combined object. Next, the union of the top 2000 variable genes for each, 10x and Drop-seq, combined objects was used to perform canonical correlation analysis (CCA) between 10x and Drop-seq datasets. Then, CCA subspaces were aligned using 1:16 CCA dimensions, which was followed by integrated t-SNE visualization for all cells. The expression of established lineage marker genes was used to assign cell types. Cells with >25% mitochondrial content (based on UMIs) was considered of poor quality and set aside in the second round of analysis. For villi, an average of 6% of cells in each lineage were poor-quality cells except for one cluster of cells comprising >50% poor-quality cells. This cluster also showed fewer UMIs than other lineages, was devoid of expression of cell type–specific genes, and therefore was removed in the second round of analysis. Next, we analyzed each individual cluster separately following similar strategy as described above. Individual cluster analysis identified ~4.5% cells to be doublets based on expression of multiple lineage markers that were removed in the final round of analysis. Subclustering analysis of only those clusters showing subpopulation was reported. In this process, we also identified proliferating subclusters within VCT, EB, FB1, and FB2 clusters of villi that showed expression of MKI67, TOP2A, TK1, and PCNA genes.

Antibody validation for immunostaining

For cloning of lentiviral overexpression constructs, first, the wild-type cDNA of REN (renin) was cloned into pLenti-C-Myc-DDK-IRES-Puro (OriGene) plasmid by digestion with Asc I and Mlu I. For the lentivirus production, 5 × 105 human embryonic kidney (HEK) 293T/17 cells were cultured in one well of a six-well plate in high-glucose Dulbecco’s modified Eagle’s medium supplemented with 10% (v/v) FBS. Twenty-four hours after seeding, the HEK293T/17 cells were transfected using Lipofectamine 2000 (11668019, Thermo Fisher Scientific) and 1.2 mg of pLenti-C-Myc-DDK-IRES-Puro-REN plasmid, 1.2 mg of psPAX2 (#12260, Addgene), and 0.8 mg of pMD2.G (#12259, Addgene) packaging plasmids. Forty-eight hours after transfection, the viral particle containing supernatant was collected, cleared by centrifugation at 500g for 10 min, and filtrated (45 mm; 09-720-005, Thermo Fisher Scientific). The viral titer was estimated using the Lenti-X GoStix Kit (#631244, Clontech). For viral transduction, Flp-In 293 T-Rex cells were infected with the virus at multiplicity of infection ~0.4 and ~1. Lentivirus-transduced cells were selected with puromycin (2 mg/ml) starting 2 days after viral infection. Expression of the C-terminal Myc-DDK–tagged REN protein was verified by Western blot using anti-MYC antibody.

Immunostaining of cell pellets and placenta tissue

Cell pellets were fixed in 4% paraformaldehyde overnight at room temperature and processed for paraffin embedding using Leica ASP6025 tissue processor (Leica Biosystems). Freshly cut 5-μm paraffin sections were stained for Renin polyclonal (5 μg/ml; #AF4090, R&D Systems) for 1 hour and using 10 min of 1:200 Tyramide Alexa Fluor 488 detection (#T20948, Life Technologies) on Leica Bond Protocol F. The sections were pretreated with Leica Bond ER2 Buffer (Leica Biosystems) for 20 min at 100°C before staining.

First-trimester placental villi tissues were fixed with 10% neutral buffered formalin (Azer Scientific) at room temperature for 24 hours. Samples were transferred to 70% ethanol, followed by paraffin embedding. Freshly cut 5-μm paraffin sections were stained for sequential double immunofluorescence on Leica Bond RX (Leica Biosystems) with Vimentin mouse monoclonal (0.05 μg/ml; #VP-V684, Vector Labs) for 1 hour and using 10 min of 1:200 Tyramide Alexa Fluor 647 detection (#T20951, Life Technologies) on Leica Bond Protocol F, followed by Renin polyclonal (0.2 μg/ml; catalog no. AF4090, R&D Systems) for 1 hour and using 10 min of 1:200 Tyramide Alexa Fluor 488 detection (#T20948, Life Technologies) on Leica Bond Protocol F. The sections were pretreated with Leica Bond ER2 Buffer (Leica Biosystems) for 20 min at 100°C before each staining. Images were recorded on the Olympus VS110 and processed using Visiopharm Integrated Systems software.

Fetal and maternal origin analysis

Only male fetus samples were used for this analysis, and female fetus samples (P3V_DS, P2D_DS, and P5D_DS) were removed to facilitate assignment of embryonic and maternal origin cell types. Gene expression matrices were used to perform CCA, Seurat alignment, and t-SNE clustering as described earlier. Expression of Y chromosome–encoded genes (DDX3Y and EIFAY) was used to assign fetal origin, whereas expression of XIST gene was used to assign maternal origin.

FB gene expression analysis

Average gene expression was calculated for each FB subtype. Next, each subtype expression was normalized to 10,000 to create TPM-like values, followed by transforming to log2(TPM + 1). Log-transformed values for the union of the top 60 genes expressed in each cell cluster were used to perform hierarchical clustering by pheatmap in R using Euclidean distance measures for clustering.

TFs analysis

We created average gene expression profiles (TPM) for each individual cell types. To identify cell type–specific TFs, we used a list of known TFs described earlier (35) and crossed it with TPM-normalized expression by cell types. Only TFs with expression greater than or equal to 0.75 TPM and expressed in more than 10% of cells were considered. For each selected TF, a Z score based on log2(TPM + 1) was calculated, and to find TFs highly specific to an individual cell type, we applied a threshold of 0.6 to the ratio of second maximal expression value to the maximal. For the graphical representation of results, we selected top TFs (maximum of 10) with the highest Z score for each cell type.

Gene ontology analysis

We performed differential gene expression analysis of all decidual cell types using Wilcoxon rank sum test inbuilt in Seurat package. Top 10 differentially expressed genes for NK1 and NK2 clusters were subjected to DAVID (https://david.ncifcrf.gov/) (36, 37) gene ontology functional annotation with default parameters. We ranked terms according to the P value [also known as EASE (Expression Analysis Systematic Explorer) score], which is derived from a modified Fisher’s exact test.

Diffusion map analysis of stromal cells of decidua

We performed Seurat clustering analysis on 1524 stromal cells (DSC, FB1, and FB2) of P6D_10x sample and identified 1751 variable genes (FindVariableGenes). We used log-transformed gene expression matrix of variable genes to generate diffusion map by destiny R package. To identify genes on pseudotime, we used a simplified approach using only the first coordinate (DC1) of diffusion map. We divided the range of DC1 coordinate into 50 bins of equal length and calculated mi,k, mean log2(TPM + 1), for each ith gene and jth cell combination falling into the kth bin. We filtered out genes for which mi,max ≤ 1, mi,max = max(mi,1, mi,2, …, mi,k). For demonstration, we selected the top 50 variable genes by sorting by σi2, variance for mi,k for ith gene among all bins. The order of genes was determined by the index of Medi (median for mi,k), sorted in ascending order. A standard score or Z score, zi = (mi,kMk)/σi, where Mk = mean(mi,1, mi,2, …, mi,k), was calculated for each gene separately.

Receptor-ligand interaction analysis

To identify potential interactions between and within villi and decidua, we used a list of ligand-receptor pairs, which was obtained and manually curated from the Database of Interacting Proteins (http://dip.doe-mbi.ucla.edu) and the IUPHAR/BPS guide to pharmacology (www.guidetopharmacology.org), as described earlier (30). The list of interactions was intersected with TPM-normalized expression table (data file S1) of genes for each cell type, and only ligands and receptors expressed by >15% cells and with expression >5 and >1.5 TPM, respectively, were used. Receptor-ligand pairs where both partners pass the filtering criteria were considered for graphical representation. Custom Perl and R scripts were used for analyses and to draw the interaction diagrams.

SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/4/10/eaau4788/DC1

Fig. S1. Experimental workflow and cross-platform data integration.

Fig. S2. Comparison of gene expression between 10x and Drop-seq datasets.

Fig. S3. Contribution of individual samples to each cluster.

Fig. S4. Validation of anti-REN antibody.

Fig. S5. Immunohistochemical staining data and lineage markers expression in placenta.

Fig. S6. Immunohistochemical staining data and lineage marker expression in decidua.

Fig. S7. Differential gene expression to identify cell type–specific gene signatures.

Fig. S8. Expression pattern of collagen genes over pseudotime trajectory of stromal cells (FB1, FB2, and DSC) of decidua sample P6D_10x.

Fig. S9. Differentially expressed genes between EVT and VCT based on Wilcoxon rank sum test.

Table S1. Sample information of villi and decidua data.

Data file S1. Average expression profiles of all cell types of placenta.

Data file S2. Expression of ligands and receptors.

Data file S3. Differential gene expression analysis by Wilcoxon rank sum test.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.

REFERENCES AND NOTES

Acknowledgments: We acknowledge expert support by the members of the Rockefeller University Genomics Resource Center and the Epigenomics Core facility of Weil Cornell Medicine. We also acknowledge K. Manova and M. Turkekul (Molecular Cytology Core, core grant P30 CA008748, MSKCC, New York) for their expert support on immunostaining experiments. We thank B. Rosenberg and S. Uhl (Icahn School of Medicine at Mount Sinai, New York) and E. Macosko, M. Doldman, and S. McCarroll for advice to establish Drop-seq. We also thank E. Der (Albert Einstein College of Medicine, New York) for discussions on clustering analysis. We also thank Tuschl lab members H. Totary-Jain (University of South Florida) and R. Clancy (NYU School of Medicine) for their feedback on the manuscript. Funding: H.S., P.M., A.S., T.T., and Z.W. were supported by NIH grants (R01HD086327 and U19CA179564). T.T. was supported by the Howard Hughes Medical Institute. Author contributions: T.T., Z.W., and H.S. conceived the study. H.S. and A.S. performed all sample dissociations. H.S. performed single-cell experiments with assistance of A.S. A.G. and H.S. performed the target protein overexpression and antibody validation using Western blots for downstream immunohistochemical staining. N.S., Z.W., and A.S. assisted with patient consent and sample acquisition. H.S. and P.M. performed the computational analysis assisted by M.K. K.E.A.M. standardized and performed the RNA isolation for bulk RNA-seq. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are available under BioProject ID PRJNA492324 or present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.
View Abstract

Navigate This Article