Loss of DNA methylation in zebrafish embryos activates retrotransposons to trigger antiviral signaling

Cytosine methylation promotes heterochromatin formation on terminally repressed regions, such as repeats (Biscotti et al., 2015; Yoder et al., 1997), imprinted genes (Li et al., 1993) and transposable elements (TEs) (Yoder et al., 1997). These canonical roles for DNA methylation are well established; however, little is known about the physiological responses to epigenetic stress when these functions fail. DNA hypomethylation is associated with diseases such as cancer and autoimmune disorders, making this topic both important and relevant. We used zebrafish embryos with mutations in the central DNA methylation machinery to identify the cellular and physiological responses to global DNA hypomethylation.

Despite intense work in the DNA methylation field, its role in regulating differentially expressed genes is controversial (Gutierrez-Arcelus et al., 2013). Many studies in organisms across the evolutionary spectrum have demonstrated that DNA methylation does not serve as a universal mechanism of repressing gene expression (Feng et al., 2010; Grow et al., 2015; Gutierrez-Arcelus et al., 2013; Jackson-Grusby et al., 2001; Jacob et al., 2015; Potok et al., 2013; Qi et al., 2015; Zemach et al., 2010; Zhang et al., 2016). In fact, most data indicate that relatively few genes are induced when methylation is lost: there is a lack of global gene induction in mammalian cells deficient in the proteins required for DNA methylation, namely ubiquitin-like with PHD and ring finger domains 1 (Uhrf1) (Qi et al., 2015) or DNA methyltransferase 1 (DNMT1) (Jackson-Grusby et al., 2001), or during the stage of early embryogenesis when the parental methylome is nearly completely erased (Grow et al., 2015; Potok et al., 2013). Similarly, we found no evidence of widespread genome activation in zebrafish uhrf1 mutants (Jacob et al., 2015). Instead, we found that genes upregulated in uhrf1 mutants fell essentially into two major categories: cell cycle and immunity (Jacob et al., 2015). Our previous work examined the functional relevance of cell cycle gene activation (Jacob et al., 2015). Here, we investigate the basis for immune gene induction in zebrafish mutants that lack DNA methylation due to uhrf1 (Sadler et al., 2007) or dnmt1 (Anderson et al., 2009) mutations.

Although DNA methylation contributes to the regulation of some genes, many studies have led to the conclusion that this function is the exception rather than the rule. Indeed, we find no evidence that the genes induced in uhrf1 mutants are directly regulated by DNA methylation (Jacob et al., 2015). What, then, is the primary purpose of DNA methylation? In most organisms, transposons are the most heavily methylated regions of the genome. Several recent examples show widespread induction of retrotransposon expression in cancer cells (Chiappinelli et al., 2015; Leonova et al., 2013; Roulois et al., 2015), embryonic stem cells (Sharif et al., 2016), neural stem cells (Ramesh et al., 2016) and in human embryos (Grow et al., 2015) when the genome becomes hypomethylated, indicating that transposon repression is the central repressive role of DNA methylation. Interestingly, in several of these cases, the same set of immune genes is upregulated in uhrf1 mutants. Thus, transcriptome data highlight the conserved function of DNA methylation to repress transposons and prevent activation of the immune system.

The relationship between TE induction and immune activation has recently been investigated. Prokaryotic and viral genomes are not methylated, and this ʻnon-self' signal triggers the host immune system to detect and clear infected cells (Singer et al., 2015; Zhong et al., 2006). Unmethylated DNA in the cytoplasm elicits a nucleotide detection signaling pathway involving DAI (ZBP1), cGAS (MB21D1), DDX41 and AIM2 and a complementary RNA-sensing pathway that involves RIG-1/MDA5 (IFIH1) (Dempsey and Bowie, 2015). The DNA- and RNA-sensing pathways activate stimulator of interferon genes (STING, or TMEM173) and mitochondrial antiviral signaling protein (MAVS), respectively, and these both converge on TANK-binding kinase 1 (TBK1) to translate the non-self signal into interferon production and activation of a systemic antiviral response (Dempsey and Bowie, 2015). The model that emerges is that loss of DNA methylation leads to derepression of endogenous retrotransposable elements, mimicking a viral infection. The response activated as a result can serve to clear these pseudo-infected cells. Thus, TE derepression might be a mechanism to flag cells with epigenetic stress (Milutinovic et al., 2003; Timp and Feinberg, 2013) for immune surveillance.

We used uhrf1 and dnmt1 mutant zebrafish with global DNA hypomethylation (Anderson et al., 2009; Jacob et al., 2015; Sadler et al., 2007) to determine how loss of genomic methylation leads to induction of the innate immune system in the developing embryo. We identified a robust activation of type I interferon, the first line of action in antiviral signaling, and an expanded population of immune cells, which could not be explained by loss of DNA methylation in the promoter of the upregulated genes. Instead, we found that Sting and phosphorylated Tbk1 (pTbk1) were induced in uhrf1 and dnmt1 mutants accompanied by widespread induction of retrotransposons. Blocking either Sting, pTbk1 or production of cytoplasmic DNA by inhibiting retrotranscription repressed the interferon response in uhrf1 mutants. We conclude that TE repression is a primary function of DNA methylation during vertebrate development and suggest that antiviral pathways serve to mark epigenetically damaged cells for clearance by the immune system.

