Zscan4c activates endogenous retrovirus MERVL and cleavage embryo genes

Abstract

Endogenous retroviruses (ERVs) contribute to ∼10 percent of the mouse genome. They are often silenced in differentiated somatic cells but differentially expressed at various embryonic developmental stages. A minority of mouse embryonic stem cells (ESCs), like 2-cell cleavage embryos, highly express ERV MERVL. However, the role of ERVs and mechanism of their activation in these cells are still poorly understood. In this study, we investigated the regulation and function of the stage-specific expressed ERVs, with a particular focus on the totipotency marker MT2/MERVL. We show that the transcription factor Zscan4c functions as an activator of MT2/MERVL and 2-cell/4-cell embryo genes. Zinc finger domains of Zscan4c play an important role in this process. In addition, Zscan4c interacts with MT2 and regulates MT2-nearby 2-cell/4-cell genes through promoting enhancer activity of MT2. Furthermore, MT2 activation is accompanied by enhanced H3K4me1, H3K27ac, and H3K14ac deposition on MT2. Zscan4c also interacts with GBAF chromatin remodelling complex through SCAN domain to further activate MT2 enhancer activity. Taken together, we delineate a previously unrecognized regulatory axis that Zscan4c interacts with and activates MT2/MERVL loci and their nearby genes through epigenetic regulation.

Figure 1.

Analysis of ERV expression dynamics in early embryogenesis by WGCNA. (A) Principal component analysis of embryo samples according to mRNA expression of protein coding genes and ERV expression. The selected ERV loci were separated from protein coding genes and were filtered by low expression (mean cpm < 1). Var., variation. (B) Hierarchical cluster of ERV modules using eigengene values. ERV modules were identified by WGCNA. Only modules with embryonic expression specificity in both datasets (28, 29) were presented. Troph, trophectoderm; ICM, inner cell mass. (C) The ERV repeat family composition of 6 stage-specified modules presented in percentage of all enriched ERVs in each module. (D) WGCNA enrichment of ERVs in six most stage-specified modules. P values were determined by Fisher's exact test.

Figure 2.

The role of Zscan4c in activating of MT2/MERVL expression. (A) Expression of MT2/MERVL, Zscan4c and Dux during early embryogenesis according to published single cell RNA-seq data (28). (B) Expression of total Zscan4 gene cluster members and Dux in MERVL+ and MERVL- ESCs according to published RNA-seq data (17). (C) qPCR analysis of the expression of Zscan4c and retroelements after Zscan4c depletion in mESCs. P1, primer set1; P2, primer set 2. Biological-triplicate data (n = 3 dishes) are presented as mean ± s.e.m. (D) Flow cytometry analysis of MERVL-gag+ population of ESCs expressing control shRNA or Zscan4c shRNA. Biological triplicate data (n = 3 dishes) are presented as mean ± s.e.m. (E) qPCR analysis of the expression of Zscan4c and retroelements after Zscan4c overexpression in mESCs. P1, primer set1; P2, primer set 2. Biological triplicate data (n = 3 dishes) are presented as mean ± s.e.m. (F) Flow cytometry analysis of MERVL-gag+ population of ESCs after Zscan4c overexpression. Biological triplicate data (n = 3 dishes) are presented as mean ± s.e.m. (G) Expression of zinc finger domains of Zscan4c and full length Zscan4c determined by qPCR after overexpression of Zscan4c zinc finger domains. (H) Expression of pluripotency markers (Prdm14, Sox2 and Nanog) and MT2/MERVL determined by qPCR after overexpression of Zscan4c zinc finger domains. P1, primer set1; P2, primer set 2. Biological-triplicate data (n = 3 extracts) are presented as mean ± s.e.m. OE: overexpression.

Figure 3.

