Major transcriptomic, epigenetic and metabolic changes underlie the pluripotency continuum in rabbit preimplantation embryos

In mammalian embryos, totipotent blastomeres become pluripotent after differentiation of the trophectoderm lineage during the morula to blastocyst transition, and form a seemingly coherent cluster of cells called the inner cell mass (ICM). When some ICM cells differentiate to the primitive endoderm (PE), pluripotency is confined to the epiblast (EPI), which will give rise to the embryo proper. Primordial germ cells (PGCs) are then specified in the epiblast. Ultimately, the posterior part of the epiblast gives rise to the primitive streak (PS) (Blakeley et al., 2015; Chazaud and Yamanaka, 2016; Hayashi et al., 2007; Nakamura et al., 2016; Petropoulos et al., 2016; Stirparo et al., 2018). The developmental window that extends from morula differentiation to gastrulation lasts 2-3 days in mice, compared with 7-10 days in humans and non-human primates (Nakamura et al., 2021; Shahbazi, 2020). Within that window, pluripotent cells are thought to transit through three main pluripotency states known as the naive, formative and primed states (Savatier et al., 2017; Smith, 2017; Takahashi et al., 2017). In mice, the ICM and epiblast cells of blastocysts, as well as the embryonic stem (ES) cell lines that are derived from them, are in the naive state. Naive pluripotency is associated with the expression of transcription factors Klf2, Klf4, Gbx2, Tfcp2l1, Esrrb, Tbx3 and Sall4 (Dunn et al., 2014), low DNA methylation (Leitch et al., 2013, 2016), high mitochondrial respiration (Carbognin et al., 2016; Sone et al., 2017) and enrichment of active histone modifications at promoter regions of developmental genes (Hayashi et al., 2008). Formative pluripotency characterizes the epiblast cells of peri-implantation embryos, as well as the corresponding pluripotent stem cell lines, called formative stem (FS) cells (Kinoshita et al., 2021; Smith, 2017). The formative state is associated with an increased expression of transcription factors Etv5, Rbpj, Tcf3 and Otx2, and intermediate DNA methylation (Kalkan et al., 2019; Kinoshita et al., 2021). Finally, the epiblast cells of early to late gastrula-stage embryos and their in vitro counterparts, epiblast stem cells (EpiSCs), are in the primed pluripotent state (Brons et al., 2007; Tesar et al., 2007). The primed state is associated with an increased expression of Sox3 and Oct6 transcriptional regulators (Corsinotti et al., 2017), high DNA methylation (Habibi et al., 2013; Hackett et al., 2013), low mitochondrial respiration (Zhou et al., 2012) and enrichment of repressive histone modifications at promoter regions of developmental genes (Hayashi et al., 2008; Smith, 2017). These three successive pluripotent states characterize the pluripotency continuum.

The characterization of the mammalian embryo transcriptome in a wide range of species is essential for unravelling the molecular mechanisms involved in pluripotency and understanding the adaptation of embryos in maintaining pluripotency across different developmental strategies. It is also important for the development of chemically defined culture media, which can aid in the preservation and maintenance of naive, formative and primed pluripotency in embryo-derived stem cell lines. Single-cell transcriptomic data from preimplantation embryos is available in a large variety of species, including mice (Argelaguet et al., 2019; Mohammed et al., 2017; Tang et al., 2010), pigs (Kong et al., 2020), cattle (Zhao et al., 2016), marmosets (Boroviak et al., 2018), rhesus and cynomolgus macaques (Liu et al., 2018; Nakamura et al., 2016), and humans (Blakeley et al., 2015; Petropoulos et al., 2016). Transcriptome data at pre- and mid-gastrula stages is available only in mice (Peng et al., 2016; Wen et al., 2017) and cynomolgus macaques (Nakamura et al., 2016). It is indeed difficult to study this developmental period, as embryos have already implanted by the gastrula stage in many species. In mice, humans and non-human primates, gastrulation starts 2-8 days after implantation, when embryos are already deeply buried in the uterine wall, making the formative- and primed-pluripotency states difficult to study (Nakamura et al., 2021; Shahbazi, 2020). In contrast, rabbit embryos do not implant until the end of day 7 of development, after the onset of gastrulation (Fischer et al., 2012; Puschel et al., 2010). This strategy of lagomorph development allows easier access to a larger window of development compared with rodents and primates. However, the transcriptome of the rabbit embryo is poorly documented (Leandri et al., 2009; Schmaltz-Panneau et al., 2014). For this reason, we conducted a transcriptome analysis of rabbit preimplantation embryos between E2.7 (morula stage) and E6.6 (early primitive streak stage) using two complementary techniques – bulk RNA-sequencing (Illumina NextSeq) and single-cell RNAseq (10x Genomics) – to study specific gene expression. In addition, we studied oxidative phosphorylation and glycolysis, and analysed epigenetic modifications. With this unique combination of tools and cellular processes, we explored in detail the major alterations that occur in pluripotent cells in vivo, along the pluripotency continuum of rabbit embryos.

