Human embryonic genome activation initiates at the one-cell stage.

In human embryos, the initiation of transcription (embryonic genome activation [EGA]) occurs by the eight-cell stage, but its exact timing and profile are unclear. To address this, we profiled gene expression at depth in human metaphase II oocytes and bipronuclear (2PN) one-cell embryos. High-resolution single-cell RNA sequencing revealed previously inaccessible oocyte-to-embryo gene expression changes. This confirmed transcript depletion following fertilization (maternal RNA degradation) but also uncovered low-magnitude upregulation of hundreds of spliced transcripts. Gene expression analysis predicted embryonic processes including cell-cycle progression and chromosome maintenance as well as transcriptional activators that included cancer-associated gene regulators. Transcription was disrupted in abnormal monopronuclear (1PN) and tripronuclear (3PN) one-cell embryos. These findings indicate that human embryonic transcription initiates at the one-cell stage, sooner than previously thought. The pattern of gene upregulation promises to illuminate processes involved at the onset of human development, with implications for epigenetic inheritance, stem-cell-derived embryos, and cancer.

Introduction

Fertilizing spermatozoa and metaphase II (mII) oocytes are transcriptionally quiescent (Zhou and Dean, 2015). The first transcription in newly formed embryos is known as embryonic genome activation (EGA), but its onset, timing, and profile are poorly understood (Jukam et al., 2017). In human embryos, EGA is held to have occurred by the eight-cell stage, up to ∼68 h (∼3 days) after fertilization (Braude et al., 1988; Leng et al., 2019; Tesarík et al., 1988; Vassena et al., 2011; Xue et al., 2013; Yan et al., 2013), but this model is likely incomplete. First, it does not accommodate hints that transcription initiates earlier (Leng et al., 2019; Wu et al., 2018; Xue et al., 2013; Yan et al., 2013), albeit previous analyses are restricted by poor signal-to-noise ratios, low embryo or donor numbers, and RNA-sequencing (RNA-seq) library preparation protocols that reflect mRNA polyadenylation status (recruitment) (Blower et al., 2013; Oh et al., 2000; Temeles and Schultz, 1997). Second, the model does not explain how the embryo genome is maintained in a transcriptionally silent state during cell proliferation to the eight-cell stage (Alpha Scientists in Reproductive Medicine and ESHRE Special Interest Group of Embryology, 2011). Third, no mechanism has been proposed that explains how maternal factor activity required for early development is regulated at different phases of the first two-to-three cell cycles in the absence of endogenous transcription. Fourth, the model does not address the cue that instigates transcription, which can be provided either in vivo or in vitro or whether it is the cumulative consequence of transcription-independent processes. We, therefore, evaluated the open possibility that gene expression is triggered after fertilization in human one-cell embryos.

Results

We sought to determine the gene expression profile of human one-cell embryos and relate it to development (Figure 1A). Human mII oocytes (n = 12, from seven donors, aged 22.5–31) and bipronuclear (2PN) one-cell embryos (n = 12, from six couples, with no oocyte donor overlap) from various ethnic backgrounds and that appeared healthy (Figure 1B) were subjected to single-cell RNA-seq (scRNA-seq) that avoided poly(A) capture and its attendant potential library bias (Blower et al., 2013; Oh et al., 2000; Temeles and Schultz, 1997). Whole-transcriptome amplification produced indistinguishable yields between oocytes and one-cell embryos (p = 0.595), and in-depth scRNA-seq yielded a mean of 66.3 million reads per cell. Principal component analysis, coupled with t-distributed stochastic neighbor embedding (t-SNE), segregated transcriptomes into discrete oocyte and embryo groups (Figure 1C), and comparison revealed 2,879 differentially expressed genes (DEGs; FDR < 0.1; Figures 1D and S1A). The number of DEGs inversely correlated with fold-change (FC, log2FC; Figure S1B), and there were 1,395 DEGs with an absolute log2FC of ≥0.5 (of which 1,081 survived an FDR cutoff of 0.05). Expression levels of 1,557 genes (54.1% of DEGs; FDR < 0.1) decreased in one-cell embryos, reflecting maternal transcript degradation (Figures 1D and S1C) (Alizadeh et al., 2005). These excluded orthologs of classically down-regulated mouse transcripts, including MOS and GDF9, consistent with relatively slow human maternal transcript degradation lasting several cell cycles (i.e., days) (Leng et al., 2019; Paynton et al., 1988; Petropoulos et al., 2016; Sha et al., 2020; Xue et al., 2013). Outputs from Ingenuity pathway analysis (IPA) (Krämer et al., 2014) corresponded to the initiation of embryogenesis and conclusion of gametogenesis (Figure 1A; Table S1).

