This study was carried out in strict accordance with the recommendations in the Guide for the Care and Use of Laboratory Animals of the National Livestock Research Institute of Korea. The protocol was approved by the Committee on the Ethics of Animal Experiments of the Korea Research Institute of Bioscience and Biotechnology.

Bovine oocytes were collected from ovaries supplied by a local slaughterhouse and matured in the paraffin oil covered in vitro maturation medium for 20 h at 38.5°C with 5% CO2. The medium for oocytes maturation was prepared by combining TCM-199 (Invitrogen) supplemented with 10% (v/v) fetal bovine serum (FBS; Invitrogen), 10 μg/ml FSH-P (Folltropin-V, Vetrepharm), 0.6 mM cysteine, 0.2 mM sodium pyruvate, and 1 μg/ml estradiol-17β together. To generate IVF embryos, the matured oocytes were fertilized by incubating with 2 × 106 sperms/ml in fertilization medium at 38.5°C in 5% CO2 for 20–22 h (Park et al., 2007). After the insemination, cumulus cells were removed by gentle pipetting, and the fertilized eggs were further cultured in CR1aa supplemented with 3 mg/ml BSA (fatty acid free). After 3 days, cleaved embryos were cultured in CR1aa (with 10% FBS) for 4 days at 38.5°C in 5% CO2 (Koo et al., 2002).

For the generation of bovine scNT blastocysts, bovine mature oocytes were manipulated as described elsewhere (Koo et al., 2002). Oocyte manipulations including enucleation were performed by using a micromanipulator equipped with an inverted microscope (Leitz, Ernst Leitz Wetzlar GmbH). The medium containing TL-Hepes with 7.5 μg/ml cytochalasin B was used for manipulation. The first polar body and a part of the cytoplasm were removed together by a micropipette, and single cells were individually transferred to the perivitelline space of the recipient oocytes. The donor cell containing oocytes were equilibrated for 10–20 s in 50 μl of cell fusion medium and transferred into a fusion medium containing 0.01% BSA, 0.1 mM CaCl2, 0.3 M mannitol, 0.5 mM Hepes, and 0.1 mM MgCl2. The donor cells were fused into the oocytes by a single pulse of direct current of 1.6 kV/cm for 20 μs each by an Electro Cell Manipulator 2001 (BTX). After 1 h, the oocytes with no visible somatic donor cells in the perivitelline space were selected, and they were activated with 5 μM Ionomycin for 5 min, followed by treatment with 2.5 mM 6-dimethyl-aminopurine (DMAP, Sigma) in CR1aa supplemented with 10% FBS for 3.5 h at 38.5°C in 5% CO2 in air. The activated reconstructed oocytes were cultured in the same conditions as IVF embryos for 7 days until they formed blastocysts. As for donor cells, ear skin fibroblasts were obtained from an adult female or male cow and passaged 3 times in the standard culture condition as described before (Koo et al., 2002).

For generation of sham NT blastocysts, the zygotes presenting both parental pronuclei were manipulated at ~15 post-IVF h. For the precise imitation of the physical damages by enucleation, zona pellucida was partially ripped, and the polar body and a part of the underlying ooplasm were carefully removed by aspiration using a micropipette without disturbing pronuclei. The oocytes were activated, after 2 h of incubation, using 5 μM ionomycin (Sigma) for 5 min, followed by treatment of 2.5 mM DMAP in CR1aa culture media supplemented with 0.3% BSA for 3.5 h. Blastocysts were generated 7 days post-IVF or -NT. The quality of each blastocyst was assessed by Hoechst staining, and only high-quality ones at mid blastocyst stage possessing 60–80 blastomeres were chosen for transcriptomic analysis.

For chemical activation of mouse zygotes which were used to examine the effect of chemical-mediated secondary activation on their ability of normal development, fertilized mouse oocytes were collected from superovulated C57BL/6 females as described previously (Hogan, 1994). Briefly, female C57BL/6 mice at 5 weeks of age were injected with 5 IU of pregnant mare serum gonadotrophin, followed by 5 IU of human chorionic gonadotropins 48 h apart, and mated with male mice. Successful mating was determined the following morning by detection of a vaginal plug. Mouse zygotes were collected from mouse oviduct and transferred to M2 medium (Sigma) containing 0.1% (w/v) hyaluronidase to remove cumulus cells and cultured in M16 medium (Sigma) at 37°C, 5% CO₂ in air (Yeo et al., 2005). The mouse zygotes were subjected to the same protocol of chemical activation as the bovine zygotes (see above) and cultured to the 2-cell stage. The cleaved 2-cell embryos were transferred to a pseudo-pregnant surrogate mother. Fetal development was observed at 13.5 days post coitem (dpc).

