Single-cell RNA-seq and single-cell bisulfite-sequencing reveal insights into yak preimplantation embryogenesis

Abstract

Extensive epigenetic reprogramming occurs during preimplantation embryonic development. However, the impact of DNA methylation in plateau yak preimplantation embryos and how epigenetic reprogramming contributes to transcriptional regulatory networks are unclear. In this study, we quantified gene expression and DNA methylation in oocytes and a series of yak embryos at different developmental stages and at single-cell resolution using single-cell bisulfite-sequencing and RNA-seq. We characterized embryonic genome activation and maternal transcript degradation and mapped epigenetic reprogramming events critical for embryonic development. Through cross-species transcriptome analysis, we identified 31 conserved maternal hub genes and 39 conserved zygotic hub genes, including SIN3A, PRC1, HDAC1/2, and HSPD1. Notably, by combining single-cell DNA methylation and transcriptome analysis, we identified 43 candidate methylation driver genes, such as AURKA, NUSAP1, CENPF, and PLK1, that may be associated with embryonic development. Finally, using functional approaches, we further determined that the epigenetic modifications associated with the histone deacetylases HDAC1/2 are essential for embryonic development and that the deubiquitinating enzyme USP7 may affect embryonic development by regulating DNA methylation. Our data represent an extensive resource on the transcriptional dynamics of yak embryonic development and DNA methylation remodeling, and provide new insights into strategies for the conservation of germplasm resources, as well as a better understanding of mammalian early embryonic development that can be applied to investigate the causes of early developmental disorders.

Keywords: embryo development; embryonic genome activation; single-cell RNA-seq; whole-genome bisulfite sequencing; yak.

Figure 1

Transcriptional profiles of yak in vitro MII oocytes and preimplantation embryos.A, representative images of transcriptional profiling. The scale bar represents 100 μm. B, total number of genes in each sample. C, venn diagram of differentially expressed genes at five stages. D, heatmap of correlation: The values in the grid represent correlation coefficients. E, heatmap after hierarchical clustering of the samples. The intensity of the red color reflects the degree of increased expression. F, principal component analysis (PCA). The same symbol indicates repetition of the sample.

Figure 2

Gene expression patterns during preimplantation developmental stages.A–C, WGCNA of each sample of yak. A, cluster dendrogram showing 10 modules corresponding to 10 colors. Height (Y-axis) indicates the degree of correlation. B, heatmap of correlations (and corresponding p-values) between coexpressed modules and various periods. C, functional terms of stage-specific genes. D, data have been partitioned into 6 clusters using k-means. E, heatmap of embryonically expressed genes. F, heatmap of maternally expressed genes. WGCNA, weighted gene coexpression network analysis.

Figure 3

Transcriptional patterns during the period of embryonic genome activation.A, number of DEGs between any two consecutive stages. B and C, gene ontology classification of DEGs. The “Percentage of genes” in the GO graph refers to the proportion of upregulated differentially expressed genes (DEGs) annotated to GO entries among the all upregulated DEGs (in red) or the proportion of downregulated DEGs annotated to GO entries among all downregulated DEGs (in green). D, significant functional pathways revealed from upregulated or downregulated genes between 2-cell and MII oocytes stages. E, enrichment plots from GSEA. F and G, significant functional pathways revealed from upregulated (F) or downregulated (G) genes between morula and 8-cell stages. GO, gene ontology; GSEA, gene set enrichment analysis.

Figure 4

Heatmaps of overlapping genes in stage-specific modules across multiple species.A–C, from left to right, the joint analysis heatmaps of yak and cattle, yak and human, and yak and mouse. Each cell includes the number of intersecting genes and the p-value (−log10). D, venn diagram of conserved maternal genes across species. E, maternal hub genes. F, venn diagram of conserved zygotic genes across species. G, zygotic hub genes.

Figure 5

Genome-wide methylation dynamics during preimplantation embryonic development in the yak.A, correlation analysis of methylation levels for all samples. B, whole genome DNA methylation profiles for all samples. C, mapping of methylation level distribution in the whole gene 2-kb window for all samples. D, distribution of methylation levels of CGI in all samples. E–G (E–G) are the methylation levels and gene densities of CG, CHG, and CHH, respectively. Circos plot of chromosome methylation levels: from outer to inner circles are, stain map scale, a (for Yak2cell), b (for Yak4cell), c (for Yak8cell), d (for Yakmorula). H, IGV visualization of methylation kinetics in the X chromosome genomic region for all samples. I, methylation levels and gene density of the X chromosome. Circos plot of chromosome methylation levels: from outer to inner circles are, stain map scale, a (for Yak2cell), b (for Yak4cell), c (for Yak8cell), and d (for Yakmorula). Expression levels of all genes are in the black line and XIST genes are in the red line. CGI is indicated by the blue line. J, IGV visualization of DNA methylation levels (BS2cell01-BSmorula3) and transcriptome expression levels (YMII01-YMorula03) of the XIST gene in all samples. CGI, CpG islands; IGV, integrative genomics viewer.