Figure 1

Human embryonic transcription initiates at the one-cell stage

(A) Schematic of human one-cell embryo development at the times after fertilization (Capmany et al., 1996). Pb1, first polar body; Pb2, second polar body; PN, pronuclei; PNMB, pronuclear membrane breakdown.

(B) Brightfield images of representative human metaphase II oocytes (mII) and bipronuclear one-cell embryos (emb). Arrowheads indicate pronuclei. Pb1, first polar body; Pb2, second polar body. Scale bar, 50 μm.

(C) t-SNE analysis (without filtering) of single-cell RNA-seq data for human mII oocytes (mII, n = 12) and bipronuclear (2PN) one-cell embryos (emb, n = 12). Circles correspond to oocytes from an African American/Hispanic donor and triangles correspond to embryos from Asian donors, with all other oocytes and embryos being of Caucasian origin and filled triangles and circles corresponding to symbols of (D).

(D) Heatmap showing changes in gene expression levels (FDR < 0.1, log2FC > 0.58) in human mII oocytes (mII) (n = 12) and bipronuclear (2PN) one-cell embryos (emb; n = 12) of (C), indicating donors (top) and the Z score scale (−4 to 4). Each patient is represented by a symbol to indicate the provenance of oocytes and embryos.

(E) Single-cell qPCR of upregulated DEG transcripts in individual human oocytes (mII; n ≥ 3 independent biological replicates) and 2PN one-cell embryos (emb; n ≥ 3 independent biological replicates) (FDR < 0.1, log2FC > 0.5). Different oocytes and embryos were used to those of (D). Corresponding log2FC values from RNA-seq are indicated beneath histograms. Primer pairs flanked exon junctions except in the cases of TIGD5 and NFKB1A (which have a single exon) and PIAS3. Values are ± SEM and normalized against mII oocytes.

(F) Pie chart showing functional classes of upregulated DEGs (FDR < 0.1, cpm > 1.0). Ψ, pseudogene; ERV, endogenous retrovirus; lncRNA, long non-coding RNA; AS, antisense; miRNA, microRNA; si/so, sense-intronic/sense-overlapping; pt, processed transcript.

(G) Pie chart showing processing classes (FDR < 0.34) of upregulated DEGs (p < 0.05; cpm > 1.0). Letters indicate Ingenuity codes.

See also Figure S1 and Table S1.

In addition to down-regulation, expression of 1,322 genes increased in one-cell embryos compared with mII oocytes (FDR < 0.1) (Figure 1D). To test this, we assessed the overlap between upregulated DEGs and 657 transcripts common to four datasets of stably expressed human genes (n = 6,040 ± 2,549 genes; Eisenberg and Levanon, 2013; Lin et al., 2019; Petropoulos et al., 2016). Of these, 542 were expressed in one-cell embryos, most (90.4%) at levels indistinguishable from those in oocytes. Single-cell qPCR of upregulated DEGs in different one-cell embryos mirrored increases revealed by scRNA-seq (log2FC > 0.5; n ≥ 3 different embryos per target; n = 14 targets; Figure 1E). Most (76.3%) upregulated DEGs encoded annotated proteins (Figure 1F) and transcriptome re-assembly de novo predicted the use of canonical transcription start sites in 60.9% (n = 357) of cases and normative splicing in 85.1% (Figure 1G). Sashimi plots demonstrated increases in the levels of mature mRNA transcripts spliced at canonical exon junctions in one-cell embryos, with evidence of embryo-specific alternative splicing (Figure 2A). Exon-flanking qPCR further corroborated mature, spliced mRNA level increases (Figure 1E). These findings reveal the onset of human EGA in one-cell embryos to produce canonically spliced mRNA transcripts.

