Cohesin facilitates zygotic genome activation in zebrafish

Zygotic genome activation (ZGA) establishes, for the first time in a zygote, a genome that is competent for transcription (Blythe and Wieschaus, 2015; Fassnacht and Ciosk, 2017; Onichtchouk and Driever, 2016; Pálfy et al., 2017; Svoboda et al., 2015). ZGA involves the transfer of maternal to zygotic control of embryonic development.

Commensurate with ZGA, maternal transcripts must be degraded at the maternal-to-zygotic transition (MZT) (Marco, 2017). In zebrafish, many maternal transcripts are targeted for degradation by miR-430, which is among the few early expressed transcripts in the embryo (Bazzini et al., 2012; Giraldez et al., 2006). Other maternal RNAs are N6-methyladenosine (m6A) modified, and are cleared by a m6A-binding protein, Ythdf2 (Zhao et al., 2017). Unique RNA-binding proteins also play a role in controlling RNA metabolism and turnover during MZT (Despic et al., 2017). Therefore, clearance of maternal RNA is essential for transition to the zygotic transcription program.

Mechanisms regulating both transcriptional activation and transcriptional repression are thought to control ZGA in the early embryo (Joseph et al., 2017; Lee et al., 2013; Onichtchouk and Driever, 2016). Evidence from Xenopus and zebrafish suggests the existence of a titratable, maternally deposited repressor that initially holds the transcription of the zygotic genome in check (Kimelman et al., 1987; Newport and Kirschner, 1982b; Nothias et al., 1995). In Xenopus, ZGA coincides with a dramatic increase in the nucleus-to-cytoplasm (N:C) volume ratio (Jevtić and Levy, 2015; Newport and Kirschner, 1982a). Increasing the N:C ratio by the addition of extra DNA (Newport and Kirschner, 1982b) or by injection of scaffolding proteins (Jevtić and Levy, 2015) accelerates ZGA. An in vitro study in Xenopus egg extracts showed that histones H3 and H4 are strong candidates for the maternal repressor activity (Amodeo et al., 2015); repression of transcription by H3/H4 could be manipulated in vitro by altering the ratio of DNA template to histone quantities alone. In zebrafish, core histones outcompete transcription factors for access to the genome, thereby regulating the onset of transcription (Joseph et al., 2017). These studies suggest that transcription is activated as the histone repressors are titrated out during successive cell divisions. Therefore, up until ZGA, repression mechanisms counteract factors that activate transcription in the embryo.

Transcriptional activation at ZGA appears to involve a combination of ‘pioneer’ transcription factor activity, and a gain in active chromatin modifications. At the zebrafish MZT, distinctive histone modifications appear (Andersen et al., 2013; Vastenhouw and Schier, 2012; Vastenhouw et al., 2010), and nucleosomes become strongly positioned at promoters (Zhang et al., 2014b). Even before ZGA, the zebrafish genome is marked with modified histones (Lindeman et al., 2011) and specific sites of DNA methylation (Jiang et al., 2013; Potok et al., 2013). However, although chromatin modifications can demarcate active regions of transcription, additional factors are usually needed for activation of transcription (Hontelez et al., 2015). Sequence-specific transcription factors operating at ZGA vary between species. In Drosophila, the zinc finger protein Zelda activates many early genes (Harrison et al., 2011; Li et al., 2014). In zebrafish, Nanog, Pou5f3 (also called Oct4) and SoxB1 regulate expression of early zygotic genes (Lee et al., 2013; Leichsenring et al., 2013; Onichtchouk and Driever, 2016). Recently, the DUX family of transcription factors was found to activate zygotic genes in mice (De Iaco et al., 2017; Hendrickson et al., 2017; Whiddon et al., 2017).

Global chromatin structure is also linked to transcription activation; for example, formation of architectural features such as topologically associated domains (TADs) mark the onset of transcription in the mouse embryo (Flyamer et al., 2017; Lu et al., 2016). In Drosophila, chromatin architecture in the form of TADs emerges at ZGA independently of gene transcription (Hug et al., 2017). This is consistent with the idea that genome structure formation precedes transcription (Krijger and de Laat, 2017). Chromatin structure in turn influences the binding of transcription factors and RNA Polymerase II (RNAPII) (Newman and Young, 2010).

Although individual players in ZGA may vary between species, a universal theme is that the spatial organization of chromosomes changes as cells commit to developmental fates (de Wit et al., 2013; Hug et al., 2017; Phillips-Cremins, 2014; Vietri Rudan et al., 2015). Spatial organization of the genome depends in part on the nuclear architectural proteins cohesin and CCCTC-binding factor (CTCF), which contribute to the three-dimensional (3D) organization of chromosomes (Vietri Rudan and Hadjur, 2015), and the formation of DNA loops within TADs (Giorgetti et al., 2014; Hug et al., 2017; Van Bortle et al., 2014; Vietri Rudan et al., 2015). Compartmentalization of active and inactive regions of the genome does not depend on CTCF (Nora et al., 2017) or cohesin (Merkenschlager and Nora, 2016). However, local spatial organization within TADs can facilitate transcription of developmental loci (Ferraiuolo et al., 2010; Narendra et al., 2015; Rousseau et al., 2014), thereby determining cell fate and driving embryo development.