Figure 6

Clustering analysis of functional region methylation annotations and differential loci.A, distribution of methylation levels in different genomic functional regions. B, distribution of methylation levels of CG, CHG, and CHH in genomic regions in different cytosine contexts. C, clustering diagram of methylation levels of CpG loci in all samples. Each row and column represents a CpG site and cell sample respectively, and the right scale values are methylation levels (0–100).

Figure 7

Differential methylation analysis at the 8-cell embryos stage and the morula stage.A, length of DMRs: hyper is high methylation and hypo is low methylation. B, functional annotation of DMRs. C, position of DMRs in the overall genome. D, GO functional enrichment analysis. Scatter plots are used to show the functionally significant pathways recognized from hyper (E) or hypo (F) DMRs between the morula and 8-cell stages. DMRs, differentially methylated region; GO, gene ontology.

Figure 8

Correlation of dynamic changes in transcriptome and methylationome.A, scatter plot showing the correlation between gene expression levels and corresponding promoter methylation levels in early embryos. From left to right, 2-cell, 4-cell, 8-cell, and morula stages. Blue numbers are correlation values. B, linear plots of the proportion of methylation and nonmethylation gene of DEGs in 8 cell and morula stage. C, linear plots of the proportion of methylation and nonmethylation gene of differentially expressed genes at the 4-cell and 8-cell stage. D, DNA methylation distribution of genomic regions in different cytosine contexts for genes with different expression levels. Gene expression levels are divided into six groups from high to low, first is the highest expression; the ordinate is the methylation level; the abscissa is the gene body and the upstream and downstream 2 Kb domains. From left to right, 2-cell, 4-cell, 8-cell, and morula stages. E, violin plot of promoter methylation level distribution of DEGs. F, the line graphs show the methylation level of the repeat elements and the expression level of the repeat elements, respectively. G, venn diagram of DEGs and DMRs in adjacent consecutive periods based on upregulated and downregulated genes, from left to right, Yak4cell-versus-Yak2cell, Yak8cell-versus-Yak4cell, YakMorula-versus-Yak8cell; network diagram for differentially methylated regions of genes and DEGs that are negatively correlated in the EGA period, that is, methylation-driven genes. H, heat map of promoter methylation levels and expression levels of DMGs at different stages. I, heatmap of methylation (methylation) and expression levels (FPKM) of 26 imprinted genes. Blue and red gene names represent paternal and maternal genes, respectively. DEGs, differentially expressed genes; DMG, DMR-related genes; DMRs, differentially methylated regions; EGA, embryonic genome activation; FPKM, fragments per kilobase of exon per million mapped fragments.

Figure 9

Validation of DNA methylation and gene expression data.A–F, methylation levels of six DMR-related genes validated by bisulfite sequencing PCR analysis between 4-cell, 8-cell, and morula. G–K, validation of transcriptome sequencing results by qRT-PCR. Circles on each horizontal chain represent CpG sites in the study sequence, whereas white and black circles represent nonmethylated and methylated CpG, respectively. Each horizontal strand represents a sequenced clone. The circles in the BS-PCR results indicate DNA methylation levels automatically generated by the BiQ Analyzer. In each sample, the number of methylation sites divided by the total number of sites indicates the level of gene methylation. BS-PCR and qRT-PCR were repeated three times for each sample. Different letters indicate samples that differ significantly (p < 0.05). BS, bisulfite; DMR, differentially methylated region; qRT-PCR, quantitative real-time PCR.

Figure 10

Effects of HDAC1/2 and USP7 on the development of yak embryos.A and B, expression analysis of HDAC1/2 in different tissues of yak. C, expression analysis of USP7 in different tissues of yak. D, expression profile of USP7 in preimplantation embryonic development in the yak. E, representative photos in the bright fields on day 5 and day 7 after fertilization of yak embryos treated with DMSO or FK228. The scale bar represents 100 μm. F and G, scatter plots show the total cell number, ICM cell number, and TE cell number of blastocysts in the FK228-treated and P5091-treated groups and the corresponding control groups, respectively. H and I, immunofluorescence (IF) staining of CDX2 in yak blastocysts after DMSO or FK228 treatment. The scale bar represents 50 μm. J and K, IF of CDX2 in yak blastocysts after DMSO or P5091 treatment. The scale bar represents 50 μm. L and M, IF of 5 mC and 5hmC in yak blastocysts after DMSO or P5091 treatment. The scale bar represents 50 μm. (n = 3; 10–20 embryos per group per replicate). Data shown as Means ± SEM. ∗ p < 0.05, ∗∗ p < 0.01. DMSO, dimethyl sulfoxide; HDAC, histone deacetylases; ICM, inner cell mass; TE, trophectoderm.