Each single blastocyst generated by IVF, scNT, or sham NT was lysed, and mRNAs were directly isolated from the lysates using Dynabeads mRNA DIRECT Kit (Invitrogen) according to the manufacturer's recommendation. The ultra-pure mRNAs were reverse transcribed using a custom pico-profiling adapter as previously described (Min et al., 2015). First-strand cDNAs were synthesized by incubating mRNAs with 100 pmol of pico-profiling adapters and 1 μl SuperScript III enzyme (Invitrogen) in 15 μl reaction volume using a following sequential program: 18°C for 10 min, 25°C for 10 min, 37°C for 30 min, 42°C for 10 min, and 70°C for 20 min. Then, 1 μl T4 DNA polymerase (NEB) were added and further incubated at 37°C for 1 h to tag both ends of cDNAs. For amplification of the cDNAs, 15–20 cycles of PCR was performed 94°C for 2 min, 70°C for 5 min using an anchor primer. The amplified cDNA libraries were digested by MlyI enzyme (NEB) overnight and then purified using Agencourt AMPure XP beads (Beckman Coulter) before used for NGS library construction.

A total of 35 RNA-seq libraries (12 IVF, 12 scNT, and 11 shNT blastocysts) were constructed using TruSeq DNA Sample Preparation kit (Illumina) (Min et al., 2015, 2017). Two hundred nano grams of pico-profiled amplicons were applied to the End-repair reactions by incubating with End Repair Mix (Illumina) at 30°C for 2 h. The reactions were purified with AMPure XP beads (Beckman). Then, A-Tailing Mix (Illumina) was added to the samples, and the mixtures were incubated at 37°C for 2 h. Next, indexed adapter was added into each sample with Ligation Mix (Illumina) and the mixtures were incubated overnight at 16°C. Ligates were purified using AMPure XP bead (Bechman). For size selection, adapter ligated DNA samples were loaded onto Pippin Prep (Sage Science), and DNA fragments between 300 and 500 bp were isolated. Finally, the isolated DNAs were enriched by 15–18 cycles of PCR reaction. The enriched DNA fragments were purified and sequenced (2 × 100) using HiSeq2500. All raw sequencing data are available at Gene Expression Omnibus (GSE95311) and our previous publication (PMID: 28443134).

Raw read sequences (2 × 100 bp) from the HiSeq2000 system were preprocessed to remove NGS artifacts such as adapter sequences and low quality bases by “trimgalore,” and the trimmed reads were aligned on the reference genome (UMD3.1) using TopHat (Kim et al., 2013) with the seed length (L) of 10 to increase specificity. We followed the cufflinks pipeline (Trapnell et al., 2012) to estimate the gene expression levels and identify the DEGs. All cufflinks suits were operated with their default options unless noted otherwise. First, transcriptomes of individual samples were assembled using “cufflinks” with –G option to obtain the expression levels of only known genes due to the incomplete genome assembly of bovine. The transcriptomes were merged into one by “cuffmerge,” and differential expression (DE) analysis was performed by “cuffdiff.” Pearson correlation coefficients between samples and groups were calculated by “cor,” an R function. Gene ontology and KEGG pathway analysis was performed using DAVID Bioinformatics Resources 6.8 (Huang da et al., 2009). Weighted gene co-expression network analysis (WGCNA) was executed as described elsewhere using the power of β = 5 (Langfelder and Horvath, 2008; Xue et al., 2013). All plots in this study were produced by R scripts, MS EXCEL, and Origin 8.

