Dissecting human trophoblast cell transcriptional heterogeneity in preeclampsia using single-cell RNA sequencing.

OBJECTIVE: PE is a pregnancy-specific syndrome that affects 3%-5% of pregnant women. It often presents as new-onset hypertension and proteinuria during the third trimester. PE progresses rapidly and may lead to serious complications, including the death of both mother and fetus. In low-income countries, PE is one of the main causes of maternal and child mortality. While the cause of PE is still debated, clinical and pathological studies suggest that the placenta plays an important role in the pathogenesis of PE.
MATERIALS AND METHODS: In this single-cell RNA-sequencing (RNA-seq) study, the placenta was taken from the designated position after cesarean section. We compared placental cell subsets and their transcriptional heterogeneity between preeclampsia and healthy pregnancies using the single-cell RNA-seq technology. A developmental trajectory of human trophoblasts was shown.
RESULTS: Gene expression in endoplasmic reticulum signaling pathways in syncytiotrophoblast was upregulated in the PE group. The villi cytotrophoblasts (VCT) and extravillous trophoblasts were mainly involved in immune responses.
CONCLUSION: The placental immune function of patients with PE was altered. Proteasomes, spliceosomes, ribosomes, and mitochondria were abnormally active in the new VCT cell type.

INTRODUCTION

The placenta plays a crucial role in anchoring the conceptus, preventing its rejection by the maternal immune system, and delivering effective nourishment for the fetus, including gaseous exchange and transporting nutrients and waste between the fetal and maternal circulatory systems (Burton & Fowden, 2015; Redman & Sargent, 2005; Sibai et al., 2005). The placenta has three main types of epithelial trophoblast cells: the villi cytotrophoblasts (VCTs), the syncytiotrophoblasts (SCTs), and extravillous trophoblasts (EVTs) (Redman & Sargent, 2005; Turco & Moffett, 2019). The monocytic VCTs are located on the basement membrane beneath the SCTs. Unlike other trophoblasts, VCTs express molecular markers related to cell proliferation, such as EGFR, Met, and some members of the Wnt signaling pathway (Jokhi et al., 1994; Leisser et al., 2006; Pijnenborg & Vercruysse, 2013), so it is considered the germinal layer of the trophoblast. SCTs envelop the whole surface area of the villous tree surrounding the pregnancy, consisting of 58 billion nuclei and having a surface area of 12–14 m2 (Napso et al., 2018). SCTs also secrete hormones to support healthy pregnancy that then enter the maternal circulation, such as human chorionic gonadotropin and placental prolactin (Evain‐Brion & Malassine, 2003; Kaufmann et al., 1977; Sonderegger et al., 2007). EVTs are derived from cytotrophoblast cell columns at the distal tips of the villi and have two differentiation pathways (Pijnenborg et al., 1980). Substances released from the placenta into the maternal blood are associated with the clinical symptoms of PE, such as elevated blood pressure and levels of urinary proteins (Mastorakos & Ilias, 2003). When under oxidative stress, SCT release a series of complex factors, including proinflammatory cytokines, exosomes, and antiangiogenic agents, into the maternal circulation. These factors interfere with maternal endothelial function and lead to a systemic inflammatory response (Redman et al., 2014; Velegrakis et al., 2017).

Many studies have examined human placental transcriptome. Tsang et al. identified diverse cellular subtypes in the human placenta and enabled the reconstruction of the trophoblast differentiation trajectory by analyzing over 24,000 nonmarker selected cells from preeclamptic placentas (Tsang et al., 2017). Liu et al. performed single‐cell RNA sequencing on sorted placental cells from first‐and second‐trimester human placentas. New subtypes of trophoblasts, Hofbauer cells, and mesenchymal stromal cells were identified (Liu et al., 2018; Pavlicev et al., 2017).

However, it is still unclear what changes have taken place in gene expression profiles of various cell types in placental tissue under the condition of elevated blood pressure. In this study, we used microfluidic single‐cell digital transcriptomic technology to comprehensively characterize the transcriptomic heterogeneity of the human placenta.

MATERIALS AND METHODS

Sample collection