Transcriptome analysis of Zscan4c-regulated retroelements and genes after Zscan4c overexpression or depletion. (A, B) Scatter diagrams of TE expression change after Zscan4c overexpression (OE) (A) and depletion by shRNA (B). TEtranscript results were used to plot these diagrams. Coloured dots indicate retroelements with significant expression change (Wald test, FDR adjusted P value < 0.05). Triangles represent ERVs with fold change >3. (C) Expression heatmap of selected retroelements after Zscan4c depletion and Zscan4c overexpression. Those selected TEs were upregulated after Zscan4c overexpression and downregulated after Zscan4c depletion (fold change > 1.5; Wald test, FDR adjusted P < 0.05). (D) Gene set enrichment analysis (GSEA) indicating that upregulated genes after Zscan4c overexpression were highly enriched in the 2-cell/4-cell embryo gene set. Red, upregulated genes; blue, downregulated genes. (E) qPCR analysis of the expression of 2-cell/4-cell embryo genes after Zscan4c overexpression in ESCs. (F) qPCR analysis of the expression of 2-cell/4-cell embryo genes after overexpression of Zscan4c zinc finger domains in ESCs. Biological-triplicate data (n = 3 extracts) are presented as mean ± s.e.m. Expression level of the control sample was set as 1. OE, overexpression.

Figure 4.

Genome-wide binding profile of overexpressed Zscan4c on ESC genome. (A) Locations of Zscan4c peaks relative to the nearest transcription units (Promoter, 2 kb around transcriptional start sites; 5′ proximal, 2–10 kb upstream of gene; 5′ distal, 10–100 kb upstream of gene; 3′ proximal, 0–10 kb downstream of gene; 3′ distal, 10–100 kb downstream of gene; Gene desert, > 100 kb away from the nearest gene). (B) Zscan4c-recognized DNA motifs identified by MEME-ChIP and enriched on MT2. MEME Suite was used to determine E value of motif. (C) Prediction of activating/repressive function of Zscan4c by BETA according to Zscan4c ChIP-seq data and RNA-seq data after Zscan4c overexpression. Red line, upregulated genes; purple line, downregulated genes; dashed line, the non-differentially expressed genes as background. Genes were cumulated by the rank on the basis of the regulatory potential score from high to low. P values that represented the significance of the distributions of upregulated or downregulated genes were compared with the non-differentially expressed genes by the Kolmogorov–Smirnov test. (D) Genome-wide association of Zscan4c regulatory potential score with gene expression change after Zscan4c overexpression. Regulatory potential score was predicted by BETA minus function. Genes were ranked by expression in log₂ (ratio of Zscan4c OE/Control OE) on heatmap on the left. The line plot (right) shows the corresponding BETA scores (moving averages) of genes. (E) ChIP-qPCR analysis of Zscan4c binding on gene promoters. ChIP-qPCR data was normalized to input and that of the control region. Biological-triplicate data (n = 3 extracts) are presented as mean ± s.e.m. (F) ChIP-seq enrichment heatmap of Zscan4c, H3K27ac, H3K4me1, H3K14ac, H3K4me3 signal on Zscan4c binding peaks. The ChIP-seq signal was calculated as the ratio of normalized reads relative to input.

Figure 5.