To obtain a first transcriptional map of the rabbit preimplantation embryo development, we first performed high-depth bulk RNA-sequencing of micro-dissected embryos between embryonic day E2.7 (morula) and E6.6 (expanded blastocyst at the early primitive streak-stage) (Fig. 1A, Fig. S1A). The 51 samples – morula, trophectoderm (TE), inner cell mass (ICM), epiblast at E6.0 (EPI), anterior epiblast at E6.3/6.6 (EPIant) and primitive endoderm (PE) – resulting from these dissections were separated according to developmental stages in a principal component analysis (PCA) (Fig. 1B). We then investigated the expression of some known lineage markers of mouse, human and cynomolgus macaque embryos (Table S1), and represented the results on a heatmap (Fig. 1C). As expected, the trophectoderm (TE; E3.5 to E6.6) and inner cell mass (ICM; E3.5 and E4.0) were enriched in trophectoderm (e.g. TFAP2C and KRT18) and pluripotency [e.g. DPPA5, ESRRB, POU5F1 (OCT4) and SOX2] transcripts, respectively. Most of the naive pluripotency genes had lower transcript levels in epiblast samples (EPI at E6.0, E6.3 and E6.6, thereafter called EPI_6.0, EPIant_6.3 and EPIant_6.6, respectively), which showed higher expression of late epiblast markers (e.g. OTX2). Finally, the primitive endoderm samples (PE at E6.0, thereafter called PE_6.0) were enriched in PE-specific transcripts such as SOX17 and GATA6. PE markers, including PDGFRA, GATA6, HNF1B, FOXA2, RSPO3 and SOX17 were detected in the E3.5 and E4.0 ICM samples, suggesting that the PE differentiation begins as early as E3.5. Note that FBXO15 and TDGF1/CRIPTO, two markers of naive pluripotency in mice (Boroviak et al., 2015; Mohammed et al., 2017), were upregulated in rabbit EPI samples. Conversely, the TE markers PLAC8 or SLC7A4 expressed in humans (Blakeley et al., 2015) were not found in rabbit TE. On the other hand, the primate-specific EPI markers ZIC3 and FGF2 (Nakamura et al., 2016), as well as the mouse-specific EPI markers FGF5 and SOX4 (Mohammed et al., 2017), were all expressed in rabbit EPI samples. Overall, these observations suggest that rabbit gene expression patterns share characteristics with both rodents and primates, although there are some notable differences.

We next used the 10x Genomics technology to refine the transcriptional map of the rabbit preimplantation embryo development at a single-cell resolution (Fig. 1A). This technology is best suited for analysing large numbers of cells, although at a lower sequencing depth than bulk RNAseq (Nowotschin et al., 2019). Note that our single-cell analysis included additional E5.0 embryos that could not be included in the bulk RNAseq study. After quality control, we retained 13,942 single-cell transcriptomes from 529 embryos, ranging from E3.0 (morula) to E6.6 (early primitive streak-stage) (Fig. S1B-D). We applied uniform manifold approximation and projection (UMAP) dimension reduction and found that cells separated according to developmental stage when projected onto the first two components (Fig. 1D, Fig. S1E-G, Table S2). The six embryonic day clusters separated further into smaller clusters as follows: E3.0 cells (morula stage) formed one cluster. Cells harvested from the E3.5 early blastocysts formed two distinct clusters. The first one was enriched in ICM/EPI genes (e.g. SOX2) whereas the other cluster was enriched in TE genes (e.g. GATA2). At E4.0, cells separated into three clusters: the presumptive ICM/EPI and TE clusters, and a third cluster corresponding to presumptive PE represented by PDGFRA expression. At E4.0, the presumptive EPI and PE clusters were still close to each other. In contrast, at E5.0, the three lineages formed well-separated clusters. Based on these observations, we were able to delineate with certainty the morula, ICM, EPI, TE and PE clusters on the UMAPs.

To have a more global view of differential gene expression between the main cell clusters, a heatmap was generated using five landmark markers (Fig. 1E, genes in colour fonts), with the 20 most overexpressed genes in TE versus all other cell types, ICM/EPI versus all other cell types, and PE versus all other cell types (Fig. 1E,- Table S2). Most of the genes represented on the heatmap had similar expression profiles to those described in mouse, cynomolgus macaque and human (Guo et al., 2010; Mohammed et al., 2017; Nakamura et al., 2016; Stirparo et al., 2018). However, two genes showed unexpected expression profile. Namely, SUSD2, a marker of ICM in primates (Bredenkamp et al., 2019), is highly expressed in rabbit TE, and CCND2, a marker of the primitive streak in mice (Wianny et al., 1998), is already expressed in E5.0 EPI cells in rabbits. Finally, we integrated our dataset with published single-cell data from mouse and cynomolgus macaque (Nakamura et al., 2016) and thus confirmed the identity of the rabbit clusters (Fig. S2). Taken together, these results indicate that rabbit preimplantation embryos typically express the same lineage markers as mouse and cynomolgus macaque embryos, but there are also noticeable differences.

To obtain a dynamic view of lineage specification, we performed pseudo-time analyses of the single-cell RNAseq dataset based on the 1000 most variable genes. A first pseudo-time analysis was performed with the morula, ICM/EPI_3.5/5.0 and TE 3.5-5.0 cells (Fig. 2A). The morula cells (E3.0) split into two main branches at E3.5, leading to either EPI or TE cells at E4.0 (Fig. 2B). At E5.0, EPI and TE cells formed two distinct clusters at the tip of their respective branches. We also identified two minor clusters located between the main clusters (Fig. 2A) and connected to the TE branch (Fig. 2B); we qualified those cells as ‘TE intermediates’. A heatmap (Fig. 2C) was generated grouping 10 landmark genes with the 20 most overexpressed genes in the ICM/EPI versus all other cell types and TE versus all other cell types (Table S2). Upregulated genes in TE versus EPI lineage included TFAP2C, GATA2, AQP3 and CDX2 (Fig. 2C). Immunostaining of rabbit embryos confirmed ICM/EPI-enriched expression of SOX2 at E3.5 and E4.0, and TE-enriched expression of TFAP2C at E4.0 (Fig. 2D).

In the second pseudo-time analysis, we used the same dataset, except that the TE cells were replaced by the PE 4.0-5.0 cells (Fig. 2E). This analysis revealed a clear split of embryonic cells into two lineages, EPI and PE, respectively, from E4.0 onwards (Fig. 2F). At E5.0, EPI and PE cells formed two distinct clusters at the tip of their respective branches. A heatmap (Fig. 2G) was generated grouping 10 landmark genes with the 20 most overexpressed genes in the ICM/EPI versus all other cell types and PE versus all other cell types (identified in Fig. 1E, Table S2). Upregulated genes in PE versus EPI lineage included GATA6, LAMA1, SPARC and PDGFRA (Fig. 2G). Immunostaining of rabbit embryos confirmed PE-enrichment of GATA6 versus EPI-specific localization of SOX2 and EPI-enrichment of OCT4 at E4.0 and E5.0 (Fig. 2H). Together, these results confirm the establishment of the three main lineages, EPI, PE and TE, at E5.0 in rabbit embryos.

