Competition for H2A.Z underlies the developmental impacts of repetitive element de-repression

Throughout the eukaryotic lineage, canonical H2A histones are replaced by variant H2A.Z histones within actively transcribed gene promoters and enhancers (Jin et al., 2009; Ku et al., 2012; Murphy et al., 2018; 2020; Subramanian et al., 2013; 2015; Weber et al., 2010; 2014). This epigenetic alternation is thought to facilitate transcription via improving chromatin competency surrounding transcription initiation sites, and H2A.Z-containing nucleosomes have been demonstrated to function as a platform for assembling additional epigenetic modifications during development, including bivalent chromatin in stem cells (Hickey et al., 2022; Hu et al., 2013; Ku et al., 2012). However, several prior studies have identified areas of H2A.Z enrichment located outside genic regions, including centromeres and repetitive elements (REs) (Belotti et al., 2020; Boyarchuk et al., 2014; Greaves et al., 2007; Hardy et al., 2009; Rangasamy, 2013; Rangasamy et al., 2003; Sarcinella et al., 2007). The function of H2A.Z at these intergenic sites remains poorly understood.

Among intergenic locations where H2A.Z has been identified are ancestral transposons (Belotti et al., 2020; Greaves et al., 2007; Martire and Banaszynski, 2020; Rangasamy, 2013), viral remnants from prior generations which can cause mutagenesis upon reactivation (Bourque et al., 2018; Cosby et al., 2019). Owing to their destructive nature, strongly repressive epigenetic features are established over the majority of REs during early development, including the histone modification tri-methylation at histone 3 lysine 9 (H3K9me3) and DNA methylation (Burton et al., 2020; Guo et al., 2014; Laue et al., 2019; Potok et al., 2013; Ross et al., 2020; Smith et al., 2012; Xu et al., 2022). The presence of these repressive features along with H2A.Z (Fan et al., 2004; Greaves et al., 2007; Rangasamy et al., 2003; Ryan and Tremethick, 2018; Sarcinella et al., 2007; Spracklin et al., 2023) is in stark contrast to the epigenetic patterns occurring over transcriptionally active genic loci, which possess high levels of H2A.Z along with activating histone modifications, such as histone H3 lysine-27 acetylation (H3K27ac) and DNA hypomethylation (Brunelle et al., 2015; Hu et al., 2013). Additional support for H2A.Z functioning in transcriptional control of REs comes from studies of invertebrates and plants. In Drosophila, H2Av (the homolog to H2A.Z) has been shown to localize at highly repetitive loci, including the centromeric heterochromatin (Swaminathan et al., 2005; Zhang and Pugh, 2011), and in Arabidopsis, H2A.W (a related variant of H2A) resides within heterochromatin, contributing to silencing of mobile transposable elements (Osakabe et al., 2021; Yelagandula et al., 2014). Interestingly, loss of the Arabidopsis chromatin remodeler DDM1 causes H2A.W to be replaced by H2A.Z at transposon-rich loci (Jamge et al., 2023), indicating that H2A.Z localization might depend on the presence of other epigenetic features. In vertebrates, it has not been determined whether H2A.Z provides a distinct function at intergenic REs or transposons compared with its function at genic regulatory loci.

Recent studies have demonstrated that REs can sometimes become de-repressed upon exposure to environmental stressors or harmful toxicants (Faulk et al., 2016; Pathak and Feil, 2018), ultimately leading to innate immune induction. Deleterious agents typically operate through epigenetic changes and/or DNA damage (Bollati and Baccarelli, 2010; Hou et al., 2012), followed by transcriptional re-activation of REs (Chen et al., 2021; Miousse et al., 2015). Among environmental toxicants functioning in this manner are bisphenol A, an estrogen agonist found in many plastic products (Dolinoy et al., 2007; Roy et al., 2012; Sowers et al., 2020), and Tris(1,3-dichloroisopropyl)phosphate (TDCIPP), a widely used chemical commonly found in flame retardants and pesticides (Liu et al., 2016; Volz et al., 2016). Genetic mutations causing DNA damage or epigenetic malfunction in humans or model organisms can also prompt this response, including mutations in the kinase ataxia-telangiectasia, mutated (ATM) (Gasser and Raulet, 2006; Gasser et al., 2005; Härtlova et al., 2015; Petersen et al., 2012), RNase H (Amon and Koshland, 2016; Reijns et al., 2012), the deoxynucleoside triphosphohydrolase (dNTPase) SAMHD1 (Coquel et al., 2018), the histone demethylase LSD1 (also known as KDM1A) (Sheng et al., 2018), the DNA methyltransferase DNMT1 (Chernyavskaya et al., 2017; Magnani et al., 2021) and several histone methyltransferases (DiNapoli et al., 2022; Fang et al., 2012; Hansen et al., 2022). Recently, separate pharmacological strategies have been developed exploiting this response to combat carcinogenesis, a process referred to as ‘viral mimicry’ (Chen et al., 2021; Chiappinelli et al., 2015; Roulois et al., 2015). These include the use of decitabine (DAC), an FDA-approved form of the DNA methyltransferase inhibitor 5-aza-deoxycytodine (Kantarjian et al., 2006; Stresemann and Lyko, 2008), and T-3775440, a recently developed inhibitor of the LSD1 histone demethylase, both of which have been used for treating acute myeloid leukemia (Takagi et al., 2017). Effective strategies for viral mimicry are thought to operate through epigenetic de-repression and transcriptional re-activation of endogenous transposons, which in turn causes the innate immune system to ‘sense’ a looming infection. Notably however, several studies have found that activation of viral mimicry can sometimes occur in the absence of significant DNA methylation loss (Bowling et al., 2021; Goel et al., 2017; Lee et al., 2020).