Figure 2

Human one-cell upregulated gene characteristics and pathways

(A) Raw scRNA-seq density plots (Sashimi plots) along exons and exon junctions. Arcs representing splice junctions connecting exons and display the number of reads split across the junction (junction depth) in mII oocytes (mII) and 2PN one-cell embryos (emb). Genomic coordinates (chrom) and gene annotation tracks are aligned beneath each respective plot. Solid black bars above plots indicate regions of potential alternative splicing.

(B) qPCR for transcripts in individual human monopronuclear (1PN) and tripronuclear (3PN) one-cell embryos (3 ≥ n ≥ 6 biologically independent oocytes or embryos per target). Values for metaphase II oocytes (mII) and bipronuclear one-cell embryos (2PN) from Figure 1E are included for comparison. Values are ± SEM and normalized against mII oocytes (∼1.0). Unpaired t tests indicate p < 0.2.

(C) Ingenuity pathway analysis (IPA) of gene networks upregulated (FDR < 0.1, log2FC > 0) in 2PN one-cell embryos.

(D–F) Upstream transcription regulators inferred by IPA of upregulated gene networks (FDR < 0.1, log2FC > 0) in 2PN one-cell embryos for E2F4 (D), MYC (E), and MYCN (F).

See also Figure S2 and Table S1.

The degree of overlap between upregulated DEGs and previous expression data for cleavage-stage human embryos (Figures S1D and S1E) possibly reflected protocol differences (e.g., poly(A) capture in library preparation), donor ethnicity, and timing of gamete and embryo collection (Leng et al., 2019; Xue et al., 2013). Only 3.8% of 3,476 human-sperm-intact RNAs (Sun et al., 2021) corresponded to upregulated DEGs (FDR < 0.1; Figure S1E), although sperm-associated RNA may not enter mII oocytes (Amanai et al., 2006; Asami et al., 2020; Zhou et al., 2019). Of the top 200 upregulated genes (FDR < 0.1), 153 were differentially expressed in two-, four-, and eight-cell embryos (Leng et al., 2019) (Figure S1F). Most (119/153; 77.8%) exhibited sustained expression throughout two-to-four-cell stages and markedly declined by the eight-cell stage (clusters II and III, Figure S1F).

We investigated expression of upregulated transcripts in tripronuclear (3PN) one-cell embryos, which have been adopted for genome editing studies (Kang et al., 2016; Li et al., 2017; Liang et al., 2015). Of 14 targets upregulated in 2PN embryos that we evaluated, levels of 12 (85.7%) were lower (p < 0.2) in 3PN embryos (Figure 2B). This suggests that EGA initiation was disrupted in most 3PN one-cell embryos, which undergo developmental failure (Kola et al., 1987; Mutia et al., 2019). Moreover, expression of 10 of the 14 targets (71.4%; p < 0.2) was lower in monopronuclear (1PN) one-cell embryos (Figure 2B).

IPA of upregulated DEGs indicated that expressed gene function corresponded to developmental processes in healthy 2PN embryos (Figures 1A and 2C; Table S1). These included ATM activation, which induces G2-/M-phase arrest in response to DNA damage in the mouse (Wang et al., 2013), protective redox-dependent systems (Perreault et al., 1988; Yanagimachi, 1994), and chromosome segregation, consistent with the elimination of aneuploid two-cell embryos (Pauerova et al., 2020; Schatten et al., 1988). Upregulated DEGs associated with early developmental arrest included WEE2 (p = 0.006; Sang et al., 2018) and BTG4 (p = 0.015; Zheng et al., 2020) as well as disease-associated nuclear genes involved in mitochondrial function (e.g., TFB2M, p = 6.71e-05; MFN1, p = 1.34e-05), consistent with important early embryonic roles played by mitochondria (Facucho-Oliveira and St John, 2009; McConnell and Petrie, 2004; Rusecka et al., 2018).