We sought to characterize the segregation of endodermal cells into visceral and parietal endoderm (VisE and ParE, respectively) using single-cell RNAseq data from PE cells (Fig. 3A). Although OTX2 is primarily described as a formative/primed marker gene, it is re-expressed in the VisE in mice (Perea-Gomez et al., 2001). In primate embryos, it is expressed in the late EPI and primitive endoderm, and is extinguished in the ParE (Boroviak et al., 2018). The expression pattern of OTX2 can therefore be used to determine the timing of VisE/ParE segregation. At E4.0 and E5.0, expression of OTX2 gene was quite homogeneous among endodermal cells. However, at E6.0 and E6.6, endodermal cells formed two closely related clusters and OTX2 expression was restricted to one of these two subgroups (Fig. 3A). To characterize these OTX2-positive cells, we performed immunofluorescence analysis of E6.0 embryos. Strong OTX2 signal was observed only in SOX2-negative VisE cells underlying the layer of SOX2/OTX2 double-positive EPI cells and was not detected in the ParE (Fig. 3B). These results suggest that segregation of visceral and parietal endoderm takes place at E6.0 in rabbit embryos, i.e. before embryo implantation.

To determine the timing of anterior visceral endoderm formation, we next examined the expression of CER1, HHEX and LEFTY2 (Fig. 3C), and identified a sub-cluster of VisE cells co-expressing HHEX and CER1 at E6.0 (Fig. 3D). These findings point to the onset of VisE polarization at E6.0 in rabbit. Remarkably, we noticed the upregulation of BMP2 and BMP4 in the parietal endoderm (Fig. S3A), as previously described using in situ hybridization (Hopf et al., 2011).

To characterize mesodermal versus endodermal segregation, we used our rabbit/mouse/cynomolgus integrated dataset (Fig. S2) and identified anterior (EPI_ant) and posterior (EPI_post) epiblast within the EPI cells (Fig. 4A). We also identified a minor cluster located between these EPI_ant and EPI_post clusters that we qualified as EPI_int (EPI_intermediate). We then performed a principal component analysis (PCA) on the EPI_post cells (n=239) using 15 markers (5 per tissue: anterior epiblast, endoderm and mesoderm; Fig. 4B). Cells formed two separate clusters, with one cluster having higher expression of mesodermal markers, including TBXT, HAND1, LEF1, PDGFRA and WNT5A, and the other cluster having higher expression of definitive endoderm markers, including CHRD, CER1, GSC, OTX2 and HHEX (Fig. 4B). To localize primordial germ cells (PGCs) in the embryos, we performed immunofluorescence analysis to identify cells positive for both TFAP2C and TBXT (Fig. S3B). These double-positive cells (arrows in Fig. S3B) were observed only in the posterior epiblast at E6.6, not at E6.0. These observations are in line with previous studies reporting segregation between mesoderm and endoderm after E6.0 in rabbit embryos, followed by the emergence of PGCs at E6.6 (Hassoun et al., 2009a,b; Hopf et al., 2011; Kobayashi et al., 2021; Viebahn et al., 2002).

To characterize the transcriptomic changes taking place in the pluripotency continuum of rabbit embryos, we focused subsequent analysis on morula, ICM and EPI cells within the 10x dataset. The resulting ‘PLURI dataset’ included 4315 cells from six developmental stages. We then examined the expression of genes, the expression of which is increased in naive, formative and primed pluripotency (Fig. 5A, Tables S1 and S3). A stronger expression of naive pluripotency genes was observed in morula, ICM and EPI_4.0/5.0 cells, whereas a stronger expression of formative/primed genes was observed in EPI_5.0, EPI_int and EPI_ant cells. Finally, gastrulation-associated genes were mostly upregulated in the EPI_post cells. These results suggest that the epiblast cells of rabbit embryos between E3.5 and E6.6 encompasses the whole pluripotency continuum.

A heatmap was generated with the 20 most overexpressed genes in each cell type. Morula, ICM and EPI_4.0 cells had a higher expression of naive pluripotency genes, including ESSRB, DPPA5, OOEP, BAG3 and FOLR1 (Fig. 5B,C, Table S3). These markers declined in EPI_5.0, concomitantly with the upregulation of a new set of markers, including genes associated with the exit of naive pluripotency ETV4 and TCF7L1. Other genes enriched in EPI_5.0 included cell-cell interaction- and metabolic pathway-related genes, including CLND10, CKB, OAT, ALDH7A1, ANXA2 and S100A6 (Fig. 5B,C). A set of genes including OTX2 was upregulated in both EPI_ant and EPI_post cells (Fig. 5B,C). Finally, several genes were upregulated only in the EPI_post cells, including MIXL1, PITX2 and EOMES, revealing the onset of gastrulation. Differential localization of DPPA5, OOEP, ESRRB, OTX2, FOLR1 and CD57 was confirmed by immunofluorescence (Fig. 5D,E, Fig. S4A); DPPA5, OOEP, ESRRB and FOLR1 were more strongly detected in the EPI at E4.0 compared with E6.0 embryos, whereas OTX2 was more strongly detected in the EPI at E6.0 than in the EPI at E4.0. The primed pluripotency marker CD57 was also strongly detected in the EPI at 6.0. Genes upregulated in epiblast cells (EPI_ant and EPI_post) also included metabolic pathways-related genes (i.e. MT1A, MT2D, FABP5 and FABP7; Fig. 5B). Together, these results suggest alterations in metabolic pathways concomitant with the transition from the naive to primed pluripotency state. This observation was corroborated by GO and KEGG pathway enrichment analysis based on the 500 most differentially expressed genes between early and late epiblast cells (Fig. S4B,C, Table S2). Taken together, these results allowed us to identify gene markers of the naive, formative and primed pluripotency states in rabbit preimplantation embryos.