We previously found that incorporation of H2A.Z at developmental stage-specific promoters and enhancers is essential for gene expression control during zebrafish stem cell formation (Hickey et al., 2022; Murphy et al., 2018) and differentiation (Akdogan-Ozdilek et al., 2022). However, mechanisms by which H2A.Z landscapes were established, whether H2A.Z functions in transcription control of REs and whether H2A.Z localization changes occur in response to extrinsic stressors, remain unknown. If indeed H2A.Z is important both for transcriptional control of REs and regulation of lineage-specific gene expression, then de-repression of REs during development would present a unique challenge for organisms. Embryos would thus be obligated to utilize H2A.Z both for response to insult and for establishing proper developmental transcription patterns. To investigate this potential conflict, we exposed developing zebrafish embryos to separate stimuli known to cause de-repression of REs (via viral mimicry), the environmental toxicant TDCIPP and the chemotherapeutic agent DAC. We further assessed impacts on H2A.Z localization, finding that distinct changes occur at promoters compared with changes at REs. Our study reveals that H2A.Z localization patterns are highly dependent on total H2A.Z availability, that REs function as a nuclear sink for sequestration of H2A.Z, and that total H2A.Z protein levels determine embryonic sensitivity to environmental stressor and toxicants.

To establish normal H2A.Z localization patterns in developing zebrafish embryos, we first profiled shield stage (early gastrula) zebrafish embryos using cleave under targets and tagmentation (CUT&Tag) (Kaya-Okur et al., 2019). We found H2A.Z to be enriched primarily at genic regulatory regions (Fig. 1A), and as anticipated, its promoter [transcription start site (TSS)±1 kb] levels were positively correlated with gene expression (Fig. S1A, Spearman correlation 0.53). High levels of H2A.Z at intergenic loci were found at many classes of REs, such as long interspersed nuclear elements (LINEs) and long terminal repeats (LTRs) (Fig. 1A; Fig. S1B). With respect to other epigenetic features, H2A.Z enrichment over REs was somewhat distinct from that at genic regulatory loci, such as CpG islands or gene promoters (Murphy et al., 2018). Unlike genic loci, REs that possessed H2A.Z also harbored moderately high levels of the repressive marks H3K9me3 and DNA methylation (Fig. 1B,C; Fig. S1C,D), as well as H3K27ac, which typically occurs at transcriptionally active loci. The combination of repressive (H3K9me3 and DNA methylation) and active marks (H3K27ac), led us to speculate that H2A.Z-containing REs might be temporarily maintained in a poised transcriptional state, such that they can become active later during development or in response to external stimuli or stressors. In support of this possibility, H2A.Z enrichment over REs correlated with transcription levels (total RNA) (Fig. 1D), further implicating a potential H2A.Z-associated priming mechanism in RE activation.

To investigate H2A.Z function in the context of REs, we utilized separate tools for de-repression, and exposed developing embryos to DAC, a well-known cancer therapeutic agent (Kantarjian et al., 2006), and TDCIPP, an environmental toxicant (Liu et al., 2016; Volz et al., 2016). Minimal exposure to these agents caused phenotypic outcomes in embryos at 6 h post fertilization (hpf) and 12 hpf that were similar for both conditions (Fig. 2A), as in prior studies (Dasgupta et al., 2018; Potok et al., 2013; Ribas et al., 2017; Volz et al., 2016), and genes involved in immune response were upregulated, including interferon-induced protein with tetratricopeptide repeats 11 (ifit11), ifit8 and interferon-stimulated gene 15 (isg15) (Fig. 2B; Fig. S2A). We next assessed genomic H2A.Z localization patterns and found that exposure to either TDCIPP or DAC led to a drastic reduction in H2A.Z levels at regions defined as ‘peaks’ (see Materials and Methods) in untreated conditions, where 98% and 83% of differential peaks significantly decreased after DAC or TDCIPP exposure, respectively (Fig. S2B-E). Additionally, we detected many loci in which H2A.Z decreases occurred for both exposure conditions, including 72% of sites affected by TDCIPP (Fig. 2C). Peaks where H2A.Z losses occurred tended to be genic, including promoters and intronic regions (Fig. S3A), whereas gains in H2A.Z occurred primarily over intergenic loci. H2A.Z levels were generally reduced across promoters (Fig. 2D), and promoters with the highest initial H2A.Z levels were affected most (Fig. S3B). Expression of genes proximal to these H2A.Z-marked promoters was also reduced (Fig. S3C), and metabolic processes were found to be associated with impacted genes (Fig. S3D,E). Peaks where H2A.Z increased (both for TDCIPP and DAC) overlapped more with REs than did peaks where H2A.Z decreased (Fig. 2E).