We found no evidence for major upregulation of LINE-1 retrotransposons (Percharde et al., 2018) or transcription factors (TFs) thought to drive cleavage-stage EGA, including OCT4 (Gao et al., 2018), LEUTX (Jouhilahti et al., 2016), and DUX4 (De Iaco et al., 2017; Hendrickson et al., 2017). However, upregulated DEGs included 63 human endogenous retrovirus (hERV) loci (Figures S2A and S2B). Trans-activators of EGA initiation predicted by IPA (FDR < 0.1) included MYC (p = 1.94e-6), MYCN (p = 4.54e-8), RABL6 (p = 0.02), FYN (p = 0.05), and E2F4 (p = 6.48e-12) (Figures 2D–2F and S2C–S2I). Transcripts for a third of these (33/97) were detected by scRNA-seq in mII oocytes and down-regulation of E2F4 correlates with one-cell arrest (Suo et al., 2018).

Discussion

These findings have several implications. First, they suggest that transcription is instated soon after fertilization, during meiotic progression and gamete reprogramming, which should illuminate mechanisms that coordinate chromatin remodeling, the cell cycle, and transcription complex assembly. Second, they provide a readout of epigenetic states that support embryonic transcription and locate genomic addresses of active chromatin. Third, they indicate that decoding gamete chromatin modifications acutely after fertilization might mediate epigenetic inheritance. Fourth, commencement of human gene expression in one-cell embryos accommodates previous hints of early transcription (Wu et al., 2018; Yan et al., 2013; Xue et al., 2013), circumvents protracted genome latency and autonomous regulation, and suggests the trigger for EGA: fertilization.

Other contexts relevant to cellular potency also manifest small magnitude transcript-level changes that may be common in cellular potency transitions. Trans-activation by MYC induces relevant target gene expression less than two-fold (Baluapuri et al., 2020; Nilsson et al., 2005) and mouse pluripotency factor; Esrrb apparently causes less than two-fold upregulation of embryonic-stem-cell-specific genes (Chronis et al., 2017). The average negative effect of the transcription regulator, Polycomb, is about two-fold (Berrozpe et al., 2017), and its associated histone modification, H3K27me3, modulates cellular potency and imprinting (Bernstein et al., 2006; Santini et al., 2021). Several lines of evidence suggest that transcription at the one-cell stage is functional. Upregulated transcripts exhibited a stereotypical pattern that corresponded to early embryonic processes (Perry and Verlhac, 2008; Yanagimachi, 1994; Zhou and Dean, 2015), encoded protein, and utilized canonical exons and transcriptional start sites. Predicted transcription regulator functions mapped to one-cell embryonic processes and the regulation of cancer (Figures 1A, 2C–2F, and S2C–S2I). The switch at fertilization that activates presumptive maternal TFs may involve phospho-relay signaling (Perry and Verlhac, 2008), as MYC, MYCN, FOXM1, E2F4, and others are regulated by kinases (Joshi et al., 2013; Morillo et al., 2012; Sjostrom et al., 2005; Vervoorts et al., 2006). Trajectory analysis (Figure S1F) implied that TFs responsible for early expression are deactivated during the four-to-eight-cell transition, coinciding with a major wave of EGA.

Transcriptional initiation in human one-cell embryos also has clinical implications. Disruption of cancer-associated genes might manifest as both impaired fertility and cancer, and female infertility is indeed associated with ∼15% elevated cancer risk (Murugappan et al., 2019). It is also possible that parentally inherited epigenetic marks (e.g., those associated with acquired obesogenic traits) affect gene expression immediately after fertilization (Huypens et al., 2016; Samata et al., 2020). Maternal factors required to initiate EGA may include polar body markers of oocyte quality (Klatsky et al., 2010; Reich et al., 2011; Swain and Pool, 2008; VerMilyea et al., 2011) and reveal processes in artificial oocytes or embryos necessary for embryonic transcription and totipotency (Yamashiro et al., 2018). Although stem-cell-derived human embryos (blastoids) may one day bypass embryonic totipotency (Liu et al., 2021; Yanagida et al., 2021; Yu et al., 2021), EGA initiation may alternatively be indispensable for normative development. Finally, evaluating genome editing in 3PN embryos (Kang et al., 2016; Liang et al., 2015; Li et al., 2017) should accommodate disrupted gene expression to have clinical utility (Perry, 2000; Yeh et al., 2019).