Uhrf1 recognizes hemi-methylated DNA (Arita et al., 2008; Avvakumov et al., 2008; Hashimoto et al., 2008; Qian et al., 2008) and recruits Dnmt1 during DNA replication (Bostick et al., 2007; Sharif et al., 2007). Thus, total cytosine methylation [5-methylcytosine (5MeC)] levels in uhrf1 mutant embryos are reduced to less than half of wild-type (WT) levels (Feng et al., 2010; Jacob et al., 2015; Tittle et al., 2011). Transcriptome analysis of uhrf1 mutants at 120 h post fertilization (hpf), a time point when the morphological phenotype of these mutants is fully evident, was carried out using both microarray and RNA-seq. These approaches revealed that genes annotated as having a function in the immune system dominated the category of upregulated genes (Fig. 1A, Tables S1 and S2), and that this pattern was mirrored in dnmt1 mutants analyzed by RNA-seq at 120 hpf (Fig. 1A, Table S1). Gene set enrichment analysis (GSEA) using the Gene Ontology (GO) database (geneontology.org) showed that immune-related genes and apoptosis are two of the most robust and significantly enriched pathways in uhrf1 mutants (Fig. 1B, Table S2). We reasoned that the apoptosis genes are largely related to the high degree of cell death observed in uhrf1 mutants (Jacob et al., 2015; Sadler et al., 2007; Tittle et al., 2011). Nearly all the immune genes identified from the microarray were confirmed as significantly upregulated in uhrf1 (Fig. 1C, Fig. S1, Table S3) and dnmt1 (Fig. S1B) mutants at 120 hpf. On the other hand, genes significantly downregulated in uhrf1 mutants correlated with reduced organs and tissues that accompany the mutant phenotype (Table S4).

Interferon signaling is the first line of defense to viral infection, and all of the interferons, their receptors and other key players in the interferon response are well conserved across vertebrates (Briolat et al., 2014; Hamming et al., 2011; Robertsen, 2006; Zou and Secombes, 2011). Genes dysregulated in uhrf1 mutants were highly correlated with those classified as involved in the type I and type II interferon response in human cells (Fig. 2A, Table S2). Many of the upregulated genes in uhrf1 mutants were also enriched in gene sets from human cells treated with interferons (Fig. 2B, Table S2). Zebrafish infected with the chikungunya virus (CHIKV) mount a type I interferon response (Briolat et al., 2014; Palha et al., 2013), whereas those infected with the hematopoietic necrosis virus (IHNV) do not (Briolat et al., 2014). We compared the previously published gene expression pattern from zebrafish infected with these two viruses (Briolat et al., 2014) with the gene expression in uhrf1 mutants and found that the pattern of expression of the top ten significantly upregulated and downregulated genes from CHIKV-infected embryos was similar to that of uhrf1 and dnmt1 mutants (Fig. 2C). However, there was not a strong correlation with the gene expression pattern induced by the IHNV, which does not elicit a robust interferon response (Briolat et al., 2014) (Fig. 2C). A similar pattern was observed in RNA-seq analysis of 120 hpf dnmt1 mutants (Fig. 2C). Furthermore, ifnphi1, the zebrafish ortholog of mammalian type I interferon (Aggad et al., 2009), was significantly upregulated in uhrf1 mutants at 120 hpf (Fig. 2D). Therefore, although the embryos in our study were not infected, the gene expression pattern indicates that uhrf1 and dnmt1 mutants mount an interferon response as a result of loss of genome-wide DNA methylation similar to that observed upon viral infection.

We next investigated the timecourse of expression of select interferon-response genes and ifnphi1 by RT-qPCR (Fig. 3A) in parallel with assessment of global DNA methylation status using slot blotting of genomic DNA for 5MeC (Fig. 3B). We find that some genes (ifnphi1, nfkb2 and irf1b) are detected as upregulated by RT-qPCR as early as 55 hpf (Fig. 3A) and became further overexpressed in older uhrf1 mutant larvae (Fig. 3A, Fig. 1C, Fig. S1, Table S3). These genes were also detected as significantly upregulated by RNA-seq analysis of uhrf1 and dnmt1 mutants at 120 hpf (Fig. 3A, Table S1). Although RNA-seq analysis of 55 hpf uhrf1 mutant embryos failed to detect robust gene expression changes (Fig. 3A, Table S1), this is likely to reflect differences in assay sensitivity between RT-qPCR and RNA-seq (Seqc/Maqc-Iii Consortium, 2014). Regardless, the finding that some interferon-response genes are upregulated as early as 55 hpf (by RT-qPCR) is significant, as this immediately follows the global depletion of uhrf1 maternal mRNA (at 48 hpf) and DNA hypomethylation (Fig. 3B) (Jacob et al., 2015). Importantly, these early responses precede any detectable phenotypic features of uhrf1 mutants (Jacob et al., 2015). Taken together, these data show that induction of type I interferon expression is the earliest detectable change in uhrf1 mutants following DNA hypomethylation.