Despite these widespread changes in H2A.Z localization, nuclear H2A.Z protein levels remained unaltered (as measured by immunofluorescence; Fig. S4A,B). Additionally, differential gene expression analysis indicated that the observed outcomes were likely not due to a developmental delay. Genes that were upregulated following exposure to TDCIPP or DAC were found to be more highly expressed at later developmental stages than at earlier stages (Fig. S4C), and genes that normally undergo developmental activation during gastrulation were found to be upregulated following DAC and TDCIPP exposure at 6 hpf (Fig. S4D). Prior measurements of DNA methylation changes following DAC and TDCIPP exposure in cleavage stage zebrafish embryos reported only minor impacts (Volz et al., 2016). Our re-analysis of published datasets confirmed these findings (Fig. S4E), and we further verified (through anchor-based bisulfite sequencing measurements; Chapin et al., 2022) that DAC exposure during the cleavage phase of zebrafish development had only minor impacts on overall DNA methylation levels (Fig. S4F).

To assess the degree to which H2A.Z enrichment changes occurred over REs, we next performed differential enrichment analyses focusing specifically on RE subclasses. Unlike promoters, which exhibited widespread H2A.Z reduction, we found that H2A.Z levels primarily increased at REs upon exposure to either DAC or TDCIPP. DAC exposure led to significant H2A.Z increases at several repeat s­­ubclasses re­­presenting more than 600,000 uniquely mapped loci, and TCDIPP exposure cause similar outcomes, with increases representing more than 200,000 uniquely mapped REs (Fig. 3A). DAC and TDCIPP exposure caused 75% and 93% of differential repeat subclasses to gain H2A.Z (respectively), and by comparison, 99.5% and 98.5% of differentially impacted promoters lost H2A.Z (Fig. 3B,C). Additionally, a significant portion of changes (both at REs and promoters) were shared between exposure conditions (Fig. S4G), suggesting that they might represent a general feature of chemically induced RE reactivation.

Rather than locus-specific DNA methylation changes, we observed a modest overall reduction in methylation throughout the genome following DAC exposure, including REs that gained and lost H2A.Z (Fig. S4H). By contrast, decreases in the repressive histone modification H3K9me3 occurred preferentially over REs that gained H2A.Z, and H3K9me3 accumulated at loci that lost H2A.Z (Fig. 3D,E). As further support for H2A.Z functioning in transcriptional activation of REs, mean chromatin accessibility was initially higher at REs that gained H2A.Z as compared with those that lost H2A.Z, and DAC exposure led to a modest additional increase (Fig. 3F). We also observed an increase in RNA transcripts over REs that gained H2A.Z, and there were no such changes over REs that lost H2A.Z (Fig. 3G).

To investigate the widespread nature of genomic H2A.Z re-localization, we next assessed changes in median H2A.Z levels over mapped REs, and found that increases occurred for the majority of RE subclasses upon DAC exposure (93%) and TDCIPP exposure (62%) (Fig. 4A,B), and, interestingly, the total amount of H2A.Z that accumulated at REs was approximately equal to the amount lost at promoters (calculated as sequence-depth normalized total reads over each feature). H2A.Z levels increased from approximately 4.5×10⁶ to 6×10⁶ over impacted REs, and reads over promoters decreased by nearly the same amount, going from ∼5×10⁶ to 3.5×10⁶ (Fig. 4C). REs which acquired H2A.Z include tandemly repeated U1 (n=384) for the DAC exposure condition, Gypsy158-LTR_DR (n=202) following exposure to TDCIPP, and the DNA transposon DNA-8-10_DR (n=1723), which acquired H2A.Z upon TDCIPP or DAC exposure (Fig. 4D). Interestingly, impacted REs possessed an abundance of binding motifs for important developmental regulatory transcription factors, such as Nanog, KLF5 and TBX6 (Fig. S5A), and innate immune factors, such as STATs, p53 (TP53) and IRFs (Fig. S5B). Several studies have demonstrated that evolutionarily younger REs can sometimes function as co-opted enhancers and recruit transcription factors to control neighboring gene expression patterns (Chuong et al., 2013; 2016; Fueyo et al., 2022; Sundaram and Wysocka, 2020). Despite finding a reduction in overall branch length (measured via the nucleotide substitutions rate) (Fig. S5C) – indicating these elements are evolutionarily younger – we did not find evidence for expression changes at genes neighboring impacted REs (Fig. S5D), suggesting that they might not function as cis regulatory elements. Overall, we find that the phenotypic outcomes of exposure to TDCIPP and DAC coincided with genome-wide re-localization of H2A.Z, with H2A.Z accumulating ectopically at REs and decreasing at gene promoters.

In prior studies of zebrafish and mouse (Murphy et al., 2018; 2020), we found that disruption of H2A.Z removal machinery caused excess H2A.Z to accumulate across promoters, and, in the current study, we found that promoters generally lost H2A.Z. These observations led us to speculate that total H2A.Z levels might contribute to genomic localization patterns in early zebrafish embryos. To investigate this possibility, we took advantage of a transgenic zebrafish stain in which H2A.Z abundance was enlarged by overexpression of a fluorescently labeled transgene driven by the h2afx promoter (H2A.Z-mCherry; McMahon et al., 2009) (Fig. 5A,B). Indeed, we found that H2A.Z enrichment was higher over promoters in embryos overexpressing transgenic H2A.Z-mCherry, compared with for promoters in wild-type (WT) embryos (Fig. 5C,D), and there were modest increases over REs (Fig. 5E). Median H2A.Z enrichment was higher in mCherry embryos for 90% of RE subclasses, but these changes were modest as compared with DAC or TDCIPP treatment.