Limitations of the study

Human 2PN one-cell embryos that may legitimately be used for research are extremely rare, in part reflecting current standard practice among assisted reproduction facilities. This study sourced archival 2PN embryos, placing constraints on sample availability that, inter alia, precluded corroborative analysis orthogonal to high-resolution scRNA-seq, including inhibitor studies (e.g., to block transcription by RNA polymerase II) and confirmatory immunofluorescence (e.g., of predicted maternal TFs). The findings are likely to be evaluated in more tractable (e.g., mouse) models in the foreseeable future, and perhaps eventually in closely related primate (e.g., baboon) models.

STAR★Methods

Key resources table

Resource availability

Lead contact

Further information and requests for resources and reagents should be directed to the lead contact, Tony Perry (perry135@aol.com).

Materials availability

This study did not generate new unique reagents.

EXPERIMENTAL MODEL AND SUBJECT DETAILS

Human oocyte and embryo sample collection

Human oocytes and single monopronuclear (1PN), bipronuclear (2PN) and tripronuclear (3PN) one-cell embryos were supplied anonymously subject to informed consent for use in research by couples who had finished family building or decided for other reasons to discontinue fertility treatment. Consents strictly adhered to guidelines of the Ethics Committee of the American Society for Reproductive Medicine. Embryos were cryopreserved and lysed on site at Ovation Fertility before being anonymized and shipped for analysis.

Method details

Human metaphase II oocytes and one-cell embryos (zygotes)

Patients underwent ovarian stimulation according to guidelines of each clinic, where protocols included agonist luteal phase and antagonist suppression. On the day of retrieval (day 0), mature, metaphase II (mII) oocytes were either cryopreserved by a slow freeze method using propanediol (PROH) (Gook et al., 1993) 3-6 h post-collection, or used to produce embryos by in vitro fertilization (IVF) or ICSI. One-cell embryos used here were morphologically assessed for pronuclear number 19-23 h post-fertilization, cryopreserved 1-5 h later using dimethylsulfoxide (Camus et al., 1989) and stored under liquid nitrogen. Some sibling embryos of morphologically normal bipronuclear (2PN) one-cell embryos (i.e., each containing two pronuclei) gave rise to children. When required, cryopreserved oocytes and one-cell embryos were thawed by rapid warming using a Vit-Warm Kit (FUJI Irvine Scientific, USA) according to the recommended protocol, and viability confirmed. All mII oocytes and one-cell embryos were washed in protein-free multi-purpose handling medium (FUJI Irvine Scientific, USA) and each placed in a 0.2 ml PCR tube containing 0.8 μl 1x single-cell lysis buffer supplemented with RNase inhibitor (Takara Clontech, USA). Oocyte and one-cell embryo (2PN; monopronuclear, 1PN; tripronuclear, 3PN) donor groups did not overlap; mII oocytes and one-cell embryos came from different individuals (Table S3). For scRNA-seq, there were seven mII oocyte donors (six were aged 22.5, 24.5, 25, <30, 27 and 31 years); six were Caucasian and one African American/Hispanic. There were six 2PN one-cell embryo donor couples; for two, male and female ages were respectively 36 and 38, and 40 and 50 (data are unavailable for the other couples) and five of the couples were Caucasian, with one Asian couple.

Single-cell RNA sequencing (scRNA-seq) of oocytes and one-cell embryos