We investigated the DNA methylation dynamics during rabbit preimplantation development using immunodetection of 5-methylcytosine (5meC) and 5-hydroxy-methylcytosine (5hmeC) in E3.0 to E6.6 embryos (Fig. 6A, Fig. S5A). 5meC fluorescence was low in morulae (E3.0) and early EPI cells (E3.5-E4.0) compared with later stages. There was also higher 5meC fluorescence in EPI cells compared with TE cells in E6.0 and E6.6 embryos. On the other hand, 5hmC signal was higher in early stages (E3.0/E3.5), suggesting active demethylation. The signal then progressively disappeared (E4.0-E5.0). Finally, 5hmC signal could again be detected at E6.0 and E6.6, predominantly in EPI cells. These 5meC and 5hmeC fluorescence patterns are consistent with the increasing expression of DNA methyltransferases (DNMTs) and the decreasing expression of ten-eleven translocation methylcytosine dioxygenases TET1 and TET2, observed during the transition from ICM to EPI (Fig. 6B, Fig. S5C). However, TET3 expression increased between the ICM and late EPI stages, which may explain why 5hmeC is still detected in the EPI at E6.0 and E6.6, despite the increased level of 5meC and increased expression of DNMT genes. Thus, the rabbit embryo shows DNA methylation dynamics similar to those previously described in mice and humans (Zhang et al., 2018; Zhu et al., 2018).

Through the study of X chromosome coating by XIST RNA, it was shown that X inactivation begins at the morula stage in the rabbit (Okamoto et al., 2011). To describe this process in more detail across the pluripotency continuum, we analysed H2AK119ub and H3K27me3 marks, two post-translational modifications of histones associated with X chromosome inactivation in mice (Chaumeil et al., 2011). At the morula stage, H2AK119ub immunostaining appears as small diffuse nuclear spots (E3.0) (Fig. 6C). From the early blastocyst stage onwards, labelling begins to form foci in half of the embryos analysed (n=34). In those, the percentage of cells with a single nuclear focus increased from 15% (E3.5) to 100% (E5.0). Immunostaining of the repressive mark H3K27me3 revealed a similar dynamic (n=12 per stage) (Fig. 6C, Fig. S5B). Consistent with these observations, expression of the H2AK119ub erasers ASXL2 and ASXL3, as well as expression of the H3K27me3 erasers KDM6A and KDM6B was down-regulated during the transition from ICM to EPI_6.0/6.3/6.6 (Fig. 6B, Fig. S5B). From these results, we conclude that X-chromosome inactivation by repressive marks begins as early as E3.5 and is established at E5.0.

The bulk RNAseq dataset was used to investigate the expression of genes associated with oxidative phosphorylation (OXPHOS) and glycolysis. A stronger expression of both nuclear (NC)- and mitochondrial (MT)-encoded genes related to OXPHOS metabolism was observed in TE and ICM cells compared with EPI_6.0/6.3/6.6 (Fig. 7A). Analysis of the single-cell ‘PLURI dataset’ confirmed this observation, showing a gradual downregulation of NDUFV2, UQCRQ, COX8A and ATP5PD gene expression between E3.0 and E6.6 (Fig. 7B). In contrast, LDHA, LDHB and PKM, which are key genes of glycolysis, were expressed at much higher levels in EPI_6.0/6.3/6.6 compared with ICM (Fig. 7A). The single-cell ‘PLURI dataset’ confirmed this finding, showing a gradual increase in the expression of these genes between E3.0 and E6.0 (Fig. 7B). These observations are consistent with a switch from OXPHOS to glycolysis for energy production between early blastocysts and pre-gastrula stage embryos in rabbits, which correlates with observations made in mice (Houghton, 2006). However, other genes involved in glucose uptake and metabolism, including SLC, HK and PDH, were already expressed at early embryo stages (ICM and E3.0 to E4.0), which suggests that the pluripotent cells of the ICM and early EPI are poised for upregulating glucose metabolism-based energy production.

To investigate mitochondrial activity during rabbit preimplantation embryo development, embryos between E3.0 and E6.6 were treated with tetramethylrhodamine ethyl ester (TMRE), which reveals mitochondrial membrane depolarization (ΔΨm). Morula and TE cells showed strong TMRE labelling, indicating high ΔΨm and strong OXPHOS activity (Fig. 7C). TMRE labelling was much lower in the ICM and EPI cells of E3.5 to E6.0 embryos, and became undetectable by E6.6. This result was confirmed by CellROX assay, which labels the reactive oxygen species (ROS) produced by OXPHOS (Fig. S5D). These results are consistent with the expression pattern of NDUF, COX, UQCR and ATP5 families of genes, which are highly expressed in TE cells, and have low expression in EPI_6.0/6.3/6.6 cells. Fluorometric determination of 2-deoxyglucose incorporation was performed to study glycolysis in rabbit embryos (Fig. 7C). Weak 2-deoxyglucose fluorescence was observed in morulae (E3.0), but gradually increased in later stages (from E3.5 onwards), first in TE cells and then in EPI cells. Strong fluorescence was observed in the primitive streak region in E6.6 embryos, correlating with the higher expression of LDHA and PKM observed in transcriptome studies (Fig. 7B).