Significantly impacted REs represented fewer than 43,000 uniquely mapped loci, as compared with more than 600,000 and 200,000 loci for DAC and TDCIPP exposure, respectively (Fig. S6A). There was no significant overlap of impacted REs with those affected by DAC or TDCIPP (Fig. S6B). Interestingly, LTR elements were enriched among REs that gained H2A.Z (Fig. S6C,D), but accumulation of H2A.Z did not lead to RE activation (Fig. S6E), and we did not find evidence for interferon response gene activation in unexposed H2A.Z-mCherry-expressing embryos (Fig. S6F). Additionally, overexpression of H2A.Z-Cherry did not result in DNA methylation changes (Fig. S6G). These results indicate that increases in H2A.Z abundance correspond with increased H2A.Z installation both widely across REs and at promoter regions, and these broad impacts are largely distinct from the specific H2A.Z localization changes that occur following exposure to DAC or TDCIPP. Taken together, these results reveal that increased total amounts of H2A.Z in early zebrafish embryos correspond with increased enrichment of chromatin-bound H2A.Z, but these changes are insufficient to elicit an innate immune response.

Having found that H2A.Z localization patterns were influenced by total H2A.Z abundance (Fig. 5), and that the amount of H2A.Z lost at promoters upon exposure to DAC or TDCIPP (in WT embryos) was approximately equal to the amount of H2A.Z gained at REs (Fig. 4C), we next speculated that limits on H2A.Z availability might have contributed to phenotypic outcomes from TDCIPP and DAC exposure. Perhaps the burden to install H2A.Z at REs during TDCIPP and DAC exposure was too high for embryos to accommodate, exhausting limited H2A.Z pools and causing failed H2A.Z installation at promoters. To investigate this hypothesis, we exposed H2A.Z-mCherry embryos (which overexpress H2A.Z) to DAC and assessed phenotypic outcomes. As in WT embryos, DAC treatment triggered the expression of interferon response genes (Fig. S7A), but remarkably, embryos with increased total amounts of H2A.Z were partially resistant to DAC exposure (Fig. 6A). At 12 hpf, growth of the head and extension of the tailbud were dramatically stunted upon WT exposure to DAC, but overexpression of H2A.Z significantly rescued these defects (Fig. 6B,C). At 24 hpf, the majority (>75%) of WT embryos exposed to DAC displayed severe morphological defects, but in the presence of H2A.Z overexpression, severe defects were far less common (<25%), and many embryos appeared phenotypically normal or were only moderately impacted (Fig. 6D,E). Similarly, H2A.Z-mCherry embryos also showed less severe morphological defects than WT embryos in response to TDCIPP treatment (Fig. 6E).

If H2A.Z abundance is limited in WT embryos, then ectopic accumulation of H2A.Z at REs during DAC or TDCIPP treatment might have occurred at the expense of promoters. To assess whether such limited conditions exist, we further profiled H2A.Z enrichment in mCherry-expressing embryos exposed to DAC (Fig. S7B), anticipating that H2A.Z availability might no longer be limited. As anticipated, reduction of H2A.Z at promoters was quite modest following DAC exposure in mCherry-expressing embryos, as compared with WT (Fig. 7A,B; Fig. S7C), and DAC exposure led to modest additional increases over REs, which occurred primarily at sites where H2A.Z levels were highest (Fig. 7C). Compared with WT embryos, the magnitude of H2A.Z changes at REs in H2A.Z-mCherry embryos was considerably reduced (Fig. 7D; Fig. S7D), and many affected REs already exhibited increased H2A.Z levels in unexposed mCherry embryos, compared with WT, potentially explaining the more modest additional increases that occurred following exposure to DAC (Fig. 7E). In support of these findings, principal component analysis (PCA) of H2A.Z enrichment over RE families indicated that untreated H2A.Z-mCherry embryos were relatively similar to DAC-treated WT embryos (along PC1, representing 86.8% variance), whereas analysis of promoters indicated both treated and untreated H2A.Z-mCherry embryos were more similar to untreated WT embryos (Fig. S7E).

To investigate the degree to which DAC-induced gene expression changes were affected by H2A.Z-mCherry expression, we next performed whole-embryo transcriptome analysis at 12 hpf (Fig. S7F,G). Once again, impacts of gene expression were consistent with mitigated phenotypic outcomes (Fig. 7F-H). Transcriptional changes for both up- and downregulated genes were partially reversed in the presence of H2A.Z-mCherry (Fig. 7F,G). Gene ontology analysis revealed that downregulated genes tended to function in developmental progression, whereas upregulated genes tended to function in stress response, apoptosis and cell death (Fig. 7H). Additionally, REs that were activated following DAC treatment in WT were also activated upon DAC treatment in H2A.Z-mCherry embryos, but to a lesser extent (Fig. S7H). In sum, these results demonstrate that total H2A.Z abundance establishes sensitivity to repetitive element de-repression (via DAC and TDCIPP exposure) in developing embryos, and changes in H2A.Z localization are largely dependent on the amount of H2A.Z that is present during exposure in embryos (Fig. 7I).