In this study, we asked whether cohesin and CTCF contribute to ZGA. We found that cohesin (but not CTCF) depletion delays ZGA, and that chromosome-bound cohesin spreads from satellite and non-coding DNA to genes when the zygotic genome becomes activated. A fraction of gene-associated cohesion-binding sites are co-occupied by ‘pioneer’ transcription factors Pou5f3 and Sox2, and enriched for active histone marks. We propose that cohesin plays a crucial role in organizing a chromatin structure that is permissive for transcription at ZGA.

As embryos progress through MZT at 3.3 h post-fertilization (hpf), the main wave of zygotic gene transcription is activated (Fig. 1A) (Heyn et al., 2014). Relatively low levels of Rad21 (∼100) and CTCF (∼150) transcripts are present pre-ZGA, with both transcript and protein levels increasing by two- to threefold in the main wave of ZGA (Fig. 1B,C). Genes encoding Rad21 and CTCF are essential for cell survival (Nasmyth and Haering, 2005) (including germ cells), which limits the genetic tools available for their manipulation. We previously bypassed homozygous lethality by using morpholino oligonucleotides (MO) to tightly titrate the levels of Rad21 and CTCF (Marsman et al., 2014; Rhodes et al., 2010; Schuster et al., 2015). Here, we were able to substantially reduce the protein levels of cohesin subunit Rad21 and CTCF in early embryos, in order to assess their effects on zygotic genome activation (Fig. 1, Fig. S1E). Rad21-depleted embryos were rescued by a transcript encoding wild-type Rad21, but not by mutant Rad21 containing the rad21^nz171 nonsense mutation (Horsfield et al., 2007) (Fig. S2).

MOs injected at the one-cell stage reduced protein levels of Rad21 and CTCF by 40-80%, even pre-ZGA (Fig. 1D, Fig. S1). By the 4.5 hpf ‘dome’ stage, Rad21-depleted embryos had just slightly fewer cells than wild type, although the difference was not significant (P=0.1138, unpaired t-test) (Fig. S3A). Analysis of cell cycle status by flow cytometry showed that both wild-type and Rad21-depleted embryos had a large majority of cells with 2N DNA content, although the Rad21-depleted embryos had an increase in cells with 4N DNA content, from 7% in wild type to 10% with Rad21 depletion (Fig. S3B). Surviving CTCF-depleted embryos displayed similar cell cycle profiles to wild type (data not shown). However, many CTCF MO-injected embryos died pre-MZT, suggesting that survivors had sub-threshold CTCF depletion.

To determine the effects of Rad21 and CTCF depletion on zygotic transcription, we used RNA-seq to analyze the transcriptome of untreated embryos (referred to here as ‘wild type’) and embryos treated with Rad21- or CTCF-targeting MOs. Five developmental stages were analyzed spanning pre-MZT (2.5 hpf), MZT (3.3 hpf) and post-MZT (4.5 and 5.3 hpf) stages up to the tailbud stage (10 hpf). The sample-to-sample distances between the expression profiles were calculated (R package DESeq2) to cluster time points and treatments. A graphical representation of the sample-to-sample distance is shown in a principal component (PCA) plot in Fig. 2A. We found that Rad21 depletion results in a complement of transcripts that appear delayed in developmental timing relative to wild type at the post-MZT stages of 4.5 and 5.3 hpf (PC1, Fig. 2A). By contrast, profiles from CTCF-depleted embryos cluster similarly to wild-type embryos from the same stage (Fig. 2A).

We identified differentially represented transcripts between wild-type, Rad21- and CTCF-depleted embryos at each of the five stages (Table S1). Rad21 and CTCF depletion most robustly affected transcript levels post ZGA at 4.5 hpf and 5.3 hpf (Fig. 2B). We next annotated (Lee et al., 2013) the origin of the differentially represented transcripts (maternal, weakly maternal or zygotic) in Rad21-depleted embryos and CTCF-depleted embryos compared with wild type at the 4.5 hpf ‘dome stage’. We found that, upon Rad21 depletion, 3285 differentially represented transcripts (FDR=0.05) were both maternal and zygotic, with maternal transcripts more abundant relative to wild-type and zygotic transcripts under-represented (Fig. 2C, Table S1), suggesting ZGA is delayed. Following CTCF depletion, there were 888 differentially represented transcripts, almost fourfold fewer than observed upon Rad21 depletion (FDR=0.05) (Fig. 2D, Table S1).