RNA sequencing libraries from 14 mII oocytes and 13 bipronuclear (2PN) one-cell embryos (Table S3) were prepared using the Clontech SMARTer Total RNA-Seq Kit Pico Input (V2) system (Takara Clontech). Briefly, total RNA was liberated by lysis of single oocytes or one-embryos in 0.8 μl lysis buffer supplemented with RNase inhibitor (both from Takara Clontech). The RNA was incubated with SMART Pico N6 primers at 72°C for 3 min and then subjected to first-strand synthesis with SMARTScribe reverse transcription using a Pico v2 SMART adapter (template-switching oligo, TSO). After first strand synthesis, cDNA amplification was performed using SeqAmp DNA polymerase with Illumina barcoded adapters for 5 cycles of 15 sec at 98°C; 15 sec at 55°C; 15 sec at 68°C, followed by final extension for 2 min at 68°C. Ribosomal cDNA was removed using the ZapR v2 and R-Probes v2, after which there was a second round of cDNA amplification with SeqAmp DNA Polymerase for 15 cycles of 15 sec at 98°C, 15 sec at 55°C, and 30 sec at 68°C, to generate the final sequencing libraries. Libraries were analyzed using High-sensitivity D1000 ScreenTape and Agilent TapeStation 4200 (Agilent, USA). All samples were subjected to the same amplification protocol and produced indistinguishable final library yields, with mean±SD respectively of 54.9±25.3 nM and 50.7±9.7 nM for oocytes and one-cell embryos; p=0.595). For next-generation sequencing, libraries were combined at equimolar concentrations before loading onto an Illumina NovaSeq 6000 instrument (Illumina, USA) for paired-end 100 (PE100) sequencing to generate an average of ∼88.2 million raw read pairs per sample.

Ratiometric real-time-PCR (qPCR)

For qPCR, total RNA was liberated by lysis of single oocytes or one-embryos in 0.8 μl lysis buffer supplemented with RNase inhibitor (Takara Clontech). Synthesis of cDNA employed a method modified from the Clontech SMARTer protocol using a total RNA-Seq Kit Pico Input (V2). In brief, total RNA was incubated with SMART Pico Oligos mix v2 at 72°C for 3 min and immediately chilled on ice for 2 min. Samples were then subjected to first-strand synthesis with SMARTSCribe RT (Takara Clontech). First step cDNA amplification was performed using SeqAmp DNA Polymerase, skipping the addition of Illumina barcode adapters and incubating for 1 min at 94°C, followed by 5 cycles of (15 sec at 98°C; 15 sec at 55°C; 15 sec at 68°C) and then for 2 min at 68°C. Second step amplification was performed with SeqAmp DNA Polymerase and PCR2 Primer 2v2, incubating for 1 min at 94°C, followed by 16 cycles of (15 sec at 98°C; 15 sec at 55°C; 30 sec at 68°C) and cDNA stored at -20°C until required. qPCR reactions were performed in a QuantStudio 7 (Thermo Fisher Scientific, UK) or ABI 7500 Real Time PCR System (Applied Biosystems, CA) in reactions (20 μl) containing 1-2 μl template cDNA, forward and reverse primers (100 nM each) and 12.5 μl of Power SYBR (ABI), using standard parameters. Data for each target were obtained from n=3 or 4 biological replicates (i.e., independent single cells) collected on at least two days, and included technical replicates of each biological replicate. Primer sequences are given in Table S2. Primer sets (Sigma-Merck) were non-dimerizing under the conditions employed. Reactions lacking input cDNA were used to verify absence of contamination in cocktail components. Steady state transcript levels were normalized with respect to internal reference, RNA18s5, or in most cases, H3F3A. H3F3A gave mean cycle threshold (Ct) values (±s.e.m.) of 24.51±0.30 for mII oocytes (n=45 replicates over all, n=5 biological replicates) and 24.56±0.33 for 2PN one-cell embryos (n=51 replicates over all, n=8 biological replicates); p=0.829 for mII oocytes vs embryos.

Quantification and statistical analysis

Bioinformatic analysis