This study demonstrates that differences in H2A.Z protein abundance significantly influence how embryos are impacted by external stimuli that trigger RE activation. Low H2A.Z levels confer sensitivity, whereas high levels confer robustness or resistance, implying that individual organisms or tissues with varying levels of H2A.Z might be more or less sensitive to external stressors. Interestingly, prior studies have found H2A.Z to be overexpressed in many forms of cancer (Kim et al., 2013), perhaps providing tumors with a means to evade cell death following innate immune activation. In this regard, targeting H2A.Z, or regulators of H2A.Z installation, might uncover effective strategies to combat cancer, especially in the context of retroviral activation, or in tumors that are resistant to viral mimicry-based therapies (Deblois et al., 2020).

Numerous cell culture-based studies have demonstrated that DAC exposure leads to genome-wide reduction in DNA methylation (Klco et al., 2013; Navada et al., 2014; Yang et al., 2012), but we did not observe such extensive impacts following exposure in early zebrafish embryos. Rather, we found that DAC exposure caused H2A.Z accumulation at REs, which was accompanied by loss of H3K9me3 and modest increases in chromatin accessibility. It remains unknown how such differences from in vivo versus in vitro cell-culture-based DAC exposure arose. Notably, several prior studies of DAC exposure, used as a clinical treatment for acute myeloid leukemia and myelodysplastic syndrome, have arrived at conclusions similar to ours, including studies in which patients responded positively to DAC treatment independently of detectible DNA methylation changes (Komashko and Farnham, 2010; Schmelz et al., 2005; Yang et al., 2006). Additionally, several prior studies have proposed that DAC treatment is effective largely through initiation of DNA damage response mechanisms (Karpf et al., 2001; Palii et al., 2008; Widodo et al., 2007) rather than epigenetic de-repression. As in our study, it has been previously reported that histone modification changes occur following DAC exposure without evidence for significant DNA methylation changes (Komashko and Farnham, 2010; Palii et al., 2008; Takebayashi et al., 2001). Evidently, mechanisms by which DAC exposure leads to RE activation and the innate immune response remain elusive, and although we find it unlikely, it remains a possibility that undetectable locus-specific DNA methylation changes did occur in our study, and these limited changes led to the observed outcomes.

RE types impacted by both DAC and TDCIPP were classified as DNA transposons, LTRs and LINEs. Interestingly, many attributes of these zebrafish elements are shared with REs that function as enhancers during the mammalian innate immune response, including enrichment for immune response transcription factor motifs (e.g. STATs, IRFs and p53), and low nucleotide substitution frequencies, suggesting that they might be subject to purifying selection, or that these REs arose in recent evolutionary history, as compared with non-impacted elements. Unlike their mammalian counterparts, immune-associated REs in zebrafish occurred across several distinct types of repeats, including many not present in mammals, and we did not detect expression changes for genes neighboring impacted REs, thus indicating that co-option of these REs in zebrafish likely occurred separately from mammals, potentially due to shared evolutionary pressures. Another possibility is that only a subset of REs function as enhancers – a prospect that warrants further investigation.

