Human embryonic genome activation initiates at the one-cell stage

1Fertilizing spermatozoa and metaphase II (mII) oocytes are transcriptionally quiescent. The first transcription in newly formed embryos is known as embryonic genome activation (EGA), but its onset, timing, and profile are poorly understood. In human embryos, EGA is held to have occurred by the eight-cell stage, up to ∼68 h (∼3 days) after fertilization, but this model is likely incomplete. First, it does not accommodate hints that transcription initiates earlier, 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). 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. 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.

2We 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. 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). 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). Outputs from Ingenuity pathway analysis (IPA) corresponded to the initiation of embryogenesis and conclusion of gametogenesis (Figure 1A; Table S1).

3In 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;;; ). 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.

4The 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. Only 3.8% of 3,476 human-sperm-intact RNAs corresponded to upregulated DEGs (FDR < 0.1; Figure S1E), although sperm-associated RNA may not enter mII oocytes. Of the top 200 upregulated genes (FDR < 0.1), 153 were differentially expressed in two-, four-, and eight-cell embryos (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).

5We investigated expression of upregulated transcripts in tripronuclear (3PN) one-cell embryos, which have been adopted for genome editing studies. 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. Moreover, expression of 10 of the 14 targets (71.4%; p < 0.2) was lower in monopronuclear (1PN) one-cell embryos (Figure 2B).

6IPA 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, protective redox-dependent systems, and chromosome segregation, consistent with the elimination of aneuploid two-cell embryos. Upregulated DEGs associated with early developmental arrest included WEE2 (p = 0.006; ) and BTG4 (p = 0.015; ) 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.

7We found no evidence for major upregulation of LINE-1 retrotransposons or transcription factors (TFs) thought to drive cleavage-stage EGA, including OCT4, LEUTX, and DUX4. 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.

8These 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, circumvents protracted genome latency and autonomous regulation, and suggests the trigger for EGA: fertilization.

9Other 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 and mouse pluripotency factor; Esrrb apparently causes less than two-fold upregulation of embryonic-stem-cell-specific genes. The average negative effect of the transcription regulator, Polycomb, is about two-fold, and its associated histone modification, H3K27me3, modulates cellular potency and imprinting. 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, 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, as MYC, MYCN, FOXM1, E2F4, and others are regulated by kinases. 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.

10Transcriptional 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. It is also possible that parentally inherited epigenetic marks (e.g., those associated with acquired obesogenic traits) affect gene expression immediately after fertilization. Maternal factors required to initiate EGA may include polar body markers of oocyte quality and reveal processes in artificial oocytes or embryos necessary for embryonic transcription and totipotency. Although stem-cell-derived human embryos (blastoids) may one day bypass embryonic totipotency, EGA initiation may alternatively be indispensable for normative development. Finally, evaluating genome editing in 3PN embryos should accommodate disrupted gene expression to have clinical utility.

11Human 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.

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

13This study did not generate new unique reagents.

14Human 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.

15Patients 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) 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 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.

16RNA 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.

17For 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.

18Raw read pairs were mapped onto the Human GRCh38 genome and Ensembl 92 gene model using STAR (2.5.0a), 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) 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). For transcript-level analysis, Stringtie was used to generate estimated counts for all detected transcripts.

19Gene- and transcript-level count tables were imported into edgeR 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 ), 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. 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:

20To 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.

21For 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. 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. Transcripts with an FDR <0.1 were considered differentially expressed.

22Pathway 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 were generated using all mapped reads accumulated from oocytes and one-cell embryos using Integrative Genome Viewer (IGV) version 2.4.19.

23For bioinformatic analysis of data from, 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.

24Each 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.