We compared the transcriptomes of three different kinds of bovine blastocysts derived by IVF, scNT, and shNT procedures. The comparisons enabled us to separately examine the effects of various NT-involved factors such as the manipulation, partial removal of oocyte/zygote cytoplasm, and donor genome reprogramming (Figure 1A). We analyzed previous Pip-seq data (Min et al., 2017) obtained from individual blastocysts after PCR-mediated amplification of embryonic transcripts by pico-profiling (Min et al., 2015, 2016). Each blastocyst group comprised 12 blastocysts, with six males and six females (with the exception of five female shNTs). We examined the gene expression profiles, particularly of shNT blastocysts, to determine how the transcriptome differed from IVFs and scNTs. The shNTs were prepared as depicted in Figure 1B. Briefly, the zygote with small pronuclei (~15 h after sperm exposure) was put through the NT procedure, in which a part of its cytoplasm was removed carefully to avoid losing genomic material. Therefore, the zygote and its descendants kept the genomic material intact in the nucleus. Meanwhile, the sham-manipulated zygotes were treated with the chemicals for artificial activation (DMAP and Ionomycin), which aimed to reiterate NT procedure equally in the sham zygotes. The in vitro developmental rates of the IVF, scNT, and shNT groups are summarized in the Figure S1A. Only high-quality blastocysts at mid blastocyst stage possessing 60–80 blastomeres were chosen for transcriptomic analysis (Figure S1B). Meanwhile, we examined whether the repeated activation (fertilization and artificial activation) could be intolerable and would be catastrophic to the bovine zygotes as it is thought. When we equally treated mouse zygotes with the same activating chemicals, we observed that mouse zygotes could normally develop as the control embryos, indicating that the repeated activation is not disastrous to embryo development (Figure S1C). Manipulating the zygotes instead of the oocytes and retaining an intact zygote genome were the only differences in the sham NT protocol compared with the standard NT protocol (Park et al., 2007; Kwon et al., 2015).

The number of genes expressed in IVFs, scNTs, and shNTs was 14,095, 13,840, and 13,471, respectively. Although slightly lower in the shNTs, the number of genes was not significantly different (p = 0.067; one-way ANOVA) among the three groups. From principal component analysis (PCA), we found that the shNTs separately grouped from the IVFs and the scNTs, while the latter two associated together (Figure 1C). Assessment of correlations by unsupervised clustering revealed that shNTs were weakly correlated with IVFs and scNTs (Figure 1D). Pearson correlation between the groups showed that the highest correlation came from the IVF-scNT match (r² = 0.895) and the next highest from IVF-shNT (0.812) and shNT-scNT matches (0.703; Figure 1E). These data indicate that the transcriptomes of shNTs are very distinct from those of IVFs and scNTs.

We identified 1,873 DEGs (fold-change > 2.0 and p-value < 0.05) between shNTs and IVFs (Figure S2). These were assumed to reflect the physical (from zygote manipulation) and biochemical (from a partial loss of the cytoplasm) damage in the shNTs (see Figure 1A). Among them, a half (937 DEGs) were underrepresented in shNTs, whereas the other half (936 DEGs) were overrepresented. This was greater than the number of DEGs obtained from the IVF-scNT comparison (256 scNT-low and 183 scNT-high DEGs) and from the shNT-scNT comparison (889 shNT-high and 819 shNT-low DEGs). These results suggest that the IVFs-shNTs comparison has the largest transcriptomic difference. The DEG information is summarized in Supplementary File 1.

The DEGs detected between IVFs and shNTs could be either manipulation-dependent or -independent. This could be discriminated by including an additional set of DEGs, the IVF-scNT DEGs, in the simultaneous comparison. Given that shNTs have the same nuclear material as IVFs, those DEGs that were commonly detected in the IVF-scNT and IVF-shNT DEG sets are likely to be manipulation-dependent. Here, we call these “manipulation-associated DEGs (MADs).” In contrast, those DEGs present in the IVF-scNT DEG set but absent from the IVF-shNT set may be donor genome-dependent. These present as a signature of inefficient reprogramming, which we call the “donor genome-associated DEGs (DADs).”

As shown in the Venn diagram in Figure 2A, 561 DEGs (p < 0.05) were identified as MADs. Heat-map analysis shows a collection of 261 overrepresented and 300 underrepresented MADs, in which the levels differed from those of the IVFs (Figure 2B). GO enrichment analysis yielded several terms in the “biological process” category, which we summarize in Table 1. Some representative terms such as “cell-cell adhesion” (p = 1.98 × 10⁻⁴), “translation” (p = 2.55 × 10⁻³), and “transcription” (p = 0.031) are shown as expression heat-maps with gene symbols in Figure 2C. MADs representing “cell-cell adhesion” and “translation” functions were generally underrepresented, whereas those representing “transcription” functions were largely overrepresented in the scNTs and shNTs.

The DADs outnumbered the MADs (Figure 2D). Seven hundred and sixty-four DEGs were identified as DADs, with 402 scNT-high and 362 scNT-low DADs (Figure 2E). GO enrichment analysis yielded various “biological process” terms, which are summarized in Table 1. Most (11/16) of the DADs belonging to “cell differentiation” functions were underrepresented in scNTs compared to that in IVFs and shNTs (p = 0.027; Figure 2F). DADs associated with “embryo implantation” functions also showed decreased expression in scNTs (p = 0.014). Further information on the MADs and DADs is presented in Supplementary File 2. Additionally, we examined the enrichment of KEGG pathways associated with MADs or DADs and found enrichment of pathways involved in protein synthesis such as “ribosome” and “ribosome biogenesis in eukaryotes” for MADs, while DADs showed significant association with various disease related pathways including Huntington's disease and Alzheimer's disease. The full list of significantly enriched pathways are presented in Supplementary File 2.