We sought to further characterize the naive pluripotency state in rabbit embryos and compare with current data in mice and primates. To address this, we analysed the differentially expressed genes between ICM and TE, ICM and PE_6.0, and ICM and EPI_6.0/6.3/6.6 in the bulk RNAseq dataset, which resulted in 1260 upregulated genes in ICM cells compared with the other lineages (P<0.01) (Fig. 8A, Fig. S6A, Table S4). We performed a similar differential analysis between ICM cells (E3.5/E4.0) compared with other cell types (E3.5-E6.6) in the 10x-genomics dataset (average log Fold change>0.25), and identified 161 upregulated genes in ICM cells compared with the other lineages (Fig. 8B, Table S5). We found 95 genes in common between those two sets of upregulated genes (Fig. 8C, Table S5). Among these 95 genes were pluripotency-associated genes BAG3, DPPA5 and SOX15, chromatin regulator-encoding genes KDM4A, KMT2C and SMARCAD1, and metabolic pathways associated genes such as ALDH7A1 (Fig. 8D, Table S5). This observation was corroborated by GO analysis that revealed enrichment for DNA metabolic processes and ‘DNA repair’ (Fig. 8E, Table S5). Notably, both a higher level of KDM4A transcripts at E4.0 and ICM-specific detection of KDM4A protein were observed (Fig. 8F,G). We thus identified key markers of naive pluripotency that might be useful to characterize embryo-derived and induced pluripotent stem cell (ESC and iPSC) lines in the rabbit.

We then further examined the JAK-STAT, NOTCH, WNT, MAPK and TGF-β signalling pathways, all of which are associated with the regulation and maintenance of pluripotency in rodents and primates (Bayerl et al., 2021; Boroviak et al., 2015, 2018; Mohammed et al., 2017; Nakamura et al., 2016). Consistent with naive pluripotency data in mice (Mohammed et al., 2017), the expression of components of the JAK-STAT signalling pathway, including IL6R, IL6ST, IL4R, STAT4 and STAT3, were increased in ICM cells versus TE, PE and late EPI cells (Fig. S6B). Notably, the expression of NOTCH and the NOTCH target gene HES1 were low, whereas the expression of the NOTCH ligand JAG1 and the NOTCH inhibitor NUMB were high in ICM cells compared with other cells (Fig. S6C). Low NOTCH activity in rabbit ICM is consistent with recent data indicating that treatment of human PSCs with the NOTCH pathway inhibitors DBZ and RIN1 enhances naive pluripotency (Bayerl et al., 2021). Overall, these results suggest that the mechanisms controlling naive pluripotency in rabbit preimplantation embryos use signalling pathways that function in mice or primates, or in both.

Finally, we asked which of the aforementioned 95 genes are more highly expressed in the ICM and early EPI cells in both mice and cynomolgus macaques, making use of the single-cell RNAseq data previously generated (Nakamura et al., 2016). As we could not find any orthologs for 14 of the 95 genes, we had to narrow down the analysis to 81 genes (Table S5). In the mouse and cynomolgus macaque datasets, ICM and early EPI cells were identified using landmark markers of naive pluripotency in mice (Fgf4, Esrrb and Klf2) (Fig. S7A) and cynomolgus macaques (FGF4, KLF17 and SOX15) (Fig. S7B). The most overexpressed genes between ICM and early EPI cells compared with all other cell types were identified and represented by heatmaps (Fig. S7A,B, Table S5). Genes more highly expressed in naive cells were represented by volcano plots, highlighting the 81 ortholog genes (Fig. 8H,I). Of these 81 genes, 10 were also upregulated in mice (average log fold-change>0.25) and eight were upregulated in cynomolgus macaques (average log fold-change>0.25). These genes are associated with ‘pluripotency regulation’, ‘chromatin organization’, ‘mRNA splicing’ and ‘DNA metabolism and repair’, as shown in an interaction map of biological processes (Fig. S7C). Finally, a comparison of the naive markers in the three species led to the identification of 11 common genes: GDF3, KDM4C, KDM5B, SMARCAD1, RIF1, MSH6, MCM6, UNG, DPPA5, OOEP and MKRN1 (Fig. 8J). An analysis of OOEP and MKRN1 by immunofluorescence confirmed their detection in the epiblast of E4.0 but not E6.0 embryos (Fig. 8K). Taken together, these results show that rabbit preimplantation embryos share some of the marker genes for naive pluripotency with those of mice, and some with those of cynomolgus macaques.

Our study describes the results of the first thorough investigation of the rabbit preimplantation embryos at the single-cell level. It combines bulk and single-cell RNA-sequencing, protein immunolabelling and fluorometric quantification to characterize the transcriptome, epigenome and metabolome of the pluripotent cells from the morula to early gastrula stage. Four main findings emerge from this study: (1) the three early lineages, EPI, PE and TE, are fully segregated between E5.0 and E6.0; (2) ICM and early (E3.5/E4.0) and late (E6.0/E6.6) EPI cells exhibit the cardinal features associated with naive and primed pluripotency, respectively; E5.0 is a transitional stage in that respect; (3) novel markers of naive pluripotency were identified, including MKRN1 and OOEP; and (4) although the rabbit is evolutionarily closer to mouse than to primates, the transcriptome of rabbit pluripotent cells shares many common features with that of human and non-human primates, including markers of naive pluripotency.

In Eutherian mammals, TE formation is triggered by the asymmetric segregation of keratins and polarization of the outer cells of the morula at the onset of compaction (Lim et al., 2020; Gerri et al., 2020). In our study, the early expression of GATA2 and GATA3 indicates an early onset of the TE program in some morula cells. In contrast, FABP3 expression increases later, at E4.0, in the TE of the expanded blastocyst. Finally, OCT4, which is initially expressed in all cells of the early blastocyst, including the TE, becomes restricted to EPI cells only at E5.0, consistent with a previous study (Canon et al., 2018). These results indicate the progressive differentiation of TE lineage. It is also from E4.0 onwards that the expression of SOX2 and GATA6, which are markers of EPI and PE, respectively, become mutually exclusive and the two cell types separate to form two distinct cell compartments at E5.0, consistent with a previous study (Piliszek et al., 2017). These results indicate that the late segregation of EPI, TE and PE lineages are similar to what has been described in human and non-human primates (Meistermann et al., 2021; Nakamura et al., 2016; Stirparo et al., 2018).