Interferons are produced by a variety of cell types to recruit and activate leukocytes (Le Page et al., 2000; Palha et al., 2013). Innate immunity is functional at the developmental stages studied here, whereas adaptive immunity is not yet fully developed (Masud et al., 2017), and previous studies found that neutrophils and hepatocytes account for most ifnphi1 production in virally infected zebrafish (Palha et al., 2013). We asked whether an expanded population of leukocytes caused the induction of interferon genes in uhrf1 mutants. Transgenic markers of macrophages [Tg(mpeg:mCherry) (Ellett et al., 2011)] and neutrophils [Tg(lysZ:dsRed) (Hall et al., 2007)] revealed expanded populations of these cell types in uhrf1 mutants as early as 80 hpf, and this further increased by 120 hpf (Fig. 4A-D) when leukocytes were spread throughout the larvae (Fig. S2). These results were confirmed using Neutral Red to stain for macrophages (Herbomel et al., 2001), showing that the head and jaw of uhrf1 (Fig. S3A,B, arrows) and dnmt1 (Fig. S3C) mutants had significantly increased Neutral Red staining. Together, these results illustrate that the population of macrophages is expanded in both uhrf1 and dnmt1 mutants.

Hypomethylated DNA from pathogens can be an immunostimulant (Hemmi et al., 2000; Yeh et al., 2013), and we asked if extracellular hypomethylated DNA could contribute to the expanded leukocyte population in uhrf1 mutants. We injected the otic vesicle of 72 hpf Tg(mpeg:mCherry) embryos with sheared genomic DNA isolated from 120 hpf uhrf1 mutants or WT siblings. Both sources of DNA increased the number of macrophages recruited to the injection site (Fig. S4A), and we reasoned that genomic DNA from whole embryos contains cells with varying amounts of DNA methylation and thus the genomic DNA from uhrf1 mutants was not fully unmethylated. To test this, Tg(mpeg:mCherry) embryos were injected with either completely unmethylated or completely (CpG) methylated oligodeoxynucleotide. The unmethylated oligo recruited a significantly greater number of macrophages (Fig. S4B). This shows that unmethylated DNA is a strong immunostimulant, as viral and microbial DNA is entirely unmethylated at CpGs and this feature of DNA can be used by the immune system to distinguish foreign nucleic acid from self. We hypothesize that dying cells in uhrf1 mutants release DNA into the extracellular space, stimulating recruitment of macrophages and neutrophils. Another intriguing possibility is that uhrf1 loss in innate immune cells promotes their activation, as has been recently suggested in humans (Yao et al., 2016).

We tested whether the expanded population of immune cells was the source of interferon and immune gene expression in uhrf1 mutants by blocking all leukocyte development using a morpholino targeting pu.1 (spi1b) (Rhodes et al., 2005). This achieved marked reductions in mpeg:mCherry-positive and lysZ:dsRed-positive cells in both WT and uhrf1 mutant embryos at 80 hpf (Fig. 4E), but surprisingly had no effect on the morphological phenotype of uhrf1 mutants (Fig. 4E) or the induction of the interferon gene panel (Fig. 4F, Fig. S5). Similarly, depleting macrophages using a splice-blocking morpholino that effectively targeted irf8 (Fig. S6A-C) had no effect on the morphological phenotype of uhrf1 mutants (Fig. S6C) or on the expression of interferon genes (Fig. S6D,E). These data demonstrate that neutrophils and macrophages are not required for robust interferon gene expression, nor for the gross morphological defects observed in uhrf1 mutants.

Next, we determined which tissues respond to interferons in uhrf1 mutants, using in situ hybridization to detect interferon target gene expression. irf1b and isg15 were not detected in WT embryos but were highly expressed in the jaw, head, eye, liver and gut of uhrf1 mutants at 80 and 120 hpf (Fig. 5). This differs from the pattern of leukocytes in uhrf1 mutants (Fig. 4A,B), further supporting the hypothesis that immune cells are not the primary source of the interferon signaling in uhrf1 mutants.

To ask directly whether loss of Uhrf1 could induce immune gene expression in the absence of immune cells, we used siRNA to knock down UHRF1 in cultured human hepatoma cells (HuH7). Both total and phosphorylated STAT1 protein (Fig. S7A) and the expression of interferon-response genes were increased in UHRF1-depleted cells (Fig. S7B). Taken together, these data demonstrate that Uhrf1 loss induces the interferon-response genes in non-immune cells.