We then plotted the expression levels of differentially represented transcripts identified from dome stage at all time points sampled. Transcripts that were under-represented in Rad21-depleted embryos normally increase over developmental time in wild type, and transcripts that were over-represented upon Rad21 depletion are reduced over time in wild-type embryos (Fig. 3A,B), with significantly more transcripts affected when compared with CTCF depletion (Fig. S4). CTCF-depleted embryos showed a similar trajectory of differentially represented transcripts over developmental time (Fig. S4A,B). Following Rad21 depletion, delay in the expression of individual zygotic genes was confirmed by quantitative PCR (Fig. 3C). Furthermore, expression of these genes was similarly delayed by depletion of a second cohesin subunit, Smc3 (Fig. 3C, Fig. S1E), suggesting that the effect of Rad21 depletion on ZGA is mediated through abolition of cohesin complex function. We conclude that even partial depletion of cohesin causes a delay in zygotic genome activation.

While transcripts that are differentially represented upon Rad21 depletion were assignable to functional pathways (Fig. 4), transcripts responding to CTCF depletion were not. At the 4.5 hpf ‘dome’ stage, transcripts under-represented upon Rad21 depletion are involved in ribosome assembly, translation and RNA metabolism functions (Fig. 4A, Table S2). Over-represented transcripts reflect the maternal RNA landscape and are involved in energy systems and mitochondrial functions (Fig. 4B, Table S2).

The RNA-seq data indicate that Rad21 (but not CTCF) depletion led to a delay in degradation of maternal mRNAs in combination with a delay in activation of zygotic genes, when compared with stage-matched embryos of equivalent morphology and cell number to wild type. Overall, our data suggest that cohesin is necessary for the timely transition to, and promotion of, maternal-to-zygotic transcription programs.