In addition to transcription factor motifs involved in the interferon response, impacted REs also possessed an abundance of motifs for transcription factors known to control transcription during early embryonic development, including Myc, Klf and Nanog. These results are reminiscent of studies in mammals, in which activation of enhancer-like endogenous retroviruses promotes embryonic genome activation and pluripotency (Hendrickson et al., 2017; Ishiuchi et al., 2015; Macfarlan et al., 2012), or of neural crest development in zebrafish, in which the neural crest marker gene crestin is derived from an ancestral RE (Rubinstein et al., 2000). Prior studies in zebrafish indicate that similar REs might even function in the establishment of higher-order chromatin structure during embryonic development (Meier et al., 2018). Interestingly, recent studies of mouse embryonic stem cells have found that de-repression of endogenous retroviral elements enables REs to ‘hijack’ RNA polymerase II (Asimi et al., 2022) or the transcription factors KLF5 and KLF9 (O'Hara and Banaszynski, 2022 preprint). Based on the results of the current study, we suspect H2A.Z might be similarly hijacked. Evidently, further research is necessary to fully establish the role of REs during vertebrate developmental transitions and/or the emergence of particular lineages.

We propose a source–sink model to describe the role of REs in buffering against environmental insults during early zebrafish development. This model reflects our results in which both strategies used for RE de-repression, DAC or TDCIPP exposure, lead to reduced H2A.Z enrichment over promoters and increased enrichment at REs, and that augmenting total H2A.Z levels coincides with mitigated phenotypic outcomes. When H2A.Z levels are low, in WT embryos, we propose that REs function as a sink, sequestering H2A.Z away from promoters during exposure. In embryos overexpressing H2A.Z-mCherry, however, REs act as a storage location for excess H2A.Z, protecting the embryos from the consequences of subsequent exposure. It is also conceivable that the high levels of H2A.Z we observed over REs in H2A.Z-mCherry embryos were the result of increased stress or inflammation specific to these transgenic animals. If true, the observed resistance to DAC and TDCIPP could have been due to an epigenetic ‘memory’, potentially allowing for resistance to future stressors. In this context, the phenotypic rescue we observed could have been partial due to less severe H2A.Z localization changes both at gene promoters and at REs upon DAC treatment in H2A.Z-mCherry embryos. Detectable changes in H2A.Z localization following exposure could have occurred by altering H2A.Z targeting mechanisms, in a manner preferential to REs, or by H2A.Z protein translocations from promoters to activated REs. In the current study, we are unable to distinguish between these possibilities. We also note how H2A.Z accumulated at a distinct set of REs following DAC exposure as compared with TDCIPP, despite considerable overlap. Mechanisms underlying this specificity remain unknown and are an interesting topic for future pursuit.

Interestingly, prior studies have found that H2A.Z enrichment decreases over promoters following immune activation in mammalian cells (Au-Yeung and Horvath, 2018), but the function of H2A.Z in this context and impacts at REs have not been addressed. As several studies have demonstrated, REs often become active during immune response and/or viral infection (Chuong et al., 2016; Macchietto et al., 2020; Roy et al., 2020; Shen et al., 2022). If H2A.Z accumulation occurs over transcriptionally active REs, as we have found in the current study, then this raises an interesting question pertaining to the evolutionary benefits of limiting H2A.Z abundance. One compelling possibility is that limits on H2A.Z availability might help to prevent the spread of pathogens. In this scenario, if a viral infection were to occur, shuttling H2A.Z toward REs and away from developmental gene promoters would inhibit embryogenesis and lead to pragmatic lethality. Maintaining low total protein levels for H2A.Z would thus provide organisms with a convenient cell death induction escape route for combating the spread of infectious pathogens. It is therefore logical that epigenetic wiring of the nucleus would occur in a manner where REs and gene promoters compete with one another for H2A.Z (Fig. 7I).

Like zebrafish, the genome of humans largely comprises REs (more than 45%) (Lander et al., 2001), and newer sequencing technologies are likely to uncover additional functions for newly mapped elements (Hoyt et al., 2022). In this context, the mechanism we present here, whereby H2A.Z is sequestered to particular REs upon their activation, might represent a general principle in eukaryotic biology that is independent of the direct influence of REs themselves. Although we focus specifically on H2A.Z, it is certainly possible that other nuclear factors are sequestered in a similar manner. Thus, our study presents a conceptual advance in the broader understanding of transcription control, whereby repetitive DNA sequences and cell-type-specific genes networks engage in an intrinsic competition for transcription machinery, in effect governing overall cellular output while simultaneously enabling response to external stimuli.

Sources for antibodies, fish strains and chemicals, and software and algorithms, are given in Table S1.

Our re-analysis of previously published data sets used the following resources: a DAC-treated RNA-Seq dataset (Potok et al., 2013; SRA: SRP020008); a DNA methylation dataset (Potok et al., 2013; SRA: SRP020008); a TDCIPP-treated RNA-Seq dataset (Dasgupta et al., 2018; BioProject ID: PRJNA475635); a DNA methylation dataset (Volz et al., 2016; BioProject ID: PRJNA330715); and H3K9me3 ChIP-seq (Laue et al., 2019; GEO: GSE113086) and H3K27ac ChIP-seq (Zhang et al., 2018; GEO: GSE114954) datasets.

Adult WT Tübingen and H2A.Z-mCherry zebrafish were maintained on a 14 h light:10 h dark cycle. Zebrafish husbandry and care were conducted in full accordance with animal care and use guidelines with ethical approval by the University Committee on Animal Resources at the University of Rochester Medical Center. Experimental samples were early zebrafish embryos ranged from 0 hpf to 24 hpf.

Newly fertilized zebrafish eggs were collected after spawning, and ∼50 embryos per well were maintained in a six-well plate within a light- and temperature-controlled incubator. As in prior studies, embryos were treated in triplicates with 5 ml of 3.12 µM TDCIPP (MedChemExpress, HY-108712, treatment started from the two-cell stage at 0.75 hpf), or 100 µM of DAC (Cayman Chemical, 11166, treatment started from the one-cell stage), with 0.1% DMSO used for control group, and incubated until 6 hpf for H2A.Z CUT&Tag and 12 hpf for RNA-seq, or 6 hpf, 12 hpf and 24 hpf for morphological characterizations.

RNA-seq reads mapping to coding sequences of mCherry, h2afva (h2az2a) and actb2, were aligned to a custom index using bowtie2 and average read counts were measured using deeptools. Read coverage was then assessed over the mCherry coding sequence and the h2afva C-terminus, which represents the portion of H2A.Z where the mCherry fusion occurred. The sequence over this portion of the h2afva cDNA is identical to sequence from h2afvb (h2az2b). Therefore, reads were distributed equally across both h2az alleles. Read coverage was then compared with coverage over actb2 (encoding β-actin), to account for differences in sequencing depth between samples. These normalized read coverage values were then used to calculate relative mRNA expression for the h2afva-mCherry transgene.

Embryos were collected at 6 hpf and fixed with 4% paraformaldehyde in 1× PBS overnight. Embryos were washed briefly several times in PBX (1× PBS containing 0.1% Triton X-100), then blocked overnight with PBX containing 0.5% bovine serum albumin, then incubated with anti-H2A.Z antibody overnight (1:250 in PBX containing 0.5% bovine serum albumin), then Cy5 anti-rabbit-IgG secondary antibody was used at 1:500. To quantify H2A.Z intensity, confocal images were processed with Fiji to measure H2A.Z level within individual nuclei.

To quantify phenotypic outcomes of DAC and TDCIPP exposure, images of treated zebrafish embryos were captured at 6 hpf, 12 hpf or 24 hpf using the levenhuklite software package with a bench-top stereo microscope. To compare phenotypes of DAC-treated WT and H2A.Z-mCherry embryos at 12 hpf, lateral images were taken for each embryo at the same setting, and a straight line was drawn and measured in ImageJ (version 2.0.1/1.53c) using the ‘Straight’ tool between the head and tail, and between edge of head and yolk sac for maximum head width. Violin and boxplots were used to present the measures for WT and H2A.Z-mCherry embryos treated with DMSO and DAC. To compare phenotypes of DAC-treated WT and H2A.Z-mCherry embryos at 24 hpf, embryos were scored into four groups based on severity of morphological defects – normal: embryos had a well-developed head and normal length of tail; mild: embryos with well-developed head and slightly shortened length of tail; moderate, embryos with poorly developed head and obviously shortened tail; severe, embryos with barely developed head and dramatically shortened tail. The scoring was carried out prior to genotyping of transgenic animals by a researcher who was aware of the experimental conditions.

Genomic DNA was extracted using phenol-chloroform from 6 hpf developing zebrafish embryos after either DMSO or DAC treatment. Extracted genomic DNA was used to generate sequencing library based on the anchor-based bisulfite sequencing protocol (Chapin et al., 2022) and was sequenced using paired-end reads. Raw reads were aligned to zebrafish genome build Zv11 using Bismark (v0.23.1) according to the published method (Chapin et al., 2022). Specifically, bismark -p 4 --parallel 6 -q --pbat -o −1 reads1.fastq.gz −2 reads2.fastq.gz and PCR duplicates were removed using deduplicate_bismark. Mapped reads were further filtered based on MAPQ score 20 using samtools view -h -q 20, then methylated cytosines were extracted using bismark_methylation_extractor --parallel 10 -p --no_overlap --gzip --bedGraph --zero_based --cytosine_report. DNA methylation genomic tracks (CpG sites with coverage <10 were filtered out) were generated using bedgraphtobigwig.

Treated embryos (∼25 embryos per replicate) were collected at 6 hpf and manually dechorionated. Nuclei were isolated and processed according to previously described protocol to generate ATAC-Seq libraries. After tagmentation using Tn5, genomic DNA were extracted using Zymo Research DNA Clean & Concentrator-5 (D4014), and ATAC-Seq libraries were generated using NEBNext High-Fidelity 2× PCR Master Mix (M-541S).

Treated embryos were collected at 6 hpf and manually dechorionated. Then, embryos were transferred into 1.5 ml Eppendorf tubes, and dissociated by pipetting. Cells were pelleted after centrifugation at 600 g for 3 min, and nuclei were extracted on ice for 10 min, and collected by centrifugation at 600 g for 3 min. Around 100,000 nuclei were used for each replicate, and CUT&Tag libraries were generated according to published CUT&Tag protocols (Kaya-Okur et al., 2019). Polyclonal anti-H2A.Z antibody (Active Motif, 39113, 1:50) and polyclonal anti-H3K9me3 antibody (Active Motif, 39062, 1:50) were used for CUT&Tag.

CUT&Tag and ATAC-Seq libraries were pooled and sequenced using the Illumina paired-end sequencing HiSeq 4000 platform. Raw sequencing reads were checked for quality using FastQC (version 0.11.9), and adapter sequences were removed using Cutadapt (version 2.7). Sequencing reads were aligned to zebrafish genome assembly (GRCz11) using Bowtie2 (version 2.2.5) with default setting. Sequence duplicates were removed using Picard MarkDuplicates (version 2.5.0), and genome browser bigwig files for H2A.Z were generated using merged replicate bam files and converted to bigwig files using bamCoverage (version 3.5.1, --normalizeUsing RPKM --binSize 10). The Integrative Genomics Viewer (version 2.8.13) was used to visualize all sequencing data and to prepare genome browser snapshots.

To generate quantile normalized bigwig files, H2A.Z CUT&Tag scores across genomic bins of 100 bp were generated using multibigwigsummary bins --binSize 100. Then the matrix of score was quantile normalized using the normalize.quantiles function in the R package preprocessCore. To generate quantile normalized genomic tracks, new bedgraphs were generated using quantile normalized H2A.Z scores and converted into bigwig files using bedgraphtobigwig.

To identify H2A.Z occupied genomic regions, peak calling of H2A.Z CUT&Tag was performed using macs2 callpeak (version 2.2.6) with the following parameters: -f BAMPE -g 1.5e9 --nomodel --broad. Genomic annotation of H2A.Z peaks was performed using homer annotatePeaks.pl to assess annotations, distance to nearest promoters and CpG content.

To identify differential H2A.Z peaks between DMSO- and DAC- or TDCIPP-treated embryos, H2A.Z broadPeak and associated bam files for each replicate were used as input files for DiffBind (version 3.2.5) with default settings. Significantly increased H2A.Z peaks were defined as having an adjusted P<0.05 and log2 fold change>0.5, whereas significantly decreased H2A.Z peaks were defined as having an adjusted P<0.05 and log2 fold change<−0.5.

To plot H2A.Z levels at REs as heatmaps, reads per kilobase per million mapped reads (RPKM)-normalized H2A.Z bigwig files of 6 hpf treated WT and H2A.Z-mCherry embryos were used to generate coverage scores at each genomic copy of each uniquely named repeat (simple repeats were excluded from the analysis). Then median H2A.Z coverage scores were calculated for each repeat subclass and used as the input file in R to generate H2A.Z heatmaps. Heatmaps were grouped based on their repeat superfamilies and ranked based on total H2A.Z scores of DMSO and DAC in 6 hpf WT embryos.

DMSO- and DAC-treated WT and mCherry-positive embryos at 12 hpf were used to extract RNA. Chorions were removed with forceps, and embryos were collected in TRIzol and grinded with a mortar and pestle. Total RNA was extracted with phenol-chloroform-isoamyl alcohol, then the aqueous layers were purified with a Zymo Direct-zol RNA Miniprep kit. Total RNA-seq libraries were prepared and sequenced with paired-end 150 bp reads by LC Sciences. Sequencing reads were mapped with STAR with the default settings (v2.7.2a), and mapped reads were used in featureCounts (v2.0.3 in subread) using GRCz11.102.gtf.

To examine changes in gene expression after DAC or TDCIPP exposure at 6 hpf, featureCounts (v2.0.3 in subread) was used to calculate reads over exons. Read counts were normalized to sum of read counts of each sample, and the log2 fold change was calculated in R software. Differentially expressed genes were identified based on log2 fold change: log2 fold change>0.5 for increased genes, and log2 fold change<−0.5 for decreased genes. To examine temporal gene expression patterns during developmental progression, a read counts normalized (TPM) gene expression table from the one-cell stage to day 5 was used (White et al., 2017). Particular gene sets were merged with count table in R software to calculate average expression profile for boxplots.

To identify partially rescued genes in DAC-treated H2A.Z-mCherry embryos at 12 hpf, differentially expressed genes in DAC-treated WT embryos were identified using DESeq2 (log2FC at 0.5 and adjusted P<0.05). Among those genes, partially rescued increased genes were those that had less than an 0.5 difference in log2FC in H2A.Z-mCherry embryos compared with WT embryos (log2FC_H2A.Z-mCherry−log2FC_wildtype<0.5). Similarly, partially rescued decreased genes were those more than 0.5 difference in log2FC in H2A.Z-mCherry embryos (log2FC_H2A.Z-mCherry−log2FC_wildtype>0.5).

To quantify transcription levels of REs, featureCounts (v2.0.3 in subread) was used to generate the raw count table for each repeat subclass from RNA-seq data. Read counts were normalized to sum of read counts of each sample, and log2 fold change was calculated in R and simple repeats were excluded for downstream analysis.

Differential H2A.Z peaks and gene groups were used as input files, and enriched gene ontology terms associated with these regions and genes were identified by g:Profiler. The significance threshold was set at 0.05 using g:SCS threshold for multiple testing.

To identify promoters with differential H2A.Z levels, bam files of H2A.Z CUT&Tag and bed files of promoters (TSS±1 kb) were used as input files to run DiffBind, based on an adjusted P<0.05 and absolute value of log2FC at 0.5. Unlike promoters, owing to the repetitive nature of REs and mapping reads assignment, it is difficult to assess H2A.Z changes at an individual RE locus. Instead, we performed an aggregate analysis and counted reads mapped to all loci belonging to the same RE classes. To identify REs with differential H2A.Z and H3K9me3 levels, featureCounts (v2.0.3 in subread) was used to generate a raw count table for each repeat subclass using paired-end H2A.Z and H3K9me3 CUT&Tag reads from 6 hpf DMSO-, TDCIPP- and DAC-treated embryos. Then count table was used in DESeq2 to identify differentially marked RE subclasses (adjusted P<0.05 and absolute value of log2FC at 0.5) and simple repeats were excluded for downstream analysis.

To quantify the percentage of REs, repeats subclasses within the same family were counted, including DNA transposons, LTRs, LINEs, short interspersed nuclear elements, tRNA, satellite DNAs and simple repeats. Pie charts were generated using R (version 4.0.3).

To identify enriched transcription factor-binding motifs in repetitive elements gaining H2A.Z upon DAC and TDCIPP exposure, bed files for repetitive elements with increased H2A.Z were used as input to HOMER findMotifsGenome.pl. ‘-size given’, and bed files for REs with decreased H2A.Z were used as background regions for motif discovery.

To perform PCA on H2A.Z CUT&Tag replicates at RE families, RPKM-normalized H2A.Z bedGraph files were used as input files to calculate counts for RE families (simple repeats were filtered out for this analysis). Homer analyzeRepeats.pl was applied to RepeatMasker GTF file, and output counts without additional normalization were used directly for PCA plotting, using R packages prcomp and ggfortify. To perform PCA on H2A.Z CUT&Tag replicates at promoter regions, multibigwigsummary output count table using RPKM normalized bigwig files was used for PCA. To perform PCA on RNA-seq replicates at 12 hpf, DEseq2 plotPCA was used.

Heatmaps and profile plots were generated using plotHeatmap and plotProfile from deepTools (version 3.5.1) or pheatmap in R, and fonts and labels were adjusted in Affinity Designer (version 1.9.1). Density scatterplots were generated using ggplot2 geom_density. All boxplots were plotted using R; the horizontal line shows the median, and the box encompasses the interquartile range. The whiskers show the 75th percentile plus 1.5× the interquartile range (upper whisker) and the 25th percentile minus 1.5× the interquartile range (lower whisker). Statistical tests used can be found in the figure legends.

We thank Drs Amanda Larracuente and Paula Vertino, as well as everyone in the Murphy lab for helpful discussions.