If DNA methylation changes in uhrf1 mutants directly contribute to the induction of immune gene expression then a difference in the methylation pattern should be observed in the regulatory regions of these genes. Whole-genome bisulfite sequencing of uhrf1 mutants revealed that all genomic elements evaluated had reduced 5MeC levels (Feng et al., 2010); however, this low-coverage approach to measure DNA methylation was insufficient to evaluate the methylation profile at specific loci. We used Sanger sequencing of bisulfite-treated DNA to examine the CpG-rich regions in the promoter (−5 kb to +2 kb relative to the transcription start site) of ifnphi1, which is upregulated in uhrf1 mutants, as compared with ifnphi4, which is not (Fig. 2D). As a positive control, we assessed KenoDr1, a long interspersed nuclear element (LINE) that is heavily methylated in control embryos (Feng et al., 2010). All CpGs analyzed in KenoDr1 were over 90% methylated in WT embryos and significantly hypomethylated in uhrf1 mutants at 120 hpf (Fig. S8A), confirming that TE methylation is significantly reduced in this model. By contrast, all CpGs analyzed in the promoters of ifnphi1 and ifnphi4 were completely unmethylated in WT embryos, and remained unmethylated in uhrf1 mutants (Fig. S8B). Methylome datasets from WT embryos (Lee et al., 2015; Zhou et al., 2014) showed that these genes lack methylation at 24 hpf as well (Fig. S8C). Thus, there is no correlation between the expression of interferon genes and methylation of their promoters. Although we cannot exclude the possibility that other regulatory regions are affected by global loss of DNA methylation, as suggested by recent studies (Lee et al., 2015), our data argue against a direct regulatory role of promoter DNA methylation in the upregulation of interferon genes in uhrf1 mutants.

All cells have the potential to sense and respond to foreign cytosolic nucleic acids, and there are multiple sensors for non-self RNA and DNA (Dempsey and Bowie, 2015). Phosphorylation of Tbk1 (pTbk1) serves as an indicator that viral sensor signaling is activated (Ma et al., 2012). pTbk1 was barely detectable in WT larvae at 120 hpf, but was significantly induced in uhrf1 mutants (Fig. 6A,B). Since multiple antiviral pathways converge on pTbk1, we reasoned that inhibiting its activity would have a greater effect on downstream signaling than inhibiting one sensor at a time. A previous study reported that 1 µM BX795 suppressed pTBK1 signaling by inhibiting its ability to phosphorylate IRF3 (Clark et al., 2009). Here, we found that uhrf1 mutant embryos treated with 1 µM BX795 from 48-120 hpf had significantly reduced expression of nearly all of the interferon-response genes in our panel (Fig. 6C), whereas BX795 treatment had no effect on immune gene expression in WT larvae (Fig. S9A).

To determine whether the RNA- or DNA-sensing arm was activated in uhrf1 mutants, we analyzed our RNA-seq data and found upregulation of gene expression in both arms. For instance, mavs was unregulated 1.4-fold (log₂) in both uhrf1 and dnmt1 mutants at 120 hpf, whereas sting (tmem173) was significantly upregulated by 3.97-fold and 2.14-fold (log₂), respectively. We investigated the contribution of Sting signaling to the immune gene response in uhrf1 mutants using a previously characterized sting morpholino (Ge et al., 2015). This significantly suppressed all the interferon genes in uhrf1 mutants to nearly WT levels (Fig. 6D, Fig. S9B). We conclude that signaling via the cytoplasmic DNA-sensing pathway is required for induction of the interferon response in uhrf1 mutants. Further work is needed to assess the contribution of the viral RNA-sensing pathway in this system.

Retrotranscribed DNA in the cytoplasm serves as a signal to activate the STING/TBK pathway. We asked whether this contributed to immune induction in uhrf1 mutants using Foscarnet (Fos), an inhibitor of reverse transcriptases and polymerases that is an effective antiviral agent against HIV and Herpes virus (Crumpacker, 1992; Das et al., 2016; Marchand et al., 2007; Vashishtha and Kuchta, 2016) and has also been shown to reduce the activity of retrotransposons in yeast (Hage et al., 2014). Consistently, Fos treatment significantly reduced the expression of the genes that were most highly upregulated in uhrf1 mutants (ifnphi1, irf1b, isg15; Fig. 6E, Fig. S9C). Importantly, Fos treatment also reduced the enrichment of macrophages in the head of uhrf1 mutants (Fig. 6F). This points to the intriguing possibility that cytoplasmic DNA produced via retrotranscription serves as the trigger signal for cytoplasmic antiviral signaling in uhrf1 mutants.

Activation of cytosolic viral pathways coupled with a strong induction of type I interferon expression suggest the presence of an infectious agent, yet the uhrf1 mutants used in this study were not selectively infected. An alternative explanation is that intragenomic parasites, i.e. TEs, had become activated and mimicked viral infection in mutants lacking DNA methylation. We tested this by first examining the expression of zferv, a TE previously reported to be expressed specifically in the thymus of WT zebrafish larvae at 5 days post fertilization (dpf) (Shen and Steiner, 2004). We discovered that 55% of uhrf1 mutants expressed zferv in the head (Fig. 7A,B), similar to where irf1b and isg15 were upregulated (Fig. 5), indicating that TE expression and immune induction occur at the same time and place in uhrf1 mutants.

To examine TE expression more broadly, we generated an RNA-seq library that was not poly(A) selected to allow us to capture TE expression from uhrf1 mutants at 55 and 120 hpf and from dnmt1 mutants at 120 hpf (Fig. 7C,D, Table S1). There were 90 and 92 transposons significantly differentially expressed in 55 and 120 hpf uhrf1 mutants, respectively (Fig. 7C, Table S1). The majority of the upregulated TEs belonged to class I retrotransposons, whereas downregulated transposons were enriched for class II DNA transposons (84%) (Fig. 7D). Genomic TE abundance in each category cannot account for this distribution, as DNA transposons (class I) dominate, occupying 38% of the zebrafish genome as compared with 10% of the genome occupied by retrotransposons (class II) (Howe et al., 2013). The pattern of TE induction was the same in dnmt1 mutants (Fig. 7D, Fig. S10B, Table S1). These data demonstrate that TE dysregulation is a transcriptional response to DNA hypomethylation during development.