A relatively large number of genes (n = 2,921, p < 0.05) were found as shNT-specific DEGs (Figure 3A), which was too large in number to run a GO engine. The output of GO analysis with the 2,921 genes that account for nearly one-third of expressed genes and are presumably implicated in that much of biological processes would be unambiguously enormous and confusing, thus definitely irrelevant to find significant and valuable terms from the result (see Supplementary File 2 for GO result using 2,921 genes). Therefore, we raised the stringency of DEG selection (p < 0.05 plus a fold-change > 2.0). Under these new criteria, 889 shNT-high and 854 shNT-low DEGs were identified (Figure 3B). They are listed in Supplementary File 2. These shNT-specific DEGs reflect the effect of zygote manipulation (Figure 1A) and therefore we call them “zygote manipulation-associated DEGs (zMADs).” GO enrichment analysis yielded a variety of terms, some of which are shown in Figure 3C. Categories most highly ranked were those implicated in protein processing such as “translation” (p = 1.40 × 10⁻⁶), “protein folding” (p = 1.58 × 10⁻⁵), and “protein stabilization” (p = 2.98 × 10⁻³). DEGs implicated in “rRNA processing” (p = 7.24 × 10⁻⁵) and “double-strand break repair” (p = 2.12 × 10⁻³) functions were characterized by higher expression levels in the shNTs. “Apoptosis pathway” (n = 35; p = 7.57 × 10⁻⁴) and “I-κB kinase- and NF-κB signaling” (n = 22, p = 7.84 × 10⁻³) functions were also detected with statistical significance.

Because of the opposite expression patterns among genes in the GO terms (Figure 3C), we re-generated GO terms separately using the shNT-high and the shNT-low zMADs (Table 2). Among the GO terms generated from the shNT-low zMADs, several were associated with protein synthesis and protein stability, such as the “translation,” “protein folding,” “protein stabilization,” and “protein ubiquitination” terms (Figure 4A). Downregulation of genes associated with these terms may disturb protein homeostasis and related cellular processes, suggesting an innate weakness of the shNTs in terms of protein manufacturing and processing. Other GO terms frequently shown from the analysis of shNT-low zMADs were related to mitochondrial functions. These included the “regulation of mitophagy,” “mitochondrial translational initiation/elongation,” “positive regulation of protein targeting to mitochondria,” and “protein import into mitochondrial matrix” (refer to Table 2 for those not shown in Figure 4A). Early embryos rely on the mitochondrial pool inherited from the oocyte and it is not until the blastocyst stage that mitochondria commence replication (Dumollard et al., 2007). Therefore, the reduced pool of mitochondria left in the recipient cytoplasm can interfere with proper development of the shNT embryo. Genes in the term “regulation of mitophagy,” a process of autophagy-mediated selective degradation of defective mitochondria (Lemasters, 2005), were expressed in shNTs at only 34% of the level in IVFs (Figure 4B). From this, we infer that shNTs may keep saving normal-looking, but functionally defective, mitochondria in the cytoplasm. A similar level of downregulation was observed with genes in the “protein targeting to mitochondria” term (Figure 4C). Given that the same sets of genes are expressed in scNTs at ~ 85% of the level of expression in IVFs, we assume that the zygote manipulation is more detrimental to mitochondrial homeostasis in early development.

GO enrichment analysis of shNT-high zMADs generated several protein transport-related terms, including “endoplasmic reticulum (ER) to Golgi vesicle-mediated transport” and “intracellular protein transport” (Figure 4D, Table 2). shNTs had a 2.7-fold higher expression level for genes that belonged to these categories than the IVFs, whereas the scNTs showed no difference (Figure 4E). The shNT-specific increase in expression of genes with the GO keywords “ER,” “Golgi,” and “vesicle transport” implicates the ER stress response (Hetz et al., 2015). In reality, the GO result from shNT-high zMADs showed the term “response to ER stress” (p = 0.02) and the corresponding genes had 3.5-fold higher expression compared to those in IVFs (Figure 4F). In line with this, ATF6, which encodes an ER transmembrane glycoprotein that acts as a transcription activator and initiates the unfolded protein response (UPR) during ER stress (Wang et al., 2000), had ~4-fold increased expression in shNTs compared to that in IVFs (p = 2.23 × 10⁻⁶). However, the expression level of ATF6B, which functions to reduce expression of ER stress proteins, was not different (p = 0.276) between shNTs and IVFs (Thuerauf et al., 2004). KEGG pathway analysis also showed that zMADs were associated with “protein process in ER” related pathways such as “protein processing in endoplasmic reticulum” and “ribosome” in shNT (Supplementary File 2). Therefore, these results suggest that shNTs suffer from ER stress.