Studying the transition from the morula stage to gastrulation is particularly challenging in many species as the developmental stage at which it occurs corresponds roughly to the time when embryos implant in the uterus. This is particularly true of primates, where access to the newly implanted embryo (between E8 and E13) is virtually impossible. In contrast, rabbit embryos are well suited for investigating this transition because they have delayed embryo implantation until after the onset of gastrulation (Fischer et al., 2012; Puschel et al., 2010). We found that the epiblast of E5.0 rabbit embryos is a transition state, characterized by the downregulation of naive pluripotency markers and the upregulation of formative and primed pluripotency markers. In addition, it is at this stage that we observed the shift to predominantly glycolic metabolism, genome re-methylation and X chromosome inactivation, all of which are considered as key events of the naive-to-primed state transition in rodents and primates (Davidson et al., 2015; Devika et al., 2019). Despite embryo staging prior to cell preparation, some overlap between naive, formative and primed markers can be observed, which may reflect asynchrony of gene expression dynamics. This heterogeneity emphasizes the transitional-stage notion associated with E5.0 rabbit embryos.

Our study identified new marker genes for naive pluripotency in rabbit, a species in which it has not yet been possible to derive truly naive pluripotent stem cell lines (Osteil et al., 2016; Tapponnier et al., 2017). Eleven of these genes are also specifically expressed in the ICM and early EPI in mouse and cynomolgus macaque embryos. GDF3 and SMARCAD1 are well-known pluripotency-associated genes. GDF3 encodes a ligand of TGFβ superfamily that blocks BMP signalling and regulates cell fate in stem cells and embryos (Levine and Brivanlou, 2006). SMARCAD1 encodes a SWI/SNF-like chromatin remodeller with ATPase activity. Its inactivation is detrimental to mouse ESC pluripotency (Sachs et al., 2019; Xiao et al., 2017). Three other genes have also functional characteristics expected of a naive pluripotency regulator. MKRN1 is a new marker of pluripotency. It is a target of OCT4 (Cassar et al., 2015) and encodes an E3 ubiquitin transferase involved in the degradation of p53 and p21, and promotes cell-cycle progression (Lee et al., 2009). DPPA5 encodes an RNA-binding protein that interacts with mRNAs encoding pluripotency and cell-cycle regulators (Tanaka et al., 2006). Moreover, DPPA5 binds NANOG and enhances its function while also preventing its degradation (Qian et al., 2016). OOEP belongs to the same family as DPPA5 (Pierre et al., 2007) and is known to be expressed in the ICM cells and naive PSCs in various species (Messmer et al., 2019; Mohammed et al., 2017; Nakamura et al., 2016; Stirparo et al., 2018). Our study also highlighted the potential role of histone lysine demethylases in naive pluripotency. Three of them are specifically expressed in the ICM and early EPI of rabbit embryos: KDM4A, KDM4C and KDM5B. These enzymes ‘erase’ methylation on histone H3K9 and H3K4, respectively. They have been shown to co-occupy promoters in mouse ESCs (Kumar et al., 2022 preprint), suggesting a functional link to gene expression programmes. Finally, upregulation of MCM6, MSH6, RIF1 and UNG also characterizes naive pluripotency based on our inter-species comparison. These four genes are related to DNA repair and replication, two pathways known to be involved in maintaining the genome integrity and identity of pluripotency cells (Mason et al., 2009).

Although many of the naive and formative/primed pluripotency markers identified in rabbits are common to mice and primates, our study also highlights some noticeable differences. The IL6 receptor-encoding gene IL6R is expressed in the ICM of rabbit and primate embryos, but is replaced by the LIF receptor-encoding gene LIFR in mice (Boroviak et al., 2018). SOX15 is a naive pluripotency marker in rabbits (this study) and primates, but not in mice (Stirparo et al., 2018; Boroviak et al., 2018). TDGF1 is expressed in ICM and early EPI cells but downregulated in late EPI cells in mice (Boroviak et al., 2015; Mohammed et al., 2017), whereas it is expressed only in late EPI in rabbit and cynomolgus macaque (Nakamura et al., 2016). Thus, these genes exhibit an opposite expression pattern in rabbits (this study) and cynomolgus macaques versus mice. These interspecies comparisons clearly suggest a similarity between rabbits and primates in the expression of pluripotency regulators across the pluripotency continuum. Arguably, this makes the rabbit a suitable species for studying the embryo colonization capacity of human naive PSCs and the generation of inter-species chimeras (Aksoy et al., 2021). Nevertheless, the rabbit seems to have some particularities with respect to both primates and rodents. For example, the expression of ESRRB is detected both in ICM and TE in rabbit blastocyst, whereas it is almost absent in mouse TE and human ICM (Blakeley et al., 2015). Moreover, SUSD2, a marker of ICM in primates (Bredenkamp et al., 2019), is highly expressed in rabbit TE. Thus, ESRRB and SUSD2 expression appears to vary considerably across mammals and therefore cannot be considered an unconditional marker of naïve pluripotency. Finally, CCND2 appears to have a different expression pattern than that previously described in the mouse (Wianny et al., 1998). Whereas in the mouse it is activated only during primitive streak formation, it is activated in the epiblast during the naïve-to-primed transition in the rabbit embryo (E5.0). This could reflect a difference in cell cycle regulation of primed pluripotent cells in the rabbit.