Raw read pairs were mapped onto the Human GRCh38 genome and Ensembl 92 gene model using STAR (2.5.0a) (Dobin et al. 2013), employing the following parameters: (--outFilterScoreMinOverLread 0.3 --outFilterMatchNminOverLread 0.3 --outSAMstrandField intronMotif --outSAMtype BAM SortedByCoordinate). Post-alignment reads (BAMs) were further processed by Stringtie (1.3.6) (Pertea et al., 2015) to remodel the transcriptome using default parameters and the Ensembl 92 gene model as base reference. Once the re-assembly was complete, a gene-level count was performed using htseq-count (0.6.1p1) (Anders et al., 2015). For transcript-level analysis, Stringtie was used to generate estimated counts for all detected transcripts.

Gene- and transcript-level count tables were imported into edgeR (Robinson et al., 2010) for downstream differential gene expression analysis. For gene-level analysis, genes expressed at low levels (<1 count per million, >13 samples) were filtered, retaining a final total of 16,625 genes. Due to minimal differences detected in the library yields (see above), samples were subjected to trim-mean of M values (TMM) normalization from the edgeR package (for details, see Robinson and Oshlack, 2010), a common normalization method where a weighted trimmed mean of the log expression ratios is used to normalize sequencing depth. This has been shown to be among the most robust methods for RNA-sequencing differential expression analysis, including single-cell studies where global gene expression differences are large (Evans et al., 2018; Soneson and Robinson, 2018). A subtle normalization factor of 1.00±0.07 (mean±SD) was applied. Normalized gene abundance was listed as read counts per million of mapped reads (CPM) and determined using the formula:

