Single-cell transcriptomic profiling reveals diverse lineages of day 9–11 human embryos developed in vitro
To lay the foundation for ploidy-specific transcriptomic analysis in human embryos beyond implantation, we generated a global transcriptomic profile from human embryos undergoing development in vitro (Fig. 1a). Specifically, we included 20 embryos with PGT-A (euploid: n = 6, aneuploid: n = 10, mosaic: n = 4), and 13 embryos without PGT-A (Supplementary Table 1). We first cultured preimplantation embryos from 5–6 d.p.f. to 9–11 d.p.f. utilizing previously reported extended embryo culture protocols18,19 (Supplementary Video 1). At the end of the culture, 25 embryos (5 euploid, 7 aneuploid, 2 mosaic and 11 untested embryos) in good condition were processed for scRNA-seq using the 10x Genomics platform. After quality control filtering, 3159 cells from 20 embryos were included for downstream bioinformatic analyses (Supplementary Table 2, Supplementary Fig. 1a-b).
To visualize the development and segregation outcome in these post-implantation embryos, immunofluorescence imaging was performed at the end of the culture and confirmed the presence of three major lineages. These included: OCT4-positive epiblast cells, GATA4-positive hypoblast cells, and CK7-positive trophoblast cells (Fig. 1b, Supplementary Video 2). Single-cell transcriptomic data identified cells from five primary lineages, both embryonic and extraembryonic, that clustered based on their specific gene expression profile (Fig. 1c). This included the epiblast, the hypoblast and three different trophoblast subgroups: cytotrophoblast, syncytiotrophoblast and extravillous trophoblast. Each lineage was identified from differential gene expression analysis, which correspondingly showed enrichment of known lineage-specific markers (Fig. 1d-e, Supplementary Fig. 2). The epiblast was marked by the expression of POU5F1, DPPA5 and KHDC3L, the hypoblast by APOA1, MARCKS, FN1 and S100A4 expression, the cytotrophoblast by expression of KRT19, PEG10, TPM1 and FABP5, the syncytiotrophoblast by expression of CGB subunit genes, ERVW-1 and HOPX, and the extravillous trophoblast by expression of MMP2, EPAS1 and ITGA1. The lineage-specific expression patterns in our data were concordant with previously reported gene expression profiles in post-implantation human embryos20–25. To confirm that the phenotype of our embryos was not unique to our lab, we compared and integrated our single cell data with that of Molè et al., which comprises 16 embryos (8 at 9 d.p.f. and 8 at 11 d.p.f.) and a total of 4820 cells22. Our embryos contributed to all major embryonic lineages previously reported (Supplementary Fig. 1c-e). These results confirmed the applicability of an embryo extended culture approach combined with single-cell transcriptomic analysis as a model system to investigate post-implantation lineage differentiation in human embryos.
Copy number variation (CNV) and aneuploid cell inference from embryo gene expression data
To delineate the potential impact of aneuploidy on the developmental transcriptome in post-implantation embryos, we inferred cell CNV from the gene expression data using R inferCNV26,27 (Fig. 2a), which showed both common and embryo-specific aberrations in chromosomes 7, 9, 12, 19, 20, 21, and 22. This approach has been validated by G&T-seq (genome and transcriptome sequencing)28 to reliably perform copy number calls in human embryos, particularly after embryonic genome activation7. For each cell, we obtained a preliminary copy number aberration score by summing the deviation of expression intensity of each gene from its average expression intensity across all cells (Fig. 2b, Supplementary Fig. 3a). We then selected 10 cells in the lower 5th percentile of aberration score (Supplementary Table 3, Supplementary Fig. 3b) to serve as the euploid reference, spanning across all embryo lineages (Fig. 2c). This allowed gene- and chromosome-level CNV state prediction in the other cells and thus aneuploid cell inference (Supplementary Fig. 3c). We called cells aneuploid if they contained at least one chromosome exceeding 70% CNV regions (Fig. 2d). In total, 1524 out of 3159 (48%) cells were annotated to be aneuploid (Fig. 2e, Supplementary Fig. 3d). We then categorized the embryos as euploid or aneuploid if they contained > 90% euploid or aneuploid cells respectively, and the rest as mosaic embryos. Our data was categorized into 11 euploid embryos, 6 aneuploid embryos and 3 diploid-aneuploid mosaic embryos (Fig. 2f, Supplementary Fig. 3e-g).
Ploidy-specific lineage dynamics beyond implantation
To explore lineage signatures from the transcriptomic data, we analyzed lineage-specific cell proportion in the three embryo groups (Fig. 3a, Supplementary Fig. 3h, Supplementary Table 4). In particular, aneuploid embryos showed more pronounced inter-embryo variations in their lineage-specific cell proportion than euploid and mosaic embryos, especially in the epiblast. Moreover, 80% of the aneuploid embryos had abnormal lineage composition - either completely lacking the epiblast or trophoblast lineages (Fig. 3b, Supplementary Fig. 3i)2,17.
Mosaic embryos contained cells from across the lineages, irrespective of ploidy except for the hypoblast (Fig. 3c-d). In the 4 extraembryonic compartments (e.g. HYPO, CTB, EVT, STB), the average proportion of aneuploid cells (73%) was higher than euploid cells (27%) (Fig. 3c). However, in the epiblast, 40% of the cells were aneuploid while 60% of the cells were euploid. We note, however, drawing conclusions on the enrichment of aneuploid cells for specific cell types must be validated with a larger sample size of mosaic embryos in extended culture.
Differential expression (DE) analyses revealed distinct patterns of mosaic and whole-embryo aneuploidy
Based on unsupervised clustering results, the primary driver of transcriptomic diversity was embryo lineage. We next aimed to unravel the global effect of aneuploidy on the transcriptome of mosaic and aneuploid embryos to better understand aneuploid cellular functional states that could affect post-implantation outcomes such as embryo viability2,16,18. Therefore, we performed differential expression (DE) analysis on different embryo cell populations using Pseudobulk comparisons, which effectively reduces false positives in single cell DE analysis29. Differentially expressed genes (DEGs) were identified using four comparison groups to evaluate ploidy-specific transcriptomic patterns: 1) aneuploid cells of aneuploid embryos (AA) versus euploid cells of euploid embryos (EE) - AA vs. EE; 2) aneuploid cells of mosaic embryos (AM) versus euploid cells of mosaic embryos (EM) - AM vs. EM; 3) aneuploid cells of mosaic embryos (AM) versus aneuploid cells of aneuploid embryos (AA) - AM vs. AA; 4) euploid cells of mosaic embryos (EM) versus euploid cells of euploid embryos (EE) – EM vs. EE (Fig. 4a-d). In addition, we noted the differences in lineage composition within two compared populations and tested for DEGs with and without normalizing by lineage amount. From the unnormalized run, markers specific to the dominant lineages were present in DEGs. Thus, in order to reduce lineage effect on DE analysis and highlight ploidy-specific transcriptomic patterns, we balanced the number of cells in each lineage of the comparison groups by downsampling (Supplementary Table 5). From the normalized run, most but not all lineage markers disappear, while the other DEGs remained largely unaltered in the following analyses, supporting our findings using a pseudobulk approach (Supplementary Table 6–9).
When comparing AA vs. EE, there were 229 DEGs which included lineage-specific markers such as GATA3, APOE and MT1G (Fig. 4a). Intriguingly, there were 5 times more DEGs when comparing AM vs. EM (Fig. 4b). This could potentially result from a higher population number in AA vs. EE than AM vs. EM, or alternatively, the larger difference suggests that cells in mosaic embryos are phenotypically different than their homogenous counterparts. The presence of two genetically distinct populations in mosaic embryos potentially altered the cellular sensitivity towards gene dosage. The average logFC value was positive in the DEGs of AM vs. AA but negative in the other three comparison groups (Fig. 4a-d), indicating more genes may be activated in aneuploid cells contained within a diploid-aneuploid mosaic environment than aneuploid cells present in a homogenous aneuploid environment.
DEGs from AA vs. EE provided insight into the global transcriptomic changes induced by aneuploidy (Fig. 4a, Supplementary Table 6). Genes highly upregulated in AA included DDIT4, a marker of response to DNA stress30. The most downregulated genes in AA were genes encoding the RPS28, which is involved in cellular proliferation31. To elucidate the molecular pathways and functions related to the DEGs, we performed gene set enrichment analysis (GSEA). We found that the enriched genes in AA were involved in endoplasmic reticulum (ER) stress responses, likely due to the changes in gene dosage. This included elevated unfolded protein response (UPR, e.g. BCL2L1132), which is known to initiate degradation of the misfolded protein; and cellular response to interleukin-1 (e.g. NLRP733), which is a defensive response during inflammation. Besides, dysregulated cell cycle (e.g. CCND134) and upregulated glycolytic process (e.g. HIF1A35) were obserevd in AA. We also observed downregulation of cellular proliferation mechanisms, including DNA replication (e.g. MCM536) and oxidative phosphorylation (e.g. ATP5ME37).
Next, we sought to scrutinize whether euploid-aneuploid coexistence led to aneuploidy-induced transcriptomic responses in the mosaic embryos. To this end, we performed DE analysis in AM vs. EM (Fig. 4b, Supplementary Table 7). Similar to AA, the highly upregulated genes in AM pointed to the presence of anti-proliferative signals (e.g. TP6338). However, some markers highly upregulated in AA became downregulated in AM, including DDIT4 and HIF1A (Supplementary Fig. 4a-b). GSEA also revealed an increase in complement activation (e.g. PHB239) and leukocyte activation (e.g. CLU40), corresponding to ongoing inflammation processes. Surprisingly, gene sets involved in replication surveillance mechanisms were downregulated in AM (Supplementary Fig. 4c), such as DNA damage repair (e.g. DUT41), and cell cycle transition regulation (e.g. EIF4EBP142). Similarly, correction functions that act to resolve ER stress-induced protein burden were downregulated (Supplementary Fig. 4d)., including UPR (e.g. ATF432) and ubiquitin-mediated proteolysis (e.g. UBE2J143). Downregulation of ULK1 and TMBIM6 was also detected. Both genes act as crucial mediators to initiate autophagy of misfolded protein through the UPR pathways51,76,77. Overall, these findings suggests that functions to maintain proteostasis and cell fidelity may be impaired in AM (Fig. 5). Furthermore, in contrast to the increased expression of genes related to glycolysis and metabolism in AA, glycolytic processes (e.g. HIF1A) were downregulated in AM. Metabolic elevation is a hallmark of aneuploidy across different models to cope with increased energy expenditure and immune activity44. Over time, this attempt to regain cellular fitness may contribute to aneuploidy tolerance44,45. In contrast, AM demonstrated distinct metabolic patterns, which may compromise their proliferation and survival in the diploid-aneuploid microenvironment. Lastly, we did not detect significant positive enrichment of apoptotic pathways in the aneuploid cells of mosaic embryos from GSEA despite the presence of pre-apoptotic transcriptomic signatures (Supplementary Fig. 4e).
To further explore whether the transcriptomic signatures in AM population were impacted by the neighboring EM population, we performed DE analysis between AM and AA (Fig. 4c, Supplementary Table 8). Genes upregulated in AM demonstrated transcriptomic patterns related to genomic imbalances, including upregulated apoptotic cell clearance (e.g. TGM246) and negative regulation of growth (e.g. SEMA4D47). On the other hand, the downregulated genes were involved in energy metabolism and cell cycle integrity, indicated by decreased glycolysis (e.g. HIF1A) and dysregulated G1/S transition (e.g. CCND1). Further, UPR (e.g. ATF4) and autophagy (e.g. CHMP2A48), which are essential mechanisms to protect cells against stress-induced apoptosis, were both downregulated in AM. These results indicate that the presence of euploid cells may further contribute to the survival disadvantage of aneuploid cells in a chimeric environment. Next, we aksed whether the presence of aneuploid cells affects the gene expression in the neighboring euploid cells by performing DE analysis between the EM and EE (Fig. 4d, Supplementary Table 9). We found enriched gene sets related to elevated TNF-alpha signaling (e.g CEBPD49) and decreased translational activity (e.g. RPS2031). These patterns indicate the fitness of EM was somewhat affected by the environmental stress due to genetic instability of the coexisting aneuploid cells. Additionally, elevated protein catabolic activity (e.g. PSMD950) and autophagy activity (e.g. ULK151) were observed in EM. This suggests the presence of aneuploid cells led to the activation of recovery mechanisms in the neighboring euploid cells of mosaic embryos52, which we hypothesized may be an adaptive response to environmental stress.
We next examined the overlapped DEGs between AA vs. EE, AM vs. EM and AM vs. AA in order to identify common transcriptomic characteristics of aneuploid cells (Fig. 6). There was no overlap in upregulated DEG between all three groups. On the other hand, C5orf66-AS1 was downregulated in all three groups, a gene essential to promote cell proliferation53. Common DEGs were the most abundant between AM vs. EM and AM vs. AA. Common upregulated DEGs included TP63 and TGM654, which were involved in response to environmental inflammation and ER stress. DIO2, which converts thyroxine (T4) to triiodothyronine (T3), is crucial for stem cell maintenance and TE lineage differentiation55. The common downregulated DEGs included DDIT4 and DUT, which are involved in key recovery mechanisms in response to cell damage. The pattern of common DEGs between AA vs. EE and AM vs. EM was similar to the common DEGs between AA vs. EE and AM vs. AA. Common upregulated DEGs were involved in the inflammatory process. Common downregulated DEGs were related to cell cycle checkpoint, as well as cell growth and survival.
Ligand-receptor interactions suggest euploid cells participate in robust cell-stress signaling in mosaic embryos
When cell competition happens, cell-cell interactions between two cell populations are essential for fitness sensing, leading to different survival outcomes13. To understand the global crosstalk between the euploid and aneuploid populations of mosaic embryos, we performed ligand-receptor (L-R) reactome analysis between EM and AM using the CellChat56. AM were involved in a higher number of L-R interactions than the euploid cells, suggesting genomic imbalance can trigger complex cellular responses (Fig. 7a). However, the overall interaction strength was significantly higher in EM, indicating consistent and effective communication with the surrounding environment. In total, 14 enriched L-R interaction pairs and 11 enriched pathways were identified between the EM and AM (Fig. 7b-c). Pathway analysis of the L-R pairs revealed enrichment in stress response pathways such as TGF-beta signaling, p53 pathway and hypoxia. We observed similarity in the L-R profiles involving AM, whereas the L-R profile between EM (Euploid -> Euploid) showed a distinct pattern. L-R interactions involving TGFB3 were highly expressed in the L-R profiles involving AM. As a key regulator of cell death, TGFβ induces the expression of the pro-apoptotic protein Bim57. Besides, GDF15-TGFBR2 and CDH5-CDH5 were two other dominant interactions in AM. Both interactions are known to be elevated in response to inflammation and oxidative stress58,59. Under oxidative stress, cell growth is inhibited, and cells can arrest cell cycles at the G1/S interface. Indeed, DE analysis between AM vs. EM also revealed downregulation of CCNE1, CDK4 and CDK6 (Fig. 4b), which are key players in the G1/S checkpoint to promote S phase entry60. Common downregulated DEGs between AM vs. EM and AM vs. AA were also heavily involved in cell cycle regulation (Fig. 6). IGF2-(ITGA6 + ITGA4) interaction was exclusively present between AM (source) and EM (target). Such interaction has been shown to reduced apoptosis and increase mitosis in the targeted cells61. On the other hand, MPZL1-MPZL1 and F11R-F11R were exclusively enriched between EM (source) and EM (target). Both interactions are key mediators of cell proliferation and migration62,63.