Zscan4c interacts with MT2 and regulates nearby genes. (A) ChIP-qPCR analysis of Zscan4c binding on different retroelements. ChIP-qPCR data was normalized to input and that of LINE1. Biological triplicate data (n = 3 extracts) are presented as mean ± s.e.m. (B, C) Zscan4c binding profile plot (B) and enrichment heatmap (C) around the centre of MT2 sequence. The ChIP-seq signal was calculated as the ratio of normalized reads relative to input. (D) Zscan4c binding signal on MT2 and position of Zscan4c motif on MT2. Red sequences highlight the sequence of Zscan4c binding motif identified on MT2. (E) DNA pull-down analysis of the DNA binding activity of Zscan4c. DNA pull-down was done with 293T cells overexpressing HA-tagged Zscan4c. Probe with scrambled sequence was used as control. (F) MT2 nearby genes were enriched in upregulated genes after Zscan4c overexpression. Genes within 10 kb from MT2 locus were considered to be near MT2. Red, upregulated genes; blue, downregulated genes. (G) Overlap between upregulated genes after Zscan4c overexpression and MT2 10 kb nearby genes. P value was calculated by Fisher's exact test. (H) Expression heatmap of overlapped genes between upregulated genes after Zscan4c overexpression and genes within 10 kb from MT2. Troph, trophectoderm; ICM, inner cell mass. Average expression level of each developmental stage was determined and presented as Z-score across different stages. (I) Density plot of Zscan4c ChIP-seq signal and RNA-seq expression of MT2/MERVL consensus sequence from GeneBank (ID: Y12713.1). Density plot was plotted by IGV with ChIP-seq and RNA-seq bam files. Yellow region, ChIP-seq signal on MERVL-int; red region, ChIP-seq signal on MT2. ChIP input was included as control. Value of individual peak was highlighted. (J) Luciferase assay analysis of enhancer activity of MERVL LTR(MT2) after Zscan4c overexpression. MT2+, MT2 at the downstream of Luciferase gene; MT2–, reporter without MT2. Zscan4c OE +, Zscan4c overexpression; Zscan4c OE–, control vector overexpression. Biological-triplicate data (n = 3 extracts) are presented as mean ± s.e.m. (K) Luciferase assay analysis enhancer activity of wild type MT2, MT2 with mutant Zscan4c binding motif, and mutant MT2 with inserted Zscan4c binding motif after Zscan4c overexpression. Biological-triplicate data (n = 3 extracts) are presented as mean ± s.e.m.

Figure 6.

Epigenetic regulation of MT2/MERVL after Zscan4c overexpression. (A–C) H3K4me1 (A), H3K27ac (B), and H3K14ac (C) profile plot around the centre of MT2 sequences in control ESCs (Control OE) and Zscan4c-overexpressed ESCs (Zscan4c OE). The ChIP-seq signal was calculated as the ratio of normalized reads relative to input. (D–F) ChIP-qPCR analysis of H3K4me1 (D), H3K27ac (E) and H3K14ac (F) enrichment on specific locus of MT2 in control ESCs and Zscan4c-overexpressed ESCs. LINE1, MTA and IAPEz were included as control. ChIP-qPCR data was normalized to input and that of LINE1. Biological-triplicate data (n = 3 extracts) are presented as mean ± s.e.m.

Figure 7.

Role of Zscan4-associated GABF complex in MT2/MERVL activation. (A) Western blot analysis of Flag-Zscan4c/Brg1/Brd9 co-immunoprecipitation in Flag-Zscan4c overexpressed ESCs or in control vector overexpressed ESCs. IP, immunoprecipitation. IP was done with anti-Flag antibody. 1.25% input was loaded as control. (B–D) Western blot analysis of Brg1/Bd9 co-immunoprecipitation with HA-zinc finger domains (B), HA-SCAN domain (C), and HA-internal domain (D) of Zscan4c. IP was done with anti-HA antibody and separated by SDS/PAGE. 1.25% input was loaded as control. (E) Enrichment heatmap of Brd9 and Brg1 around the binding regions of Zscan4c. The ChIP-seq signal was calculated as the ratio of normalized reads relative to input. (F) ChIP-qPCR analysis of Brg1 and Brd9 enrichment on MT2 in control ESCs (Control OE) and Zscan4c-overexpressed (Zscan4c OE) ESCs. Biological triplicate data were presented as mean ± s.e.m. (G) Relative Luciferase activity of vector containing MT2 in control ESCs (Control OE) and Zscan4c-overexpressed (Zscan4c OE) ESCs treated with inhibitor of Brd9 (BI7273). DMSO treated samples were included as control. Biological triplicate data were presented as mean ± s.e.m. (H) Schematic plot of Zscan4c function in activating MT2/MERVL and 2-cell/4-cell cleavage genes. In the absence of Zscan4c, MT2 enhancer activity is not activated and MT2 nearby 2-cell/4-cell genes are silenced. In the presence of Zscan4c, Zscan4c binds to MT2, interacts with GBAF complex member Brg1 and Brd9, and augments H3K4me1, H3K27ac, and H3K14ac deposition on MT2/MERVL, thus activates MT2 which drives the expression of MERVL and enhances the expression of 2-cell/4-cell genes.