If the model that DNA methylation causes derepression of TEs leading to the induction of immune gene expression is valid, then these events should occur in that chronological sequence. RNA-seq analysis of uhrf1 mutants at 55 hpf revealed a marked upregulation of TEs: out of the 90 TEs differentially expressed at this time point, the majority were upregulated (69%), with nearly all (95%) of the upregulated TEs annotated as retrotransposons (Fig. 7D), mirroring findings in uhrf1 and dnmt1 mutants at 120 hpf (Fig. 7D, Table S1). Although the induction of retrotransposons in uhrf1 mutants at 55 hpf (Fig. 7E) is lower than that detected at 80 hpf (Fig. S10A) or 120 hpf (Fig. 7E, Fig. S10B), it occurs soon after global DNA hypomethylation is detected (Fig. 3B) and is accompanied by a modest interferon response that becomes amplified as time progresses. This suggests that TE induction is the trigger for immune gene induction in embryos with global DNA hypomethylation.

We further evaluated the relationship between TE methylation and expression through bisulfite Sanger sequencing. Analysis of locus-specific transposon methylation of many TEs was confounded by the high divergence between transposon sequences across their multiple sites of integration in the genome, divergence in TE populations between individuals and the divergence in the TE sequences between the strains of zebrafish used in this study (AB, TAB14 and TAB5), as compared with the reference genome that we used to design primers and analyze the sequencing results (Tubingen). Nevertheless, we were able to amplify and sequence two TEs that were overexpressed in uhrf1 mutants at different time points – Gypsy-21 at 55 hpf and Gypsy-10 at 120 hpf – and found both to have significantly reduced methylation in uhrf1 mutants (Fig. 7F,G). Since retrotransposons are both hypomethylated and overexpressed at 55 hpf, preceding the full activation of the interferon gene panel in uhrf1 mutants (Fig. 8A), and since preventing retrotransposon expression via inhibition of retrotranscription blocked interferon gene induction, we conclude that retrotransposon derepression and retrotranscription trigger an interferon response in embryos with DNA hypomethylation (Fig. 8B).

The data here show that loss of DNA methylation in developing embryos is the first step in a series of events (Fig. 8A) that leads to the derepression of retrotransposons, which become retrotranscribed and activate the cytoplasmic DNA sensor Sting and its partner pTbk1 (Fig. 8B) to induce interferon expression. The association between DNA hypomethylation and the induction of the interferon response has recently been reported in several systems (Chiappinelli et al., 2015; Jackson-Grusby et al., 2001; Leonova et al., 2013; Matthews et al., 2011; Ramesh et al., 2016; Roulois et al., 2015; Sharif et al., 2016); our work showing that both retrotranscription and Sting are required to induce interferon genes by DNA hypomethylation are novel findings. We speculate that global DNA hypomethylation activates the interferon response as a mechanism to prune cells with widespread epigenetic aberrations from the developing embryo. As such, the immune system might function both to flag and eliminate cells with foreign pathogens and those with other alterations that pose a danger to the organism.

The relationship between the methylation status of self DNA and immune activation is not well understood. We identified TE derepression and the activation of cytosolic antiviral signaling as a response mechanism to stimulate interferon production in cells with hypomethylated self DNA. Transcriptome analysis of uhrf1 and dnmt1 mutants uncovered a gene expression profile highly reminiscent of those caused by viral and bacterial infections (Briolat et al., 2014; Levraud et al., 2007; Meijer and Spaink, 2011; Palha et al., 2013). This was surprising, because nearly all of the work in the field linking hypomethylated DNA to the interferon response has focused on exogenous DNA from either an infectious or experimental source or has suggested that DNA methylation directly regulates the expression of interferon genes. In the first scenario, engulfment of hypomethylated DNA through the endosomal pathway triggers TLR9, which activates the immune system and interferon production (Kužnik et al., 2011; Yasuda et al., 2009). Our attempts to block Tlr9 signaling in uhrf1 mutants have not been successful and thus it remains possible that hypomethylated self DNA, perhaps from dead cells, is engulfed and activates Tlr9 signaling in neighboring cells. Although the expansion of leukocytes is a prominent phenotype of embryos with DNA hypomethylation, this does not account for the induction of immune genes in uhrf1 mutants. Instead, we speculate that leukocytes may be recruited by hypomethylated DNA and debris released from dead cells. This is consistent with recent findings that uhrf1 and dnmt1 mutants have increased neutrophils and tnfa induction in the intestine (Marjoram et al., 2015), where there is an abundance of dead and dying cells. Indeed, our finding that blocking reverse transcription reduces macrophage expansion suggests that the macrophages are activated in response to TE activation.