Considering the importance of Rad21/cohesin for progression to the zygotic transcription program (Figs 2–4), we decided to further investigate Rad21 function during ZGA. To determine the distribution of Rad21 on chromosomes in early development, we conducted chromatin immunoprecipitation followed by high-throughput sequencing (ChIP-seq) in wild-type embryos at 2.5, 4.5 and 10 hpf with custom antibodies against zebrafish Rad21 (Rhodes et al., 2010) (Fig. S5). At 2.5 hpf pre-ZGA, 2011 enriched Rad21 peaks were detected on chromosomes. After ZGA, there was significant recruitment of Rad21 to chromosomes that increased over developmental time (Fig. 5A). By the 4.5 hpf ‘dome’ stage, wild-type embryos had accumulated 7144 significant Rad21 binding peaks, and by 10 hpf, there were 18,075 peaks in total. During ZGA and early development, Rad21 peak distribution shifts closer to the transcription start sites (TSS) of genes (Fig. 5A). The data suggest that pre-ZGA, Rad21/cohesin binds to few loci and is mostly excluded from genes, whereas post-ZGA, Rad21 binding accumulates at gene-dense regions. We performed ATAC-seq (Buenrostro et al., 2015) on wild-type embryos at 2.5 hpf to determine whether Rad21 binds to open chromatin regions at pre-ZGA. About half the accessible chromatin sites at 2.5 hpf also recruit Rad21 (Fisher's exact test, right tail: P<1.00⁻²⁰) (Fig. 5B, Table S4), indicating that a subset of cohesin binding sites is located in the very few regions, 337 in total, that are accessible pre-ZGA. About 90% of the overlapping accessible regions correspond with satellite DNAs (Table S3; specific examples are shown in Fig. S6).

The remarkable redistribution of Rad21 binding to genes post-ZGA is exemplified by its recruitment at chromosome 4 (Fig. 5C). The long arm of chromosome 4 is gene-poor, has extensive heterochromatin and replicates late. High densities of 5S ribosomal DNA (rDNA), small nuclear RNAs (snRNAs), half of all tRNAs and 30% of all zinc-finger domain genes are present on the long arm of chromosome 4 (Howe et al., 2013). Prior to ZGA, many of these loci were enriched for Rad21 (right side), whereas the RNAPII gene-rich region of chromosome 4 (left side) excluded Rad21 binding (Fig. 5C). Post-ZGA, Rad21 binding became increasingly enriched at the RNAPII gene-rich region of chromosome 4, and some of its pre-ZGA binding sites were lost (Fig. 5C).

Genes that recruit Rad21 pre-ZGA (within a 20 kb window) included hsp70l, sox2, gata2a and the miR-430 complex (Fig. S7), and multiple zinc finger domain-encoding proteins located on chromosome 4 (Fig. S8). Transcripts from miR-430 and zinc-finger domain-encoding genes are expressed as early as the 64-cell stage, prior to the main wave of zygotic transcription (Heyn et al., 2014). The mature miR-430 microRNAs mark a substantial amount of maternally deposited transcripts for degradation (Giraldez et al., 2006). Interestingly, many zinc finger-encoding genes marked by Rad21 binding are not expressed until post-ZGA (Fig. S8).

Besides association with regions on the long arm of chromosome 4, 41% (824/2011) of the Rad21 peaks were found at satellite elements (satDNAs) located at pericentromeric regions of the genome (Fig. 5D,E; Fig. S6; Table 1). Further classification of these elements shows that BRSATI and SAT-1 were among the most- enriched members of satDNAs found at Rad21 sites pre-ZGA, representing over 70% of satDNAs identified. satDNAs represent less than 1% of the genome in zebrafish and are therefore significantly enriched (P value <1.0⁻²⁰, Fisher's exact test) in the 2.5 hpf Rad21 peaks, whereas DNA transposons are relatively abundant; accounting for 33% of the genome and are not significantly enriched. Long terminal repeats (LTRs) and rRNAs are also significantly enriched in pre ZGA Rad21 peaks (Table 1).

After ZGA, there was significant recruitment of Rad21 to chromosomes that increased over developmental time (Fig. 5A). By the 4.5 hpf ‘dome’ stage, wild-type embryos had accumulated 7144 significant Rad21 binding peaks, including ∼3000 that were gene associated; by 10 hpf, there were 18,075 peaks in total, 5937 of which were gene associated (Figs 5A, 6A, Table S3). Rad21 binding was significantly over-represented in coding regions after ZGA. Furthermore, Rad21 binding was particularly enriched at promoters and 5′ untranslated regions (5′ UTRs), as well as in exons, transcription termination sites (TTS) and 3′ UTRs (Fig. 6A).

Overall, 12% (293/2371) of over-represented transcripts and 15% (179/1185) of under-represented transcripts at dome stage were derived from genes that recruited Rad21 (Fig. 6B). However, the small overlap of differential expression with bound genes was non-significant. We used k-means clustering (k=2) to visualize Rad21 binding profiles over two subsequent developmental stages (4.5 and 10 hpf). About half of the regulated genes that contained Rad21 binding (58% for under-represented transcripts and 52% for over-represented) had lost that binding by 10 hpf (Fig. 6C,D), indicating that Rad21 is likely to be specifically associated with those genes during ZGA, and potentially involved in their direct regulation at that time. Following Rad21 depletion, when compared with over-represented transcripts, genes with under-represented transcripts at 4.5 hpf had higher transcription levels in wild-type embryos by 10 hpf, irrespective of Rad21 binding (Fig. 6E). Only a small fraction (5-8%) of genes with over-represented transcripts at 4.5 hpf also showed altered expression at 10 hpf. Genes found to be downregulated at 4.5 hpf were more likely to also be differentially expressed at 10 hpf (14%), but there was no difference between genes that gain or lose Rad21 binding at 10 hpf (Fig. 6F, clusters I and II, respectively). This indicates that a small subset of genes bound by Rad21 during ZGA are affected later in development by Rad21 depletion. However, although it is likely that some of the bound genes may be regulated directly, a larger fraction appears to be regulated indirectly.

Our results indicate that Rad21 is present at repetitive sequences and ncRNA genes prior to ZGA, with a transition to RNAPII genes at ZGA, once transcription starts. The marked enrichment of Rad21 at genes through developmental time suggests that cohesin may facilitate their expression. However, because many more genes are regulated by Rad21 depletion than are bound by Rad21, it is unlikely that direct gene regulation by cohesin explains the delay in ZGA.

To further investigate a possible role for Rad21 in ZGA, we sought to determine whether Rad21 binding coincides with other hallmarks of gene activation, including H3K4me1 and H3K27ac enhancer modifications, H3K4me3 marks associated with active gene promoters, and H3K27me3 modification of polycomb-repressed genes (Vastenhouw and Schier, 2012). For this analysis, we surveyed defined regions centered on Rad21 binding sites for enrichment of these modified histones by comparing Rad21 ChIP-seq data with publicly available histone ChIP-seq data (Bogdanovic et al., 2012; Zhang et al., 2014b) at peri-ZGA time points (Fig. 7A). Forty-eight percent of the Rad21 peaks (3416/7144) overlapped with at least one of the enhancer- and promoter-associated marks, H3K4me1, H3K4me3 and H3K27ac (Fisher's exact test, right tail, P<1.0⁻²⁰) (Fig. 7A, Table S4). A smaller set of 244 peaks was significantly associated with H3K27me3 (Fisher's exact test, right tail, P=2.72⁻³⁴). Fewer than half of overlapping peaks were at TSSs (Table S4). Therefore, cohesin binding significantly coincides with histone marks that are associated with active chromatin at ZGA.

The transcription factors Nanog-like, Pou5f3 and the SoxB1 family (including Sox2, Sox3, Sox19a and Sox19b) are homologs of mammalian pluripotency factors, and are thought to act as pioneering factors in zebrafish ZGA (Lee et al., 2013; Leichsenring et al., 2013). Because Rad21 depletion delayed ZGA (Fig. 2), we were interested to determine whether Rad21/cohesin binding coincides with genomic locations of these activators of early gene expression. Publicly available ChIP-seq data for Nanog-like (Xu et al., 2012), Pou5f3 and Sox2 (Leichsenring et al., 2013) were obtained (Table S4) and compared with Rad21 binding at 4.5 hpf. There was a small but significant overlap between Rad21 and pluripotency factor binding sites (Table S4 and Fig. 7B-D). Regions with overlap of Rad21 and pluripotency factors were enriched for developmental, chromatin assembly and pattern specification ontologies (Fig. 7B-D). The coincidence of a subset of cohesin protein binding with these known transcriptional activators suggests that cohesin may be involved in regulating selected Sox2 and Pou5f3 targets.

The combined data above point to a generalized role for Rad21 in early zygotic transcription. Given the known role of cohesin in the local spatial organization of chromatin (Merkenschlager and Nora, 2016), we addressed the possibility that Rad21/cohesin might contribute to ZGA through global organization of chromatin architecture.

We used antibodies to nucleolin and RNAPII to visualize nucleoli and RNAPII clustering, respectively, immediately post-ZGA. Immunofluorescence analysis in 4.5 hpf ‘dome’ stage embryos revealed the formation of nucleoli (Fig. 8A,A′) and discrete RNAPII clusters that may represent transcription foci (Fig. 8D,D′). Strikingly, depletion of Rad21 severely and significantly disrupted the formation of nucleoli (Fig. 8B-C - P <4.7⁻¹⁴) and RNAPII clusters (Fig. 8E-F - P<8.8⁻⁷) at this developmental stage, with these markers exhibiting a more fragmented appearance.

Ribosomal DNA (rDNA) is contained within the nucleolus and active rDNA interacts with nucleolin (Cong et al., 2012). In support of a direct role for Rad21 in nucleolar organization, we found that Rad21 is enriched near 5S ribosomal DNA repeats, on chromosomes 4 (Fig. 5C), 18 and 22 (Table S3).

Thus, depletion of Rad21 dramatically affects nuclear organization by the time of ZGA. Our results further suggest that Rad21 recruitment to genes at this crucial developmental stage influences the formation of RNAPII foci that could represent early sites of transcription. Transcripts from many genes are affected by Rad21 depletion in a manner consistent with an overall delay in ZGA (Figs 2 and 3), although few of these genes are directly bound by Rad21 (Fig. 6). Global disruption of chromosome organization by Rad21 depletion provides a possible mechanism for this observation.

Altogether, our results point to a global role for Rad21/cohesin in facilitating ZGA in zebrafish embryos. Rad21 locates to active regions of the genome, including genes expressed at ZGA, whereas Rad21 depletion interferes with gene expression. Rad21 depletion also affects nuclear integrity and RNAPII clusters, raising the possibility that cohesin plays a role in organizing a chromatin structure that is permissive for transcription around the time of ZGA.

Owing to their important combinatorial roles in genome organization (Vietri Rudan and Hadjur, 2015), we expected Rad21/cohesin and CTCF depletion to have similar effects on activation of the zygotic genome. We were surprised to find that, upon depletion, their effects were quite different. We previously showed that a modest depletion of Rad21 by morpholino (Schuster et al., 2015) or mutation (Horsfield et al., 2007) can have striking effects on the transcription of specific genes. By contrast, CTCF had to be dramatically depleted to affect transcription, and doing so resulted in high levels of mortality (data not shown) (Marsman et al., 2014). CTCF is essential to the integrity of the nucleus; our data suggest that a small amount of CTCF may be sufficient for this function and that CTCF depletion beyond this level is lethal in embryos in which cells proliferate rapidly. Consistent with this, maternal and zygotic depletion of CTCF leads to apoptosis and is lethal at preimplantation stages in mice (Moore et al., 2012; Wan et al., 2008). In contrast to CTCF, partial Rad21 depletion generated multiple robust biological effects at zebrafish ZGA, and this tractability encouraged us to focus our study on Rad21/cohesin.

Rad21 depletion dramatically altered the transcript complement in embryos just post-ZGA, reflecting an overall delay in transition from the maternal to zygotic transcription program. Among the top downregulated gene ontology categories were ribosome assembly, RNA processing and translation. These processes are also compromised by disrupting cohesin function in yeast and mammalian cells (Bose et al., 2012). Xu et al. previously demonstrated that translational defects in zebrafish and mammalian cell cohesin mutants were chemically rescued by L-Leucine stimulation of the TOR pathway (Xu et al., 2015, 2016). Our data are consistent with these observations that translational mechanisms require normal cohesin function. Moreover, cell proliferation dominates early zebrafish development and requires high levels of translation, consistent with the emergence of these biological pathways as the most significantly affected by cohesin depletion just post-ZGA. Compromising this aspect of the normal gene expression program will almost certainly affect embryogenesis, consistent with mutations in cohesin causing the human developmental disorder Cornelia de Lange syndrome (CdLS).

Although cell proliferation is central to early development, in Drosophila and zebrafish, ZGA is independent of, or upstream of, cell cycle number and checkpoint regulators (Blythe and Wieschaus, 2015; Zhang et al., 2014a, 2017). Consistent with these observations, we found that delay in ZGA occurred in Rad21-depleted embryos that had the same number of cells as controls. However, in zebrafish, ZGA does reflect replication timing in the early embryo (Siefert et al., 2017), and we cannot rule out the possibility that replication timing is affected in our experiments.

Rad21 was generally excluded from genes pre-ZGA with some notable exceptions. Genes that recruited Rad21 pre-ZGA included hsp70l, sox2, gata2a and the miR-430 complex. Of these, gata2a and miR-430 are expressed pre-ZGA (Heyn et al., 2014), and cohesin binding to these locations was reduced once embryos transited through ZGA. Interestingly, miR-430 is responsible for targeting maternal transcripts for clearance (Bazzini et al., 2012; Giraldez et al., 2006), and this raises the possibility that a proportion of maternal transcripts with delayed degradation upon Rad21 depletion could be accounted for by dysregulated miR-430.

Other genes that recruit Rad21 pre-ZGA are generally not expressed at that time. For example, sox2 mRNA is maternally provided, and is involved in transcription of early-expressed zygotic genes in zebrafish (Lee et al., 2013). The zygotic sox2 gene is expressed post-ZGA (Heyn et al., 2014). Significantly, the timing of sox2 expression post-ZGA coincides with a redistribution of Rad21 peaks at the sox2 gene. In addition, Rad21 is recruited to several zinc-finger protein-encoding genes pre-ZGA that are expressed at post-ZGA stages of development. Recruitment of Rad21 to silent loci pre-ZGA could indicate that cohesin has a function there, perhaps to mark their later expression.

Prior to ZGA, cohesin is highly enriched at satellite DNAs found at pericentromeric regions, which represent fewer than 1% of the genome, as well as ncRNA genes. Various satellite sequences in somatic cells are packaged into constitutive heterochromatin, which is characterized by high compaction, enrichment of repressive histone modifications, transcriptional quiescence and late replication. Most of these attributes are absent in pre-ZGA embryos, and the satellite sequences seem to take on these features successively as the embryo develops (Borsos and Torres-Padilla, 2016). ATAC-seq indicates that the satellite DNAs are highly accessible pre-ZGA. It is possible that cohesin is sequestered there merely because this chromatin is open, and thus satellite DNA serves to keep cohesin away from RNAPII genes prior to genome activation. Cohesin depletion results in organizational changes of nucleoli, which could interfere with satellite-dependent heterochromatin formation at ZGA.

At post-ZGA stages, thousands of genomic locations recruit Rad21, and markedly include genic features such as promoters, TSSs, termination sites, and 3′ and 5′ UTRs and exons. Subsequently, Rad21 is increasingly enriched at TSSs, notably at sites co-enriched in histone modifications that are indicative of active promoters and enhancers (Vastenhouw and Schier, 2012; Vastenhouw et al., 2010). Gene-associated Rad21 significantly overlapped with occupancy of the pluripotency factors Pou5f3 and Sox2 at similar time points. This raises the possibility that the pluripotency factors pioneer sites of zygotic transcription (Lee et al., 2013) and recruit cohesin at a subset of these to keep these regions in an ‘open’ configuration. In support of this idea, nucleosome density increases following cohesin loss (Yan et al., 2013), suggesting that cohesin acts to keep chromatin open. Moreover, nucleosome organization is a key feature of ZGA; in zebrafish, nucleosomes are strongly positioned at promoters (Zhang et al., 2014b) at a stage that is coincident with cohesin binding.

Cohesin may also operate at a nuclear structural level to regulate ZGA. Consistent with previous observations in yeast and human cells (Bose et al., 2012; Harris et al., 2014), Rad21/cohesin is essential for the formation of nucleoli in post-ZGA zebrafish embryos. Loss of nucleoli could have global effects on zygotic transcription and translation, as was observed in this study and others (Xu et al., 2015). In addition, reduction in Rad21 just post-ZGA resulted in dispersion of RNAPII foci that could represent transcription factories. Loss of chromosome architecture owing to cohesin depletion could lead to an inability to assemble a transcription-competent genome structure. A combination of the factors described above could lead to global dysregulation of the zygotic transcription program, and these factors are indicative of roles for cohesin at multiple levels at ZGA (Fig. 9).

A spectrum of multifactorial human developmental disorders known as the ‘cohesinopathies’ arise from mutations in cohesin regulators or cohesin subunits (Ball et al., 2014; Bose and Gerton, 2010; Horsfield et al., 2012; Skibbens et al., 2013). Our study raises the possibility that germline cohesinopathy mutations could lead to global alteration of the zygotic transcription program right from the start of development, perhaps explaining the diversity of phenotypes observed in individuals with cohesinopathy.

Zebrafish were maintained under standard conditions (Westerfield, 1995). The University of Otago Animal Ethics Committee approved all zebrafish research. Morpholino oligonucleotides (MOs) were obtained from GeneTools and diluted in water. MO sequences were: Rad21, 5′-AGGACGAAGTGGGCGTAAAACATTG-3′; CTCF, 5′-CATGGGTAATACCTACATTGGTTAA-3′ (targeting the ATG), 5′-CCAAAACAGATCACAAACCTGAAAG-3′ (targeting the splice site of intron 2); and Smc3 5′-TGTACATGGCGGTTTATGC-3′ (targeting the ATG), as described previously (Marsman et al., 2014; Rhodes et al., 2010; Schuster et al., 2015). For microinjection, 1 nl containing 1.0 pmol (for embryos up to MZT) or 0.25-0.5 pmol (for embryos grown post-MZT) of each MO was injected at the one-cell stage. CTCF MOs were combined in an equimolar ratio. For mRNA rescue of the Rad21 MO, embryos were injected with MO from one needle and rescue mRNA (200 pg) from a second needle. Mutant rad21^nz171 mRNA (Horsfield et al., 2007) was used as a control.

Wild-type embryos were collected at the one-cell stage, synchronized and either morpholino injected or kept as a control and allowed to develop to the desired stage (2.5, 3.3, 4.5, 5.3, 10 hpf) at 28°C. Three biological replicates each containing total RNA from 100 pooled embryos were isolated using the NucleoSpin RNAII Kit (Macherey-Nagel). The quality of the RNA was confirmed using the Agilent 2100 Bioanalyzer; in all samples, RIN>9.

Triplicate RNA samples from morphologically stage-matched embryos were sequenced to compare expression profiles over time. Strand-specific libraries were prepared using the TruSeq stranded total RNA-ribozero kit (Illumina) and 100 bp paired-end sequencing was performed to depth of 10 million reads per library on an Illumina HiSeq 2000. On average, 19 million 100 bp paired-end reads per library were generated. These were then adapter and quality trimmed using cutadapt (Martin, 2011) and SolexaQA (Cox et al., 2010). Each sequencing data set was independently mapped to the zebrafish genome with a bowtie2 index generated from Danio_rerio.Zv9.70 (Ensembl) downloaded from Illumina's iGenomes collection. Zebrafish genome danRer7[Zv9] was used to provide known transcript annotations from Ensembl using TopHat2 (version 2.0.9) (Kim et al., 2013) with the following options: ‘tophat2 --GTF genes.gtf --library-type fr-firststrand -p 24 --mate-inner-dist -8 --mate-std-dev 6 zv9’ (on average, 75.38% reads mapped uniquely to the genome). Transcriptomes were assembled with Cufflinks (version 2.2.0) (Trapnell et al., 2010) using options: ‘cufflinks -p 32 --GTF genes.gtf’ and differential expression analysis between control and knockdown embryos was performed using Cuffdiff. A FDR corrected P-value of 0.05 was applied as the cut-off to identify differentially regulated transcripts. The R package DESeq2 was used to compare expression profiles over time. The R package clusterProfiler (Yu et al., 2012) was used to identify enriched Gene ontology terms in upregulated and downregulated gene lists using a cut-off of 0.05 FDR corrected P-value.

From RNA-seq data, five candidate genes were selected for confirmation by quantitative PCR following Rad21 and Smc3 depletion. Embryos were collected at four stages (2.5, 3.3, 4.5 and 5.3 hpf), RNA extracted as above and cDNA synthesized (qScript). Quantitative RT-PCR was performed with primers designed to each of the five candidates (Table S5). Primers were designed to span exon-exon junctions to amplify only processed mRNA transcripts. Expression was normalized to the mitochondrial gene nd3 (Table S5).

Anti-Rad21 (Rhodes et al., 2010) and anti-CTCF (Marsman et al., 2014) were raised in rabbit against a 15 amino acid peptide of each of the zebrafish proteins (GenScript). Commercial primary antibodies were: mouse anti-γ-tubulin (T5326, Sigma-Aldrich; 1:5000), anti-nucleolin (ab22758, Abcam; 1:1000), anti-RNA polymerase II CTD repeat YSPTSPS (phospho S2) (ab5095, Abcam; 1:1000), anti-SMC3 rabbit mAb (D47B5, Cell Signaling Technology; 1:200). Secondary antibodies were goat anti-rabbit IgG (whole IgG heavy and light chains) (#A-11008, Thermo Fisher Scientific) and IRDye-conjugated antibodies (926-68070 and 926-32211, LiCor).

Following dechorionation and deyolking, zebrafish embryos were lysed in RIPA buffer and equal amounts of protein were separated by electrophoresis on 10% polyacrylamide gels. Proteins were transferred to nitrocellulose (Thermoscientific) and incubated with mouse anti-γ-tubulin (1:5000) and rabbit anti-Rad21 (1:500), secondary antibodies were the IRDye-conjugated antibodies (1:15,000). Blots were visualized with the Odyssey CLx Infrared imaging system (LiCor). Band intensities were quantified using Image Studio 4.0 Software (LiCor).

Chromatin was prepared from two independent collections of pooled embryos (n=2000) for 2.5 hpf stage embryos and n=1000 for 4.5 hpf and 10 hpf embryos as described previously (Lindeman et al., 2009). Briefly, embryos were dechorionated using a syringe with a 21G needle, fixed in 1% formaldehyde for 10 min at room temperature. Fixation was stopped by adding glycine to a final concentration to 0.125 M and incubation on ice for 5 min. Fixed embryos were then washed three times in ice-cold 1× PBS, snap frozen and stored at −80°C until use. After cell lysis, chromatin was sheared to 200-500 bp using a S220 Focused-ultrasonicator (Covaris) with the following settings per cycle: peak power=70, duty factor=5, cycles of bursts=200, time=30 s. Individual cycle numbers were optimized for each stage. Chromatin from pre-MZT embryos needed six cycles of sonication to reach the desired 200-500 bp range, whereas chromatin isolated from 4.5 hpf and 10 hpf stages required 10 cycles. Cell debris was removed by centrifugation. To provide standardized input for each ChIP experiment, chromatin was diluted to A₂₆₀=0.25. For each ChIP, 6 µg of Rad21 antibody per 10 µl dynabeads and 100 µl chromatin were incubated overnight at 4°C. After elution, ChIP DNA and input controls were purified and precipitated with ethanol.

The ThruPLEX DNA-seq Kit (Rubicon Genomics) was used to prepare the 2.5 hpf sample libraries for sequencing. 125 bp paired-end sequencing was performed to a depth of 20-50 million reads per library on an Illumina HiSeq 2500 by New Zealand Genomics. Libraries for the 4.5 hpf and 10 hpf samples were constructed and sequenced at the Beijing Genomics Institute (China), yielding 20 million 50 bp single-end reads per sample. After adapter and quality trimming using cutadapt (Martin, 2011) and SolexaQA (Cox et al., 2010), reads were aligned to the Zv9 genome assembly using bowtie2 (Langmead and Salzberg, 2012) (version 2.2.1.) with default settings. 85% of the raw reads could be aligned except for Rad21 ChIP samples from 2.5 hpf, which had a 40% mapping rate. As an alternative to bowtie2, we used the aligner SHRiMP2 (David et al., 2011), increasing the mapping rate to 79% for the 2.5 hpf Rad21 IP samples. Peak finding and downstream data analysis were performed using HOMER (Nagy et al., 2013) and MACS2 (Zhang et al., 2008). Peaks were defined at a 0.1% estimated false discovery rate. Repetitive elements were obtained from the Repeatmasker database (Tarailo-Graovac and Chen, 2009) and overlapped with Rad21 peaks using Bedtools2 (Quinlan and Hall, 2010). Heat maps were generated using the log₂ ratios of the binned reads comparing ChIP input and IP samples using deepTools2 (Galaxy version 2.2.3.0) bamCompare (Galaxy version 2.2.3), computeMatrix (Galaxy version 2.2.5) and plotHeatmap (Galaxy version 2.2.5) (Ramirez et al., 2016).

The ATAC-seq libraries from zebrafish embryos were prepared as previously described (Buenrostro et al., 2015) with some modification. Embryos were collected at the 256-cell stage and dechorionated using pronase. Yolk was removed using deyolking buffer as previously described (Link et al., 2006). 75,000 cells were used to prepare the libraries. Libraries were pooled and sequenced on an Illumina HiSeq. Data were aligned to Zv9 using bowtie2 and peaks for each replicate were called using MACS2 (Zhang et al., 2008). Peaks identified in both replicates were used for downstream analysis.

Around 100 embryos at 4.5 hpf (wild type and MO injected) were dechorionated and deyolked. Cells were fixed with 100% ethanol overnight. Fxcycle PI/RNase staining solution (Thermo Fisher Scientific) was used to stain DNA in cells. Flow cytometric acquisitions were performed on a FACSCALIBUR (BD). Analyses were performed using FlowJo software (Treestar).

Embryos were fixed, dehydrated and stored in 100% methanol at −20°C. For staining, embryos were rehydrated in methanol/PBT, incubated in 150 mM Tris-HCl (pH 9), followed by heating at 70°C for 15 min (Inoue and Wittbrodt, 2011). Embryos were washed and blocked in 5% sheep serum, 2 mg/ml BSA in PBT. Primary antibodies were used at 1:1000 and the secondary antibody at 1:1000 dilution together with 1:1000 Hoechst stain (Thermo Fisher Scientific) in blocking buffer. Embryos were washed in PBT and stored in DAKO mounting media until image acquisition. Confocal immunofluorescence images were acquired using a confocal microscope (Nikon C2, Nikon) with a CFI Plan Fluor NA 0.3/10× objective and a NA 1.4/60× oil immersion objective.

The Imaris software package with default parameters was used to quantify numbers of nuclei per embryo from z-stacks of Hoechst stained embryos. For quantifying nucleolar integrity and RNAPII foci, a single focal plane was obtained through the center of the nucleus. To quantify the area of immunodetected nucleolin or RNAPII relative to the size of the nucleus, a constant threshold setting in NIS-elements imaging software (Nikon) was used for wild-type and Rad21-depleted embryos. The area of the nucleus outlined by the Hoechst staining defined regions of interest for which the pixel area of nucleolin or RNAPII was measured.

Statistical tests were performed using R (www.r-project.org/). Contingency tables and source data for statistical analyses are provided in Table S4.

The authors would thank Dale Dorsett for advice on ChIP-seq analysis and Noel Jhinku for expert management of the zebrafish facility.