In this single‐cell RNA‐seqencing (RNA‐seq) study, samples were collected from three pregnant women complicated by PE and three healthy pregnant women from the Weifang Traditional Chinese Hospital, Shandong, China between December 2019 and April 2020. An overview of their clinical data is shown in Table 1. PE was defined as blood pressure ≥140/90 mmHg on at least two occasions, 4 h apart, and developing after a 20‐wk gestation with proteinuria of ≥300 mg in 24 h (Brown et al., 2018). The exclusion criteria consisted of multiple pregnancies, chronic inflammatory, autoimmune, hemolytic and renal diseases, hepatitis, placental ablation, pregnancies with chromosomal aberrations, and maternal diabetes. The study was approved by the Weifang Traditional Chinese Hospital Research and Ethics Committees. Informed consent was obtained from all participants before collection of placental samples. For placental parenchymal biopsy, 1 cm3 of placental tissue was dissected immediately after delivery from a region 2 cm deep and 5 cm away from the umbilical cord insertion point after peeling off the fetal membranes. This study was authorized by the Ethic committe of Weifang Traditional Chinese Hospital.

TABLE 1

Clinical characteristics of patients

Sample	Group	Gestation age, week	Mode of delivery (VD:CS)	Birth weight (g)
2019P1	PE	34+5	CS	3100
2019P2	PE	35+3	CS	2950
2019P3	PE	35+1	CS	3020
2019N1	HP	38	CS	3500
2019N2	HP	38+2	CS	3750
2019N3	HP	38+5	CS	3620

Abbreviations: CS, cesarean section; HP, healthy pregnancy; PE, preeclampsia; VD, vaginal delivery.

Tissue dissociation and preparation

GEXSCOPE™ Tissue Preservation Solution (Singleron) was used to store fresh placental tissue. The specimens were washed thrice with Hanks balanced salt solution and cut into 1–2 mm pieces. Tissue pieces were then digested for 15 min with 2 ml of GEXSCOPE™ Tissue dissociation solution (Singleron) at 37°C in a 15 ml centrifuge tube with uninterrupted agitation. After digestion, samples were filtered using 40‐micron sterile strainers and centrifuged at 350 g for 5 min. Then, the supernatant solution was removed and the sediment was resuspended in 1 ml of PBS (HyClone).

Single‐cell RNA sequencing (scRNA‐seq)

Single‐cell suspensions with 1 × 105 cells/ml in PBS (HyClone) were prepared. The suspensions were then loaded onto microfluidic devices and scRNA‐seq libraries were constructed according to the Singleron GEXSCOPETM protocol using the GEXSCOPE™ Single‐Cell RNA Library Kit (Singleron) (Dura et al., 2019). Individual libraries were diluted to 4 nM and pooled for sequencing. Pools were sequenced on an Illumina HiSeq X with 150 bp paired end reads.

scRNA‐seq quantifications and statistical analysis