Our findings indicate that DNA methylation causes induction of immune gene transcription in uhrf1 mutant cells independently of any direct regulatory function of DNA methylation in the expression of these genes. Indeed, there are few, if any, CpG islands in proximity to any interferon genes in zebrafish, mouse and human; those that are present are largely unmethylated, consistent with the consensus in the field that CpG islands are protected from methylation (Jones, 2012). The possibility exists that methylation-sensitive putative enhancers (Lee et al., 2015) exert some level of control over the genes that we assayed. Although we cannot rule out that enhancers or other cis-acting regulatory elements of ifnphi1 outside of what we examined here are differentially methylated in uhrf1 mutants, thus far there is no support for this idea in our model. Instead, our data agree with studies showing that most genes are not repressed by DNA methylation, since fewer than 1% of all genes were significantly induced in uhrf1 or dnmt1 mutants, or in other systems in which DNA methylation is globally reduced (Gutierrez-Arcelus et al., 2013; Jackson-Grusby et al., 2001; Xie et al., 2011).

By contrast, our data are consistent with the overwhelming evidence that repression of potentially mutagenic and genotoxic retrotransposons is a central and conserved function of DNA methylation (Chiappinelli et al., 2015; Coit et al., 2013; Hutnick et al., 2010; Roulois et al., 2015; Wen et al., 2007; Yeh et al., 2013; Yoder et al., 1997; Grow et al., 2015). Although transcription of endogenous retrovirus (ERV) long terminal repeats (LTRs) has previously been reported in association with DNA hypomethylation, only functional retrotransposons would be able to reverse transcribe their RNA in the cytoplasm. We propose that uhrf1 mutants undergo reactivation of latent TEs, which contributes to the expansion of macrophages and potentially could also contribute to other phenotypes. For instance, it is possible that once active and retrotranscribed, intact retrotransposons could reintegrate into the genome, causing genomic instability.

A similar mechanism of ERV activation has recently been reported to cause an interferon response via viral mimicry in cancer cells; however, that study did not explore activation of the DNA-sensing arm of the cytoplasmic antiviral cascade nor the possibility of active ERV reverse transcription (Chiappinelli et al., 2015; Roulois et al., 2015). Our study confirms and extends this by demonstrating that an analogous mechanism operates in non-transformed cells and in a whole organism, both features that are essential for translating these findings to clinically relevant fields. Moreover, its possible that some of the immune gene induction we observe in uhrf1 mutants is in response to translated viral proteins encoded by the overexpressed TEs. The finding that actively expressed and translated ERVs in human blastocysts serve as a protective mechanism against viral infection (Grow et al., 2015) suggests the potential to utilize this response to promote immunoprotection in normal cells.

How much methylation is required to repress TE expression is a key unanswered question in the field. We observe 5MeC reduction of as little as 8.4%, as in the case of Gypsy-21, and conclude that this is sufficient to derepress this transposon. However, our method of bulk methylation analysis cannot address the cell-specific changes in 5MeC that could be impacting TE expression. Thus, the level of hypomethylation can vary greatly from cell to cell. For example, in one study that used a transgenerationally silenced CpG-rich UAS promoter in zebrafish, a reduction to 69-79% from the original 90% CpG methylation was correlated with derepression (Goll et al., 2009). Given that TE expression also depends on the age of the TE, which is inversely correlated with both methylation (Hutnick et al., 2010; Jackson-Grusby et al., 2001) and the capacity for activation (Huang et al., 2012), the presence of functional regulatory sequences, which are largely uncharacterized in zebrafish, and the presence of other repressive epigenetic marks, a simple correlation between methylation and expression is difficult to decipher.

This work has relevance beyond infection, as autoimmune diseases and cancer are characterized by heightened immune responses. Interestingly, both of these diseases are also characterized by DNA hypomethylation: T-cells from human systemic lupus erythematosus patients have extensively hypomethylated DNA and a hyperactive type I interferon response (Absher et al., 2013; Brooks et al., 2010; Volkman and Stetson, 2014; Wen et al., 2007) and treating T-cells in vitro with DNA demethylating agents phenocopies some lupus symptoms (Wen et al., 2007). In cancer, the entire genome becomes hypomethylated, which both promotes malignancy and induces an immune response against cancer cells. We identified DNA hypomethylation as a driver of liver cancer in zebrafish with hepatocyte-specific UHRF1 overexpression (Mudbhary et al., 2014). Since chromosomal instability is a hallmark of many cancers, it is possible that reactivation of functional TEs leads to an increased rate of mutation and contributes to the development of aggressive cancers. Indeed several studies point in this direction (Honda, 2016; Kemp and Longworth, 2015; Solyom and Kazazian, 2012). Therefore, the canonical role of DNA methylation to suppress TE activation is essential for both normal vertebrate development and for preventing disease.

Fish were raised in accordance with the policies of the Mount Sinai Institutional Animal Care and Use Committee (IACUC) on a 14:10 h light:dark cycle at 28°C. Embryos carrying the hi272 allele of uhrf1 mutants were used as described (Jacob et al., 2015; Sadler et al., 2007). Embryos homozygous for the hi272 allele (referred to as uhrf1 mutants) (Amsterdam et al., 2004) were generated from an incross of uhrf1^hi272/+ parents. All animal research involving zebrafish was approved by the IACUC of Icahn School of Medicine and of New York University. The Tg(mpeg:mCherry) and Tg(lysZ:dsRed) lines were described previously (Ellett et al., 2011; Hall et al., 2007). dnmt1^s904/s904 mutant embryos were obtained from crossing heterozygous parents (Anderson et al., 2009) and sorted based on morphology. All morphological, histological and gene expression assays for dnmt1 mutants were performed in the same manner as for uhrf1 mutants.