Of particular interest was whether the expression of those genes important for early development were affected in shNTs. We examined the expression of epi-driver genes, stemness genes, and genes associated with trophectoderm (TE) development (Min et al., 2015; Park et al., 2017). Figure 5A shows the heat-map of expression levels of epi-driver genes in bovine blastocysts. The median expression levels in the three blastocyst groups were very similar (~17) in fragments per kilobase of exon per million mapped fragments (FPKM), showing no significant difference between the groups (p > 0.822, paired sample t-test; Figure 5B). When we looked closely at the heat-map, the upper part containing genes overrepresented in shNTs was occupied by genes implicated in heterochromatin establishment and maintenance (KDM5A, KDM5B, SETDB1, SUV420H1, HDAC8, DNMT3A, and KDM2A). This suggests that shNTs treat heterochromatin favorably. In addition, several protein arginine methyltransferase genes (PRMTs; PRMT9/10, PRMT5, and PRMT7) and sirtuin genes (SIRT4 and SIRT7) were found in that section of the heat-map. Interestingly, of the seven PRMTs that were expressed in cow blastocysts, six were overrepresented in shNTs (Figure S3), although the implication and consequence of this finding to the genome are unknown. In fact, those epi-driver proteins have diverse functions in chromatin modification as writers, erasers, and keepers of the epigenome. Different combinations of the type and amount of epi-drivers can affect a variety of synergic or opposing changes in the configuration and structure of chromatin, reshaping chromatin by either packing or relaxing it. Hence, we could not conjecture to which direction the altered expression of epi-driver genes would lead; we simply know that the altered genes expression of epi-driver indicates altered chromatin compaction or relaxation in shNTs.

Some TE genes, which are important for implantation and placental development and are frequently abnormal in NT embryos (Kang et al., 2002; Koo et al., 2002), were overrepresented (DKKL1 and CDH1; p < 0.01) or underrepresented (ESRRB, GATA2, PITPNC1, and CD9) in shNTs compared to that in IVFs (Figure 5C). In agreement, the GO enrichment analysis using shNT-low zMADs yielded the “embryonic placenta development” term (GATA2, HIF1A, EPAS1, CITED2; Table 2). In the IVF-scNT comparison, fewer genes (GATA2, PITPNC1, EPCAM, and CDX2) were detected even at a lower stringency (p < 0.05). Thus, GATA2 and PITPNC1 were shared as DEGs in common between shNTs and scNTs. Of the four scNT DEGs, the latter three coincided with those previously identified as TE-related DEGs in scNT blastocysts (Min et al., 2015).

The majority of stemness genes were underrepresented in the shNTs compared to that in the IVFs (Figure 5D). The median expression level of stemness genes was significantly different (p = 0.0020) between IVFs and shNTs, but not in other comparisons (Figure 5B). For example, OCT4 (or POU5F1) and KLF5 expression levels were significantly reduced in shNTs but not in scNTs (Figure 5E). Figure S4 shows the coverage plots for representative stemness genes. In the same context, among the stemness genes that were differently expressed in the shNTs (p < 0.05) compared to that in the IVFs, about 75% (24/32) were underrepresented, while only 43% (16/37) of epi-driver genes and 34% (60/178) of chromosome 1 (chr1) genes were underrepresented (Figure 5F). TE genes showed a pattern similar to that shown by the stemness genes, with 78% (7/9) of TE DEGs (p < 0.05) underrepresented in shNTs. Because of the small number of DEGs, a comparison of gene expression between scNTs and shNTs was not possible. Meanwhile, there were some stemness genes, such as BMPR1A (p = 2.43 × 10⁻⁸) and DPPA3 (p = 1.61 × 10⁻⁴), that were significantly overrepresented in shNTs (Figure S4).

We additionally performed WGCNA to check the blastocyst-type-specific co-expression networks (modules) and summarized in Figure S5 and Supplementary File 3.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.