All procedures in rabbits were approved by the French ethics committee CELYNE (approval numbers APAFIS #6438 and APAFIS #2180-2015112615371038v2) and COMETHEA (number 45, registered under numbers 12/107 and 15/59). Sexually mature New Zealand white rabbits were injected with follicle-stimulating hormone and gonadotropin-releasing hormone, followed by artificial insemination as previously described (Teixeira et al., 2018). Embryos were flushed from explanted oviducts 65-159 h after insemination.

For bulk RNAseq, morulae were collected at embryonic day 2.7 (E2.7), pooled into groups of 10 and immediately dry-frozen. Blastocysts were collected at E3.5 and E4.0. Mucin coat and zona pellucida (ZP) were removed by protease digestion (Sigma P8811-100MG). Zona pellucida-free embryos were subsequently cultured in TCM199 supplemented with 10% foetal bovine serum (FBS, Sigma M4530) until they regain a normal morphological aspect. Inner cell masses (ICMs) were separated from the trophectoderm (TE) by moderate immune surgery: briefly, blastocysts were incubated in anti-rabbit whole goat serum (Sigma R-5131) at 37°C for 90 min, washed thoroughly and then incubated with guinea pig complement serum (Sigma S-1639) for 5 min. After washing in PBS supplemented with 10% FBS, ICMs were mechanically dissociated from the TE by gentle pipetting with a glass pipette. Samples (ICM and TE) were pooled into groups of 10 and immediately dry-frozen. At E6.0, the mucoprotein layer was removed mechanically using glass microcapillaries. Zona pellucida-free embryos were opened and flattened on a plastic dish in FHM medium to expose the embryoblast with the primitive endoderm on top. The primitive endoderm (visceral part) was first dissociated by careful scratching with a glass needle, and the epiblast was then separated from the trophectoderm with a microscalpel. Samples were dry-frozen. At E6.3 and E6.6, embryos were processed as described for E6.0 embryos. The anterior part of the epiblast was isolated from the posterior epiblast by manual micro-dissection before dry-freezing. For E6.3 embryos, the procedure for isolating the anterior epiblast was validated a posteriori by analysing TBXT expression by RT-qPCR.

For 10x single-cell RNAseq, embryos collected at E3.0, E3.5, E4.0 and E5.0 were treated with protease from Streptomyces griseus (Sigma-Aldrich P8811) at 37°C followed by mechanical dissociation with glass microcapillaries to remove the mucin coat and zona pellucida. For embryos collected at E6.0 and E6.6, mucoprotein layers were removed mechanically with forceps. The embryonic disks were then dissected mechanically with forceps and pooled before single-cell dissociation. For cell singularization, E3.0, E3.5, E4.0 and E5.0 embryos were treated with 0.05-0.1% trypsin for 5-10 min at 37°C. For E6.0 and E6.6 embryos, epiblast cells were mechanically dissociated to obtain small cell clusters (<10 cells), which were then treated with TryPLE for 5 min at 37°C, and singularized by gentle mechanical dissociation. Enzymatic activities were stopped by adding 10% fetal bovine serum (Gibco 11563397). Cell suspensions were run through a 50 µm filter to remove any remaining cell clumps and were loaded on a Chromium controller.

Total RNA was isolated from batches of embryos, ICM or TE (n=20) using PicoPur Arcturus (Excilone, France) with a DNase I (Qiagen, Germany) treatment as recommended by the supplier. For E6.0, E6.3 and E6.6 stages, RNAs were extracted from single embryo samples. Three nanograms of total RNA were used for amplification using the SMART-Seq V4 Ultra Low Input RNA kit (Clontech) according to the manufacturer's recommendations (10 PCR cycles were performed). cDNA quality was assessed on an Agilent Bioanalyzer 2100, using an Agilent High Sensitivity DNA Kit. Libraries were prepared from 0.15 ng cDNA using the Nextera XT Illumina library preparation kit. Libraries were pooled in equimolar proportions and sequenced (Paired-end 50–34 bp) on an Illumina NextSeq500 instrument, using a NextSeq 500 High Output 75 cycles kit. Demultiplexing was performed (bcl2fastq2V2.2.18.12) and adapters were trimmed with Cutadapt1.15, so that only reads longer than 10 bp were kept. The number of reads ranged from 52 to 137 million. Reads were mapped to the rabbit genome (Oryctolagus cuniculus 2.0). 81.8 to 85.7% (depending on samples) of the pair fragments could be aligned; 70.3 to 78% of these fragments passed the mapping filter; 56.4 to 64.7% of them were assigned to a gene.

Cell suspensions were loaded on a Chromium controller (10x Genomics, Pleasanton, CA, USA) to generate single-cell Gel Beads-in-Emulsion (GEMs). Single-cell RNAseq libraries were prepared using Chromium Single cell 3′RNA Gel Bead and Library Kit (P/N 1000075, 1000153, 10x Genomics). GEM-RT was performed in a C1000 Touch Thermal cycler with 96-Deep Well Reaction Module (Bio-Rad; P/N 1851197): 53°C for 45 min, 85°C for 5 min; held at 4°C. After reverse transcription, GEMs were broken and the single-strand cDNA was cleaned up with DynaBeads MyOne Silane Beads (Thermo Fisher Scientific; P/N 37002D). cDNA was amplified using the C1000 Touch Thermal cycler with 96-DeepWell Reaction Module: 98°C for 3 min; 11 cycles of 98°C for 15 s, 63°C for 20 s and 72°C for 1 min; 72°C for 1 min; held at 4°C. Amplified cDNA product was cleaned up with the SPRIselect beads (SPRI P/N B23318). Indexed sequencing libraries were constructed with SPRIselect following these steps: (1) fragmentation end repair and A-tailing and size selection; (2) adapter ligation and cleanup; (3) sample index PCR and size selection. The barcode sequencing libraries were quantified by quantitative PCR (KAPA Biosystems Library Quantification Kit for Illumina platforms P/N KK4824). Sequencing libraries were loaded at 18 pM on an Illumina HiSeq2500 using the following read length: 28 bp Read1, 8 bp I7 Index and 91 bp Read2.

For bulk RNAseq, reads were mapped to the rabbit genome (Oryctolagus cuniculus 2.0) using the splice junction mapper TopHat (version 2.0.14) associated with the short-read aligner Bowtie2 (version 2.1.0). Finally, FeatureCounts (version 1.4.5) was then used to establish a gene count table. Hierarchical clustering was computed using hclust R package with Euclidean distance metric and Ward linkage. Principal component analysis (PCA) was made using FactoMineR (Lê et al., 2008). Data normalization and single-gene level analysis of differentially expression genes were performed using the DESeq2 (Anders and Huber, 2010). Differences were considered significant when adjusted P<0.01 (Benjamini-Hochberg method) and absolute fold change ≥2. Log2 transformed read counts (Rlog) (after normalization by library size) were obtained with DEseq2 package and used for heatmaps.

For single-cell RNAseq, data were analysed using the R software (version 3.6.3) and RStudio Desktop integrated development environment (IDE) for R (version 1.4, Open Source Edition). As the rabbit genome annotation is not complete, some genes were annotated manually using ENSG identity. The Ensembl release 104 (Howe et al., 2021), the National Center for Biotechnology Information (NCBI database) (Schoch et al., 2020) and the g:Profiler online software (Raudvere et al., 2019) were used to convert ENSG annotations into gene symbols. The Seurat package (version 3.1.2) was used to filter, normalize and analyse the data as described previously (Stuart et al., 2019). Briefly, the data were filtered by eliminating all cells with mitochondrial genome expression above 0.7%, and selecting cells in which the total number of expressed genes is between 300 and 2500. A fraction of the cells expressed mitochondrial genes at very low levels. We hypothesized that these were nuclei that had been stripped of their cytoplasm during embryo dissection, single-cell dissociation or single-cell capture in 10x droplets, and were removed from the analyses. In a previous study by Pijuan-Sala and collaborators, similar artefacts were observed and the corresponding events were eliminated (Pijuan-Sala et al., 2019). For the remaining cells, the novelty scores were above 0.80 (Fig. S1D), as expected for good quality cells.

Expression levels were then log-normalized and the 3000 most variable genes were identified. Principal component analysis (PCA) and uniform manifold approximation and projection (UMAP) dimension reductions were applied to the dataset to identify cell clusters. These clusters were identified using graph-based clustering algorithm on the first 20 components and annotated according to expression of lineage-specific genes and embryonic stage. Differential expression analysis was performed and genes with a log fold-change >0.25 or <−0.25 were considered significantly differentially expressed. The same workflow was applied for published data.

Dataset integration was carried out by retrieving the data from Nakamura et al. (2016) on mouse (GSE63266) and cynomolgus monkey embryos (GSE74767) and merging them with our own 10x dataset. The original annotations for the mouse and cynomolgus monkey cells were used. All three datasets were first processed using Seurat (Stuart et al., 2019), as previously described. After normalization, the data were integrated using 4000 genes and 20 dimensions with the Find Integration Anchors function. PCA and UMAP were then used on the merged dataset in order to localize the different embryonic lineages.

Temporal trajectories were created using the following additional packages: monocle (version 2.12.0), cellranger (version 1.1.0) and viridislite (version 0.4.0). Data analysis was performed as described previously (Trapnell et al., 2014). Briefly, differential analysis was performed on cells isolated and annotated using Seurat package, and the top 1000 most significantly differentially expressed genes were used to order the samples in pseudo-time. Stage E3.0 was set as the starting point.

Statistical enrichment analysis of Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) terms was performed using the g:Profiler online software (Raudvere et al., 2019). Cytoscape software (Shannon et al., 2003; Smoot et al., 2011) and enrichment map plug-ins were used to build and visualize networks. The node size is scaled by the number of genes contributing to each term; edge width is proportional to the overlap between each gene set, and each node is coloured by their enrichment score.