Raw reads were processed to generate gene expression profiles using an internal pipeline. Briefly, after filtering read one without poly T tails, cell barcodes and unique molecular identifiers (UMI) were extracted. Adapters and poly A tails were trimmed (fastp V1) before aligning read two to GRCh38 using Ensemble v.92 gene annotation (fastp 2.5.3a and featureCounts 1.6.2). Reads with the same cell barcode, UMI, and gene were grouped together to calculate the number of UMIs per gene per cell. The UMI count tables of each cellular barcode were used for further analysis. Cell type identification and clustering analysis was performed using the Seurat program. The Seurat program (http://satijalab.org/seurat/, R package, v.3.0.1) was also used for analysis of RNA‐sequencing data. UMI count tables were loaded into R using the read.table function. The parameter resolution was set to 0.8 for the FindClusters function for clustering analyses. The differentially expressed genes (DEGs) were generated by the function FindMarkers of Seurat with logfc.threshold 0.25 and min.pct 0.1; Genes with q_value <0.05 in differential genes list were selected and used for Gene ontology (GO) enrichment analysis by clusterProfiler (Yu et al., 2012); For (Gene Set Enrichment Analysis) GSEA analysis, the differential genes were sorted by the avg_logFC decreasingly. Based on the rearranged genes list, GSEA analysis was performed according to the KEGG database by the ‘gseapy’ python package; Monocle 2 was used for trajectory analysis. 2000 high variable genes generated by FindVariableFeatures function of Seurat were used for genes ordering in trajectory analysis. DDRTree was used for dimensionality reduction of the dataset. And differential genes of pseudotime were screened out by differentialGenesTest function in monocle 2 (Qiu et al., 2017).

RESULTS

Single‐cell transcriptome profiles distinguished placental cell types

A total of 11,518 cells were obtained from the placental tissues of six women. A graph‐based clustering method was used to identify 13 major cell clusters by uniform manifold approximation and projection (UMAP; Figure 1a). In the PE group, 6434 placental cells included immune cells (3562 cells, 55.36%), trophoblast cells (2507 cells, 38.96%), fibroblasts (157 cells, 2.4%), erythroblasts (142 cells, 2.2%), and endothelial cells (66 cells, 1.0%) (Data file S1; Figure 1b). In the healthy pregnancy group, a similar cell heterogeneity classification was retrieved (Data file S1). Large numbers of macrophages and trophoblast cells were detected in both the PE (GT01) and healthy pregnancy (XG01) groups, with natural killer cells and T cells being higher in the PE group. Trophoblast cells could be further subclassified into EVTs, VCTs, and SCTs by employing sublineage markers such as PARP1, CGA, CYP19A1, HLA‐G, and PAPPA2 (Lee et al., 2016) (Figure 1c). EVTs (1216 cells, 48.5%) were the most abundant type of trophoblast cell, followed by VCTs (1000 cells 39.9%), and SCTs (291 cells 11.6%) in the PE group. (Data file S1). VCTs express the PARP1 gene which encodes a chromatin‐associated enzyme, Poly(ADP‐ribosyl) transferase that modifies various nuclear proteins by Poly(ADP‐ribosyl)ation. CGA, which encodes for a common subunit of the four human glycoprotein hormones chorionic gonadotropin (CG), follicle‐stimulating hormone (FSH), luteinizing hormone (LH), and thyroid‐stimulating hormone (TSH), is primarily expressed by SCTs. EVTs are characterized by the expression of HLA‐G (Fiddes & Goodman, 1981). Diseases associated with HLA‐G include severe preeclampsia and asthma (Hylenius et al., 2004; Nicolae et al., 2005; Sageshima et al., 2003). As presented in Figure 1c, those clusters exhibiting high expression of CD68, APOE, and SPP1 were cells from macrophages. The specific CD37 and NKG7 genes were expressed by T cells and natural killer cells respectively. Clusters with a high expression of IGKC and JCHAIN were plasma cells that participate in the humoral immune response. Clusters with a high expression of PECAM1 were endothelial cells. Erythroblasts specifically expressed hemoglobin subunit genes, such as HBB and HBA1.

FIGURE 1

Single‐cell transcriptomic profiling and dissection of the cellular heterogeneity of the human placenta. (a) UMAP distribution of single cells from the 13 defined cell types containing both preeclampsia and normal placental tissue. (b) The proportion of cell clusters in preeclampsia (GT01) and healthy pregnancy (XG01). (c) Biaxial scatter plots showing the expression pattern of specific genes among different subgroups of placental cells

Identification of differentially expressed genes in trophoblast cells

We compared and analyzed the genes of three types of trophoblast cells in the PE and healthy pregnancy groups. SCTs revealed 610 DEGs (Data file S2). Using the average expression profiles of SCTs, VCTs, and EVTs, we performed unsupervised clustering analysis of the union of the top 40 DEGs including down‐regulated and up‐regulated genes that were expressed (Figure 2a, Figure S1a,b). Using all DEGs, we examined the main signaling pathways involved in SCTs by GSEA. Among them, genes related to protein processing in the endoplasmic reticulum pathway were upregulated in the PE group (Figure 2b). Using upregulated genes, we performed functional annotation with GO analysis. The enriched items of the biological process focused on the regulation of interleukin‐8 production and interleukin‐8 production (Figure 2c). A total of 347 DEGs were identified in VCTs (Data file S3). All these genes were analyzed by GSEA analysis. The expression of genes related to the ribosome pathway was downregulated in the PE group (Figure S1c). The GO enrichment analysis was performed on those genes whose differential gene expression was up‐regulated. The biological process was mainly enriched in the innate immune response and cell chemotaxis (Figure S1d). EVTs revealed 283 DEGs (Data file S4) and GSEA analysis was performed on all these differential genes. Genes related to the herpes simplex infection pathway were upregulated in the PE group (Figure S1e). The GO terms enriched in these cells were associated with the innate immune response and response to cytokines (Figure S1f).

FIGURE 2

SCT differentially expressed genes and functional changes between the hypertensive group and the healthy pregnancy group. (a) Heat map showing the TOP 40 DEGs, including downregulated and upregulated genes, in SCTs. (b) GSEA result confirming enrichment of the protein processing signaling pathway found in endoplasmic reticulum expression in preeclampsia. (c) GO enrichment term of upregulated genes in SCTs

Reconstruction of the developmental relationship of the three types of trophoblasts

Our dataset captured a significant number of trophoblastic cells. We observed that VCTs, SCTs, and EVTs were clustered in a continuum (Figure 1a). In fact, it has been proposed that these trophoblast subtypes are developmentally connected (Tsang et al., 2017; Vento‐Tormo et al., 2018). To delineate this relationship and to study the genes that regulate this process, we first identified the highly variable genes from these three clusters using UMAP reclustering analysis. Through the analysis of cell trajectories, we found that VCTs were in the early stage of cell development. We found that both the PE group and the healthy pregnancy group showed a continuous “Y” shaped trajectory of trophoblast development. (Figure 3a,b). FLT1, FN1, and TIMP3 were specifically upregulated along the differentiation pathway of SCTs (Figure 3c).

FIGURE 3

Reconstruction of the developmental relationship of trophoblast cells using pseudotime analysis. (a) The trophoblastic cell developmental trajectory visualization in a biaxial scatter plot. Dark colors indicate early development. (b) Distribution of trophoblast subtypes in trajectory. Trophoblast differentiation tracks are shown separately in the preeclampsia group and the healthy pregnancy group. (c) Dot plot showing the expression of the indicated genes in the indicated subtypes of cells

Identification of new subtypes of trophoblast cells

Further analysis revealed that the PARP1 positive cells were clustered into three subtypes: VCT‐1, VCT‐2, and VCT‐3 (Figure 4a). The cluster exhibiting the highest expression of CDK1, CCNB1, and RRM2 was VCT‐1. This cluster also had the highest expression of IFIT2 and IFIT3, which were the marker genes identifying VCT‐2 (Figure 4b).

FIGURE 4

Trophoblast subtypes are present. (a) Biaxial scatter plot showing single‐cell transcriptomic clustering of the VCT subtype using UMAP analysis. (b) Biaxial scatter plots showing the expression pattern of specific genes among different trophoblast cell subgroups. (c) The histogram shows the proportion of VCT subtypes in the preeclampsia group and the healthy pregnancy group. VCT‐2 is mainly distributed in the preeclampsia group, while VCT‐3 is mainly distributed in the healthy pregnancy group

VCT‐2 cells were mainly distributed in the PE group, while VCT‐3 cells were mainly distributed in the healthy pregnancy group (Figure 4c). To further understand the gene functions of VCT‐2 and VCT‐3, correlation analysis was performed. The expression of VCT‐2 was compared with other trophoblast genes, and 1151 DEGs were obtained (Data file S5). GSEA analysis of all differential genes showed that the expression of genes related to the proteasome, spliceosome, and ribosome pathways was upregulated, but the genes related to antigen processing and presentation were downregulated (Figure 5a, Figure S2a). We found that proteasome related genes, such as PSMA7 and PSMC1, were upregulated to varying degrees in the VCT‐2 differential geneset (Data file S6). The enriched GO terms in VCT‐2 were associated with cellular respiration, energy derivation by oxidation of organic compounds, and generation of precursor metabolites and energy (Figure 5b). Genes whose expression levels were upregulated included mitochondrial coding genes related to the cellular respiratory chain (Data file S5).

FIGURE 5

Functional enrichment and signal pathway analysis of new subtypes of trophoblast cells. (a) GSEA result confirming enrichment of the proteasome signaling pathway in VCT‐2. Gene sets from VCT‐2 are compared with other trophoblast cells. (b) Dot plots show the GO enrichment term of upregulated genes in VCT‐2. (c) GSEA result confirming enrichment of metabolic pathways in VCT‐3. Gene sets from VCT‐3 are compared with other trophoblast cells. (d) GO enrichment term of upregulated genes in VCT‐3

VCT‐3 was compared with other trophoblast cells to obtain 1353 DEGs (Data file S7). GSEA analysis of all the differential genes showed that the expression of genes related to metabolic pathways and the Alzheimer's disease pathway were upregulated (Figure 5c, Figure S2b, Data file S8). The expression of genes related to natural killer cell‐mediated cytotoxicity and rheumatoid arthritis were downregulated (Figure S2b). Based on GO analysis it appears from the highly enriched terms that VCT‐3 participates in the nuclear‐transcribed mRNA catabolic process, and cotranslational protein targeting to the membrane (Figure 5d). Related genes are: RPL23A and RPS14, etc. (Data file S7).

DISCUSSION

By analyzing the transcripts of 11,518 cells, the difference between the PE and healthy pregnancy groups was revealed at the single‐cell level. We obtained 1240 DEGs in total. Through further analysis of the signal pathway of SCTs, we found that the expression of protein processing in endothelial reticula was upregulated, indicating that endoplasmic reticulum function was active in the diseased state. This was thought to be related to the intermittent hypoxia of the placenta caused by narrow spiral arteries. During hypoxia, protein folding by the endoplasmic reticula is suspended, resulting in a cessation of cell proliferation and subsequent severe apoptosis (Huppertz et al., 2006). Our findings were consistent with previous studies. The release of apoptotic trophoblast particles into the maternal blood will stimulate inflammation in blood vessels (Redman & Sargent, 2007). The Go enrichment analysis of VCTs and EVTs showed that the cellular function of the PE group focused on the immune response, suggesting that the placental immune function of patients with PE was altered. On the one hand, EVTs need to remodel the uterine spiral arteries, and at the same time to invade the decidual tissue of the uterus and anchor the placenta to the uterine wall. In these processes, different types of maternal cells, such as decidual macrophages, uterine NK cells, and stromal cells, are encountered. To adapt to each other, reasonable immune regulation is necessary (Moffett et al., 2017; Pollheimer et al., 2018). Through our experiments, and for the first time, we identified three subtypes of VCTs. After GSEA analysis of VCT‐2, we found that the proteasomes, spliceosomes, ribosomes, and mitochondria were abnormally active in this cell. Proteasomes are a type of protein complex; their main function is to degrade proteins that are not needed or that have been damaged in cells through the breaking of peptide bonds. The proteasome degradation pathway is essential for many cell processes, including the cell cycle, gene expression regulation, and oxidative stress response, etc. (Shang & Taylor, 2011). Plasma proteasome and immunoproteasome are significantly upregulated in PE and hemolysis, elevated liver enzymes, and thrombocytopenia syndrome (Berryman et al., 2019). Sun et al. analyzed the gene expression profile of peripheral blood of patients with PE using gene microarray technology, and found that the expression of PSMC5 was increased (Sun et al., 2008). Our findings confirmed many previous conclusions on PE, such as stress to endoplasmic reticula and mitochondrial dysfunction. Different from previous transcriptome studies, our results were accurate to the single‐cell level.

There are some limitations of the current study. First, only three pregnant women were reruited in each group. Our results should be validated in large‐scale studies. Secondly, the interactions among these active signaling pathways in PE warrant further investigation.

CONCLUSION

This study has enriched the understanding of human placenta and its molecular functions in the diseased state. These findings will contribute to the future discovery of new molecular markers and the discovery of drugs for the treatment of the disease.

CONFLICT OF INTEREST

The authors have no conflicts of interest to declare.

AUTHOR CONTRIBUTIONS

T Z, H Z, QQ B, and YC C were involved in study design. SH L, SW Y, P J, L L, and SJ D prepared the samples for the experiments. T Z and XL W carried out data analysis. QQ B and M S drafted the manuscript. All the authors contributed to critical discussions and approved the final version of the manuscript.

ETHICS STATEMENT

This study was authorized by the Ethic committee of Weifang Traditional Chinese Hospital.

Fig S1‐S2

Click here for additional data file.

Data S1‐S8

Click here for additional data file.