Individual embryos were genotyped at 48-60 hpf, which is before uhrf1 mutants can be identified based on phenotype, using DNA from whole embryos if they had previously been fixed in 4% paraformaldehyde (PFA) or on tails if further tissue collection was required for RNA isolation. Tissue was collected into 30 µl genotyping DNA lysis buffer (10 mM Tris-HCl pH 7.5, 50 mM KCl, 0.3% Tween 20, 0.3% NP-40) per embryo and genotyped using the primers listed in Table S5 as described (Amsterdam et al., 2004).

RNA was extracted from three to ten embryos using Trizol (Thermo Fisher, 15596026) and precipitated with ethanol. For embryos younger than 5 dpf, only the anterior half was used for RNA extraction to enrich for tissue most affected by uhrf1 loss. RNA used to assess expression of TEs was DNase treated and re-extracted in Trizol as above. For DNA preparation, five to ten embryos were incubated overnight at 55°C in DNA lysis buffer (10 mM Tris-HCl pH 7.5, 5 mM EDTA, 1% SDS) and extracted with phenol/chloroform followed by ethanol precipitation. 1 μg RNA was used with qScript (QuantaBio, 95048-025) to generate cDNA.

RT-qPCR was performed on cDNA from embryos at 1-5 dpf using primers (Table S5) designed for immune genes or to verify some of the top 50 upregulated genes identified from the microarray. Expression in all samples was normalized to rplp0. All genes were analyzed in at least three clutches.

Genome-wide expression profiling was performed using the Zebrafish Genome Array (Affymetrix) according to the manufacturer's instructions. Scanned raw data were normalized using robust multiarray analysis (RMA) algorithm implemented in the GenePattern genomic analysis toolkit (broadinstitute.org/cancer/software/genepattern/) (Reich et al., 2006). Multiple probes corresponding to a single gene were collapsed into the official gene symbol provided by NCBI by extracting a probe with maximal variation. Orthologous conversion to human genes was performed based on the mapping table provided by Ensemble Biomart (http://www.ensembl.org/biomart/martview/da9e5e05cf511b144e8e81ed57fcfe80). Differentially expressed genes were determined using Bayesian t-test implemented in Cyber-T software from the top 1000 genes with the largest coefficient of variation. Posterior probability of differential expression (PPDE) >0.95 was regarded as statistically significant. As described by Jacob et al. (2015), 117 upregulated and 131 downregulated genes were detected in uhrf1 mutants compared with WT siblings. Induced or suppressed molecular pathways from the microarray dataset described by Jacob et al. (2015) (GSE55339) were determined using GSEA (Subramanian et al., 2005) implemented in the GenePattern genomic analysis toolkit and Molecular Signature Database (MSigDB, broadinstitute.org/cancer/software/genepattern/) (Jacob et al., 2015).

RNA-seq libraries were prepared according to the Illumina TruSeq RNA sample preparation version 2 protocol with Ribo-Zero Gold. RNA from a pool of embryos (between 10 and 20) were used to generate libraries, which were analyzed on an Agilent 2100 Bioanalyzer. We analyzed three clutches of uhrf1 mutants and their phenotypically WT siblings at 55 hpf, two clutches at 120 hpf and three clutches of dnmt1 mutants and siblings at 120 hpf. cDNA libraries for 120 hpf uhrf1 mutants were sequenced on the Illumina NextSeq500 platform to obtain 75 bp single-end reads, while 55 hpf uhrf1 and 120 hpf dnmt1 mutants were run on the HiSeq platform to obtain 100 bp paired-end reads. Sequencing quality was assessed using FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc) and the reads were quality trimmed using Trimmomatic (Bolger et al., 2014) to remove low Q scores, adapter contamination and systematic sequencing errors. Reads were aligned to the Danio rerio GRCz10 reference genome assembly with TopHat 2.0.9 (Trapnell et al., 2009). To estimate gene expression, both FPKM and read counts were calculated with Ensembl annotation (Aken et al., 2016; Anders et al., 2015; Trapnell et al., 2009). Annotation of TEs is according to RepeatMasker track in the UCSC table browser. Reads of each TE were counted by HTseq (Anders et al., 2015) and normalized to specific sample size. Test of differential expression uses a generalized linear model. Counts of each TE were modeled as negative binomial distributions, which were both implemented in DESeq2 in Bioconductor (Gentleman et al., 2004). An adjusted P-value [false discovery rate (FDR)] <0.05 was considered as significantly different expression. Transposon families were identified via RepBase (Jurka et al., 2005).

We treated embryos with 20 µg/ml Neutral Red stain in embryo water for 1 h followed by a 1 h washout in embryo water to allow the background dye to diffuse out. We increased the concentration of dye and decreased the incubation time from that previously reported in order to obtain maximum intensity of macrophage staining while only staining the most active population of macrophages, since that population would take up the dye more quickly.