Embryos were fixed in 4% paraformaldehyde (PFA) for 20 min at room temperature. After three washes in phosphate-buffered saline (PBS), embryos were permeabilized in (1) PBS-0.5% Tween20 for 15 min (DPPA5 and OOEP detection), (2) PBS-0.5% Triton X100 for 15 min (FOLR1 and CD57 detection) or (3) PBS-0.5% Triton X100 for 30 min (other markers). For anti-5′hydroxymethylcytosine and 5′methylcytosine antibodies, embryos were permeabilized in PBS-0.5% Triton for 15 min, washed in PBS for 20 min, then incubated in 2 M HCl for 30 min. Embryos were subsequently incubated in 100 mM Tris-HCl (pH 9) for 10 min, washed in PBS for 15 min and incubated in blocking solution (PBS supplemented with 5% bovine serum albumin) for 1 h at room temperature. Embryos were then incubated with primary antibodies diluted in blocking solution overnight at 4°C (Table S6). After two washes (2×15 min) in PBS, embryos were incubated in secondary antibodies diluted in blocking solution at a dilution of 1:100 for 1 h at room temperature. Finally, embryos were transferred through several washes of PBS before staining DNA with 4′,6-diamidino-2-phenylindole (DAPI; 0.5 μg/ml) or propidium iodide (PI; 1 μg/ml) for 10 min at room temperature. Embryos were analysed by confocal imaging (DM 6000 CS SP5; Leica). Acquisitions were performed using an oil immersion objective (40×/1.25 0.75, PL APO HCX; Leica).

After gentle removal of the mucin coat with protease from Streptomyces griseus (5 mg/ml), embryos were incubated in pre-warmed pre-equilibrated RDH medium containing 4 µM Hoechst 33342 for DNA staining and either 5 µM CellROX Deep Red reagent, 5 µM Image-iT Red Hypoxia reagent or 50 nM TMRE (all from ThermoFisher Scientific) for 30-45 min at 38°C. For glucose uptake assays, embryos were incubated for 30 min in the glucose uptake mix, according to the manufacturer's instructions (Abcam ab204702). Embryos were then collected and washed twice with the analysis buffer. After staining, embryos were transferred to new drops of RDH before mounting in M2 medium (Sigma) containing 20% Optiprep (Stem Cell Technologies) and ProLong Live antifade reagent for live cell imaging (1/100 dilution, ThermoFisher). Embryos were then imaged on a Leica SP5 confocal microscope with a 25× water immersion objective.

We thank the members of UE1298 SAAJ Science de l'Animal et de l'Aliment, who are responsible for INRAE Jouy-en-Josas rabbit facility, and the high-throughput sequencing facility of the Institute for Integrative Biology of the Cell.

The peer review history is available online at https://journals.biologists.com/dev/lookup/doi/10.1242/dev.200538.reviewer-comments.pdf