To generate a visual overview of transcriptome profiles, we employed the Rtsne package (https://cran.r-project.org/web/packages/Rtsne/Rtsne.pdf) to perform dimensionality reduction via the principal component analysis (PCA, ndims=50) coupled with t-distributed stochastic neighbor embedding (t-SNE), using the normalized abundance (CPM) of all genes as input with the default parameters. The advantage of using t-SNE over PCA is that it provides superior resolution with which to differentiate samples in fewer dimensions (two projections were used here) and it is commonly used in single-cell RNA sequencing studies. Two mII oocytes were excluded as outliers, and one 2PN one-cell embryo failed to amplify and had to be abandoned (Table S3). The remaining 24 samples (n=12 each oocytes and one-cell embryos) were re-normalized (mean normalization factor ±SD: 1.03±0.02 for mII oocytes and 0.97±0.03 for one-cell embryos) to give a final average of 66.3±9.5 (±SD) million usable read pairs. t-SNE was also re-performed using the same parameters. For differential gene expression, a generalized linear model (GLM) was applied to determine the common, trend and gene-wise dispersions, and likelihood-ratio tests were employed to detect differential gene expression. Genes with a false discovery rate (FDR) of <0.1 were considered differentially expressed.

For transcript-level analysis, transcripts from the 24 samples with low expression were first removed and estimated transcript level counts TMM-normalized (normalization factor = 1.00±0.05) and CPM determined by limma-voom (Law et al., 2014). The total number of transcripts that remained in the analysis was 42,230. We then calculated a gene-variability statistic to adjust for the mean-variance relationship using limma-voom and differential expression was determined using empirical Bayesian t-test (eBayes) from the limma package (Ritchie et al., 2015). Transcripts with an FDR <0.1 were considered differentially expressed.

Pathway and upstream regulator analyes were performed using Qiagen Ingenuity Pathway Analysis (IPA) software using HGNC gene symbols and a raw cut-off of p<0.05. Loci encompassing multiple genes were split in the analysis. Sashimi plots (Katz et al., 2010; Wang et al., 2008) were generated using all mapped reads accumulated from oocytes and one-cell embryos using Integrative Genome Viewer (IGV) version 2.4.19 (Robinson et al., 2011).

For bioinformatic analysis of data from Leng et al. (2019), 65 raw fastq files for individual blastomeres of biparental two-, four- and eight-cell embryo were downloaded via the Gene Expression Omnibus (accession number GSE133854) using the fastq-dump command from the SRA toolkit, with the --split-files argument to split reads 1 & 2. The fastq files were then mapped onto the Human GRCh38 genome, and gene abundance (Ensembl V97) counted using standard STAR 2.5.0a, as described above. The average usable read pair counts (uniquely mapped read pairs to gene) per sample was 24.47 ± 1.2 (s.e.m.) million. Count data were imported into edgeR (3.32.1) and normalized via TMM. Samples were quality checked using the plotMDS() function and samples SRR9645989, SRR9645990, SRR9645991, SRR9645992, SRR9645993, SRR9645994, SRR9645995, SRR9645996 and SRR9645997 were removed, leaving a total of 56 samples. Genes with low expression (<1 CPM in ≥ 46 samples) were removed and counts re-normalized with an average normalization factor of 1.012 ± 0.020. Differential gene expression across all stages was determined by fitting a generalized log linear model followed by likelihood ratio test with a false discovery rate (FDR) cutoff of <0.1; 11,083 (out of 16,810) genes were detected as being differentially expressed. The intersection of the differentially expressed genes in this dataset and the top 200 up-regulated genes (ranked by FDR) from our mII vs one-cell embryo dataset were extracted and we found the majority of the genes (153/200; 76.5%) were differentially expressed in both datasets. A heatmap of the 153 genes was drawn using the pHeatmap package (1.0.12), employing z-scoring to reveal clusters of gene expression trajectories.

Statistical analysis

Each experiment used n≥3 experimental samples, as indicated in the text and figure legends, and are presented as mean ± s.e.m. Statistical differences between pairs of qPCR datasets were analyzed by two-tailed unpaired t-tests. Values of p <0.05 were considered statistically significant unless stated otherwise.

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Biological samples
Human metaphase II (mII) oocytes	Ovation Fertility Austin, Embryology and Andrology Laboratories, Austin, TX 78731, USA	https://www.ovationfertility.com
Human bipronuclear (2PN) embryos	Ovation Fertility Austin, Embryology and Andrology Laboratories, Austin, TX 78731, USA	https://www.ovationfertility.com
Human monopronuclear (1PN) embryos	Ovation Fertility Austin, Embryology and Andrology Laboratories, Austin, TX 78731, USA	https://www.ovationfertility.com
Human tripronuclear (3PN) embryos	Ovation Fertility Austin, Embryology and Andrology Laboratories, Austin, TX 78731, USA	https://www.ovationfertility.com
Critical commercial assays
Clontech SMARTer Total RNA-Seq Kit Pico Input (V2) system	Takara Clontech	Cat # 634412
Deposited data
RNA-seq data	This paper	GEO: GSE157834
RNA-seq data	Leng et al., 2019	GEO: GSE133856
RNA-seq data	Xue et al., 2013	GEO: GSE44183
RNA-seq data	Wu et al., 2018	GEO: GSE101571
RNA-seq data	Sun et al., 2021	GEO: GSE137490
Experimental models: Organisms/strains
Human (Homo sapiens)		N/A
Oligonucleotides
Primers for qPCR: see Table S2	This paper	N/A
Software and algorithms
Human GRCh38 genome and Ensembl 92 gene model using STAR (2.5.0a)	Dobin et al., 2013	N/A
Stringtie (90 1.3.6)	Pertea et al., 2015	N/A
htseq-count (0.6.1p1)	Anders et al., 2015	N/A
edgeR	Robinson et al., 2010	N/A
trim-mean of M values (TMM) normalization from the edgeR package	Robinson and Oshlack, 2010	N/A
Rtsne package		https://cran.r-project.org/web/packages/Rtsne/Rtsne.pdf
limma package	Ritchie et al., 2015	N/A
Qiagen Ingenuity Pathway Analysis (IPA) software		https://digitalinsights.qiagen.com/products/qiagen-ipa/latest-improvements/current-line/
Sashimi plots	Katz et al., 2010; Wang et al., 2008	N/A
Integrative Genome Viewer (IGV) version 2.4.19	Robinson et al., 2011	N/A