Bisulfite (BS) conversion of 500 ng genomic DNA was carried out using the EpiJET Bisulfite Conversion Kit (Thermo Scientific, K1461) as per manufacturer's instructions. As is the case with BS conversion, extensive degradation and loss of DNA prevents accurate concentration measurements using conventional spectrophotometers. Instead, 1 μl of each BS-DNA sample was pipetted onto a thin layer (0.5 cm) of 0.8% agarose containing a 1:1000 dilution of GelStar (Lonza, 50535) and allowed to be absorbed. The fluorescence from each sample was compared with serially diluted DNA standards pipetted onto the same plate and the amount of each sample used in PCR was adjusted accordingly. Bisulfite-specific PCR (BS-PCR) primers (Table S5) were designed using MethPrimer software (Li and Dahiya, 2002).

Antisense probes recognizing a 350 bp region of the irf1b transcript and 500 bp of the isg15 transcript were generated by PCR from 80 hpf uhrf1 mutant cDNA using gene-specific primers (Table S5) with a T7 promoter sequence at the 5′ end and a T3 promoter at the 3′ end. A 1.9 kb fragment of zferv was amplified from the zferv PCRII plasmid (Addgene, 22399) and in vitro transcribed from the T7 promoter. ISH was carried out as previously described (Thisse et al., 2004).

uhrf1 mutants were identified based on phenotypic features after 80 hpf or were individually genotyped at earlier stages and the DNA extracted using a standard protocol as above and probed for 5MeC by slot blot essentially as described (Jacob et al., 2015; Mudbhary et al., 2014). Briefly, 3 ng DNA was denatured in 0.4 M NaOH/10 mM EDTA, neutralized with 2 mM ammonium acetate and loaded in duplicate onto a nitrocellulose membrane using a slot blot apparatus. Membranes were baked at 80°C, blocked with 5% milk followed by incubation in either anti-5MeC (Eurogentec, BI-MECY-100; 1:2000) or anti-dsDNA (Abcam, ab27156; 1:8000) overnight, washed in TBST (37 mM NaCl, 20 mM Tris pH 7.5, 0.1% Tween 20) and probed with anti-mouse HRP secondary antibody (Promega; 1:5000) for 1 h at room temperature followed by development in ECL (Thermo Scientific). Total 5MeC and dsDNA was averaged between duplicates, and the 5MeC:dsDNA ratio was calculated for at least three clutches at each time point and averaged.

At 72 hpf, Tg(mpeg:mCherry) larvae were anesthetized in 0.016% tricaine and microinjected in the left otic vesicle with 1 nl PBS containing either CpG oligodeoxynucleotide (15 µM, 100 ng/µl), methylated CpG (mCpG) oligodeoxynucleotide (15 µM) or no oligodeoxynucleotide (control). The sequence of the oligos is 5′-TCGTCGTTGTCGTTTTGTCGTT-3′, with the CpG oligo unmethylated and the mCpG oligo methylated at all cytosine residues (Yeh et al., 2013). Macrophage presence in the otic vesicle was assessed by fluorescent confocal microscopy at 2 and 4 h post injection in agarose-mounted larvae.

WT and uhrf1 mutant embryos were incubated in 5 ml 1 µM BX795, 500 ng/ml Foscarnet or 0.1% DMSO from 48-120 hpf. Owing to instability, a fresh 1 µM BX795 solution was used each day. At 120 hpf, five to ten larvae were collected for RNA extraction.

Morpholinos (Table S5) were obtained from GeneTools and injected into 1-cell embryos at an average of 4 nl per embryo using stock concentrations of 0.6 mM, 0.5 mM or 0.3 mM.

Embryos were collected in 10 µl protein lysis buffer per embryo and homogenized by sonication. Samples were run on an 8% SDS-PAGE gel, transferred to PVDF membrane and blotted with anti-pTBK1 (Cell Signaling, 5483; 1:1000). Anti-tubulin (Developmental Studies Hybridoma Bank; 1:5000) and anti-β-actin (Sigma, A2228; 1:5000) were used as loading controls.

Whole-mount embryos were anesthetized and mounted in 3% methyl-cellulose. Imaging was carried out using a Nikon SMZ25 stereomicroscope.

Analyses of continuous and categorical variables between groups were compared using the Wilcoxon rank-sum test and Fisher's exact test, respectively. Two-tailed P<0.05 was regarded as statistically significant. In experiments involving treatment of embryos to abrogate the immune response, statistical significance was determined based on the percentage residual gene induction between the treated (either drug or morphant) and control embryos. Multiple hypothesis testing was adjusted using FDR at P<0.05 or Bonferroni correction as appropriate. Band intensity and macrophage activity quantifications were performed using GelAnalyzer (http://www.gelanalyzer.com/) and ImageJ (NIH), respectively. CpG methylation analysis was carried out using QUMA (Kumaki et al., 2008). Prism6 software (GraphPad) was used for figure generation and analysis. Venn diagrams and analysis of overlapping genes between the microarray and RNA-seq datasets were generated using BioVenn (Hulsen et al., 2008).

We thank Patrick Bradley and Matthew Nash for expert fish care; Amaia Lujambio for research insight and editorial critique; Jessica Lau for technical contributions; and Jill Gregory for scientific illustration.