Zebrafish work was conducted at the Garvan Institute of Medical Research in accordance with the Animal Ethics Committee AEC approval and with the Australian Code of Practice for Care and Use of Animals for Scientific Purposes. Adult wild type (AB/Tübingen) Danio rerio (zebrafish) were bred in an equal ratio of males and females. Embryos were collected 0 h post-fertilization (hpf) and incubated in 1X E3 medium (5 mM NaCl, 0.33 mM CaCl2, 0.17 mM KCl, 0.33mM H14MgO11S) for 4 days at 28.5°C before being transferred onto a filtered system.

Whole brains were dissected from zebrafish larvae and adults before being snap-frozen in liquid nitrogen and stored at −80°C. Genomic DNA (gDNA) was extracted from brains using the QIAGEN DNeasy Blood & Tissue Kit (QIAGEN, Chadstone, VIC, Australia) according to manufacturer’s instructions. For RNA extraction, half of the lysate from the first step of DNA extraction from the QIAGEN DNeasy Blood & Tissue Kit was added to TRIsure (Bioline) and purified following manufacturer’s instructions. All experiments in this study were performed in two biological replicates.

Guide RNAs (gRNA) targeting dnmt3aa and dnmt3ab loci were designed with CRISPRscan (Moreno-Mateos et al., 2015). gRNAs for both loci were synthesized and co-injected into 1-cell stage embryos as previously described (Ross et al., 2020). CRISPR/Cas9 knockouts (cKO) fish were grown to 4 weeks of age before their brains were harvested for DNA and RNA extraction. Amplicons surrounding the CRISPR/Cas9 cut sites were PCR-amplified from genomic DNA, ligated to NEXTFLEX Bisulfite-Seq barcodes (PerkinElmer, Waltham, MA, United States), and spiked into libraries that were sequenced on the Illumina HiSeq X platform. Knockout efficiencies were calculated from the sequenced amplicons using CRISPResso (Pinello et al., 2016). RNA was reverse transcribed to cDNA using the SensiFAST^TM cDNA Synthesis Kit (Bioline), following the manufacturer’s protocol and subjected to qPCR analysis. Relative expression levels were calculated using the 2−ΔΔCT method and bactin gene as the control transcript. Two sample t-tests were performed on CT values. All oligos used in this study can be found in Supplementary Table 1.

WGBS libraries were prepared from 500 ng of zebrafish brain gDNA, spiked with 0.025 ng of unmethylated lambda phage DNA (Promega, Madison, WI, United States) and sequenced on the Illumina HiSeq X platform (2 × 150 bp) as previously described (Ross et al., 2020).

RRBS libraries were prepared from 500 ng of zebrafish brain gDNA spiked with 0.025 ng of unmethylated lambda phage DNA (Promega, Madison, WI, United States). gDNA was digested with 10 U BccI (CCATC(N)₄) and 10 U SspI (AATATT) for 2 h. A separate aliquot was digested with 20 U MspI (New England BioLabs, Ipswich, MA, United States). RRBS libraries were constructed as previously described (Ross et al., 2020), sequenced, and the BAM files corresponding to different digestion reactions (BccI/SspI or MspI) were merged before downstream analysis.

WGBS reads were trimmed with Trimmomatic: ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 SLIDINGWINDOW: 5:20 LEADING:3 TRAILING:3 MINLEN:60 (Bolger et al., 2014), and mapped using WALT (-m 5 -t 20 -N 10000000) (Chen et al., 2016) onto the GRCz11 reference genome (UCSC), containing the λ genome. BAM files, containing only uniquely mapped reads, were deduplicated using sambamba markdup (Tarasov et al., 2015). RRBS data were processed as above with the additional option of: HEADCROP:5 CROP:140 added during trimming and without deduplication. BAM files were made FLAG-compatible and processed with CGmapTools (Chen et al., 2016; Guo et al., 2018) (convert bam2cgmap) to obtain ATCGmap files, which were corrected for CH positions that showed evidence of CG SNPs (de Mendoza et al., 2021). A summary of library statistics can be found in Supplementary Table 2. Genomic data were visualized in UCSC (Kent et al., 2002) and IGV (Robinson et al., 2011) browsers.

BED file coordinates of the 10,000 most highly methylated mCH and mCAC sites, with a minimal depth of 10, were extended by 4 base pairs upstream and downstream. The resulting files were used as input for HOMER “findMotifsGenome.pl” function (Heinz et al., 2010), establishing the search for de novo motifs of length 9 (-len 9 -size given) with the GRCz11 genome used as the background sequence. Motifs were visualized using the “ggseqlogo” package in R (Wagih, 2017) and motif positions in the genome were called using the scanMotifGenomeWide.pl function (with and without -mask option checked).

Bedgraphs were generated from corrected CGmapTools outputs and converted to bigWig using bedGraphToBigwig script from Kent utils. Average mCH levels were determined from bedGraph files and calculated by dividing the sum of reads supporting a methylated cytosine by the sum of all reads mapping to that position. mCH levels in genomic features and gene bodies were calculated using BEDtools map (Quinlan and Hall, 2010). mCH levels, TPMs, and gene length were plotted using the boxplot function in R (outline = FALSE). Heatmaps were generated using deepTools (Ramirez et al., 2014) computeMatrix with the following parameters: “computeMatrix scale-regions -m 650 -b 500 -a 500 -bs 25.” NAN values were replaced with 0. Heatmaps were plotted with the plotHeatmap function, sorted, and clustered based on methylation levels. Scatterplots were generated using the geom_bin2d function in ggplot2 [(bins = 50) + geom_smooth(method = l m)]. Pearson correlations were calculated using the rcorr function in R.

Differentially methylated regions (DMRs) were called using DSS (delta = 0.1, p.threshold = 0.05, minlen = 100, minCG = 5, dis.merge = 500, pct.sig = 0.5) (Feng et al., 2014).

Repeatmsker tracks were obtained from UCSC. The percentage of repeat subfamilies overlapping the top-methylated CAC motifs was determined with BEDtools (intersectBed).

Sequenced ChIP-seq reads (Kaaij et al., 2016; Yang et al., 2020) were trimmed with Trimmomatic (ILLUMINACLIP:TruSeq3-SE.fa:2:30:10 SLIDINGWINDOW:5:20 LEADING:3 TRAILING:3 MINLEN:20) (Bolger et al., 2014) before being mapped to the GRCz11 genome using Bowtie2 with default settings (Langmead and Salzberg, 2012). BAM files were deduplicated using sambamba markdup (Tarasov et al., 2015). Peaks were called using MACS2 (Zhang et al., 2008). bigWigs were generated using deepTools (Ramirez et al., 2014) bamCompare (-e 300 -p 20 –normalizeUsing RPKM –centerReads), and bedGraphs using UCSC bigWigToBedGraph. Spearman correlations were calculated using BEDtools map (Quinlan and Hall, 2010) and rcorr(type = c(“spearman”) from the Hmisc package in R.

RNA-seq libraries were prepared with 1,000 ng of input RNA using the KAPA mRNA HyperPrep Kit, according to manufacturer’s instructions and sequenced on the Illumina HiSeq X platform (2 × 150 bp).

RNA-seq reads were trimmed using Trimmomatic: ILLUMINACLIP: TruSeq3-PE.fa:2:30:10 SLIDING-WINDOW:5:20 LEADING:3 TRAILING:3 MINLEN:60 (Bolger et al., 2014) and aligned to the GRCz11 genome using STAR (Dobin et al., 2013). Differential gene expression analysis was performed using edgeR (Robinson et al., 2010). Robinson et al., 2010 with genes selected based on a minimum ± 1.5 logFC (FDR < 0.05) between any of the analyzed time points. Z-scores were calculated based on the log2 transformations of TPM values and plotted with the pheatmap package in R. Analysis of published RNA-seq data was performed based on the provided read count tables with TPM values calculated from the average of 5–6 month old brain datasets (Aramillo Irizar et al., 2018) or from collated counts of 30 neurons (Lange et al., 2020).

To investigate mCH in the zebrafish nervous system, we analyzed WGBS data (bisulfite conversion > 99.5%) of adult brain (mean coverage = 9.9X), as well as of adult liver (mean coverage = 7.6X), to use as a non-neural control tissue (Bogdanovic et al., 2016; Skvortsova et al., 2019). We employed stringent genotype correction (de Mendoza et al., 2021) to allow for more sensitive interrogation of mCH patterns. To better understand the sequence context of mCH deposition in the zebrafish brain and how it compares to non-neural tissues (liver), we performed a motif search on the 10,000 most highly methylated CH sites. Expectedly, we recovered the TGCT motif associated with mosaic satellite repeats in both brain and liver (Ross et al., 2020), and the previously characterized TACAC-containing motif (Lister et al., 2013; de Mendoza et al., 2021), which was specific to the brain sample (Supplementary Figure 1A). Furthermore, the brain sample showed a notable enrichment in CAC trinucleotide methylation (∼1.5%) with a threefold increase compared to the unmethylated lambda genome spike-in control (∼0.5%) (Supplementary Figure 1B). This mCAC enrichment was not evident in the liver sample, thus eliminating the possibility of sequence-specific biases or artifacts pertaining to bisulfite conversion (Olova et al., 2018). We next investigated the genomic distribution of mCAC in zebrafish brains to assess if the depletion in regulatory elements and enrichment in gene bodies previously described in mammals (Lister et al., 2013; Guo et al., 2014a; He and Ecker, 2015) is evolutionarily conserved. To achieve this, we annotated transcription start sites (TSS), exons, introns, 5′UTRs, 3′UTRs, intergenic regions, and sites of H3K27ac enrichment, which correspond to active gene-regulatory elements (Kaaij et al., 2016). We found that intronic and intergenic regions are the only regions enriched in mCAC and mCH (Supplementary Figures 1C–E), whereas a notable depletion, similar to the one described in mammals, was observed at H3K27ac peaks (Supplementary Figure 1F). These results support the previous observations of mCA presence in gene bodies of vertebrate brains (de Mendoza et al., 2021) and demonstrate a conserved depletion of mCH in active regulatory regions.

Further analyses of sequence motifs associated with the mCAC context unraveled two novel sequences in addition to the previously described vertebrate-conserved TACAC motif (Figure 1A). Due to previous associations of mCH with repetitive DNA in zebrafish (Ross et al., 2020), we wanted to assess if any of the most significantly methylated CAC motifs were enriched in repetitive elements. The top TACAC motif displayed comparable methylation levels in repetitive elements and in the repeat-masked genomic fraction, with an average methylation level nearly three-fold higher (∼4%) than the average global mCAC levels (Figure 1B). However, for the remaining two motifs we found a robust increase in average methylation levels at repetitive elements when compared to the repeat-masked genome (∼6.5% and 6%, Figure 1B). Further analysis revealed that the TACAC motif is broadly distributed in the genome whereas the second and third motif were mainly located in TDR and TC1DR3 repetitive elements, respectively (Figure 1C). Both TDR and TC1DR3 elements belong to the Tc1-mariner superfamily, which is found across eukaryotes. These elements are characterized by two inverted terminal repeats, and an open reading frame (ORF) (Wicker et al., 2007). The identified motifs are found on average once per element at positions 205/541 and 1005/1071 in TDR18 and TC1DR3 model sequences, respectively (Storer et al., 2021). To further validate these associations, we interrogated previously published heart and forebrain WGBS DNA methylomes (Ross et al., 2020; de Mendoza et al., 2021). These analyses again revealed brain-specific enrichment of mCH specifically at these elements (Supplementary Figure 1G), thus excluding the possibility of library preparation and sequencing artifacts causing the observed enrichment.

We also found that average methylation of the top three methylated CAC motifs correlated strongly with overall mCH in gene bodies (R = 0.52) (Figure 1D), suggestive of a significant contribution of CAC methylation to gene body mCH. This correlation was stronger than the one observed between average gene mCH and gene length (R = 0.36), which was previously described in mammals (Gabel et al., 2015; Boxer et al., 2020; Supplementary Figure 1H). Additionally, when Gene ontology (GO) analysis (Raudvere et al., 2019) of genes containing methylated TDR and TC1DR3 CAC motifs was performed, we found overrepresentation in developmental and neural development terms (Figure 1E). Due to its widespread genomic abundance the TACAC motif was omitted from the GO analysis. Overall, zebrafish brain is enriched in mCH, particularly in the CAC trinucleotide context, predominantly in introns and intergenic regions, as well as in members of the Tc1-mariner transposon family.

To explore the relationship between mCH and gene expression in the zebrafish brain, we plotted average gene mCH and mCAC values against gene expression (transcripts per million-TPM) levels from available datasets (Bogdanovic et al., 2016; Aramillo Irizar et al., 2018; Figure 2A). We revealed a gene cluster (n = 1,860) with higher-than-average mCH levels, which displayed lower expression than the bulk transcriptome (Figures 2A,B). To provide more genomic context to these findings, we first interrogated whether this elevated mCH was driven by a higher proportion of intron sequence in these genes. To that end, we plotted gene length, gene body mCAC, intron mCAC and exon mCAC for mCH-enriched genes as well as for genes that did not display any notable mCH enrichment (Figure 2C). While the mCH-enriched genes were significantly longer (Wilcoxon test, ^∗∗∗ P < 0.001), in line with observations in mammals (Gabel et al., 2015; Boxer et al., 2020), the elevation in mCH was not driven exclusively by intron contribution as both introns and exons located within these genes had significantly higher levels of mCAC (Wilcoxon test, ^∗∗∗P < 0.001) (Figure 2C). Interestingly, genes with higher levels of mCAC also contained a higher percentage of Tc1-like elements in relation to total gene length (Figure 2D). To confirm the observation of poorly expressed genes being marked by mCH, we plotted average TPM levels from total RNA-seq data from adult brains (Aramillo Irizar et al., 2018) and combined single cell data from neurons (Lange et al., 2020; Figure 2E). Genes with higher levels of mCH (cluster 1) had significantly lower average TPM when compared to genes with moderate/low levels of mCH (cluster 2) (Wilcoxon test, ^∗∗∗P < 0.001), or to a randomly selected subset of genes (n = 1,860, Figure 2E) (Wilcoxon test, ^∗∗∗P < 0.001). This difference in expression levels between the two clusters was even more pronounced in neurons where mCH is expected to be the highest based on mammalian data (Lister et al., 2013; Figure 2E). GO analysis of these highly mCH-methylated genes again revealed enrichment for terms associated with embryonic and neural development (Supplementary Figure 2A). This enrichment of developmental genes for mCH is also in line with observations in other vertebrate brains (de Mendoza et al., 2021).

As mCH-enriched genes are on average poorly expressed, longer, and associated with neuronal and development terms, it is yet unclear which features are most important or predictive for mCH enrichment. This is further complicated by neuronal terms and gene lengths being tightly associated in zebrafish, as discussed in our previous study (Ross et al., 2020). To further explore and rank features that may be associated with genomic mCH, we analyzed histone modification ChIP-seq data from the zebrafish brain (Yang et al., 2020) and assessed the correlations between these diverse histone marks and mCH levels in gene bodies (Supplementary Figure 2B). We also performed correlation analyses for mCH and various other genomic features (Supplementary Figure 2B). These analyses revealed that gene length has the strongest positive correlation with gene body mCH levels and that H3K4me3 at transcription start sites has the strongest negative correlation with gene body mCH. Overall, these results further demonstrate that mCH-enriched genes are longer and associated with a repressive chromatin environment. Therefore, like in mammals, mCH in the zebrafish adult brain is associated with transcriptional repression which is particularly evident in long genes.

As mCH has previously been shown to accumulate during mammalian brain development (Lister et al., 2013), we next investigated whether comparable mCH dynamics could be observed in zebrafish. We generated RRBS libraries (Meissner et al., 2005) using a combination of enzymes, which were selected based on virtual digestion, to enrich for regions containing highly methylated CAC motifs (Figure 1A) identified in adult brains. We assayed zebrafish brains starting from 1 week (1W), where brain structures such as the cerebrum are not identifiable, up until 6 weeks (6W) and adulthood, where all structures are discernable (Maeyama and Nakayasu, 2000). This analysis revealed a gradual increase in mCAC in the brains of 1-week to 6-week-old zebrafish followed by a more notable increase in adult brains (Figure 3A). RNA-seq analysis across this period also revealed a gradual decrease in the expression of components of DNA methylation machinery as cells presumably become more differentiated (Figure 3B). Moreover, mecp2 expression increased during brain development coinciding with increase in mCH. This observation supports a conserved role for MeCP2 in the regulation of genes marked by mCH in vertebrates (Guo et al., 2014a; Chen et al., 2015; Gabel et al., 2015; de Mendoza et al., 2021). Differential expression analysis of all genes across brain development revealed two major gene clusters which were either consistently upregulated (cluster 1) or downregulated (cluster 2) during brain development (Figure 3C). Interestingly, 18% of the downregulated genes belonged to the mCH-enriched gene cluster, compared to only 4% of the upregulated genes (χ²-test ^∗∗∗P < 0.001) (Figure 3C). Downregulated genes were also associated with terms related to cell division and development while genes that were upregulated were enriched in terms associated with adaptive immunity (Figure 3D). This are consistent with ongoing developmental processes in the larval brain and with the notion that the adaptive immune system of zebrafish does not develop until 4 to 6 weeks post fertilization (Lam et al., 2004).

Since we determined that mCH dynamics in the developing zebrafish brain are comparable to developmental mCH dynamics in the postnatal mouse frontal cortex (Lister et al., 2013), we next wanted to explore if mCG dynamics similar to the one observed in mice could also be detected in zebrafish. Pairwise analyses of differentially methylated regions (DMRs) from RRBS data identified 1723 DMRs. The strongest difference in mCG was observed between 1-week-old brain and adult brain (Supplementary Figure 3A). This analysis also revealed that the 1-week-old zebrafish brain is most similar to the fetal mouse sample, as both the 1-week-old zebrafish brain and fetal mouse brain displayed most obvious differences when compared to other time points (Lister et al., 2013). The number of DMRs identified in zebrafish (n = 1,723) significantly differs from the number of DMRs identified in developing mouse brains (n > 142,835). This discrepancy is likely caused by different 5mC detection approaches (RRBS instead of WGBS), however other contributing factors could include: the use of whole brains instead of sorted neuron and glial populations, developmental stages not being directly comparable between zebrafish and mouse, adult zebrafish brains retaining more “juvenile features” such as radial glial cells and neurogenic capabilities (Schmidt et al., 2013), and different glia/neuron percentages of the samples, as mCG levels in glia have been reported to be more stable during development (Lister et al., 2013).

To understand better how methylation and gene expression dynamics track over developmental time, we generated WGBS datasets for 4-week-old brain tissue and compared these data against adult brain WGBS and RNA-seq data. Analysis of mCAC levels of differentially expressed genes revealed that developmentally downregulated genes accumulate more mCH when compared to upregulated ones (Figure 3E). This trend was also observed when visualizing mCAC levels and gene expression levels across all genes (Supplementary Figures 3B,C). Furthermore, quantification of mCAC levels at all mCAC trinucleotides and at the highly methylated motifs (all three combined), confirmed a significant increase in the methylation of developmentally downregulated genes (Wilcoxon test, ^∗∗∗P < 0.001) (Figure 3F). This increase in mCH in adult brains is uncoupled from global mCG changes, as DMR analysis of these WGBS datasets revealed 23,992 DMRs with the majority (n = 18,522) becoming hypomethylated during nervous system development (Supplementary Figures 3D,E). Altogether, these results demonstrate robust anticorrelation between mCH and gene expression during brain development in zebrafish, as well as developmental mCH accumulation, in line with observations in mammals.

Finally, to investigate if Dnmt3a-dependent methylation of CAC trinucleotides is evolutionarily conserved in zebrafish, we generated dnmt3aa/dnmt3ab CRISPR/Cas9 double knockouts (cKO). qPCR analysis of cDNA extracted from brains of 4-week-old cKOs revealed their partial knockout status with a 50% reduction in expression levels for both dnmt3aa and dnmt3ab (Figure 4A). Sequencing of amplicons surrounding the CRISPR/Cas9 cut sites demonstrated comparable estimates of genome editing efficiency (Figure 4B). Finally, WGBS analysis of these samples demonstrated that depletion of dnmt3aa/dnmt3ab resulted in a significant (P < 0.05, Wilcoxon test) reduction in global mCAC levels as well as in a notable reduction (43%) of mCH at top mCAC motifs (P < 0.01, Wilcoxon test) (Figure 4C). This significant reduction in mCH observed already in these partial cKOs, which could have undergone possible selection for more wild type cells, suggests that Dnmt3aa and Dnmt3ab are the major enzymes responsible for neural mCAC deposition. The reduction in mCAC levels in these cKOs can be observed across the majority of gene bodies (Figure 4D) and on a genome-wide and locus-specific scale (Figure 4E). Notably, these perturbations in mCH did not result in any obvious morphological defects in the cKO fish (data not shown). Finer analysis of brain morphology, behavior, and potential inflammatory processes in these animals will be a focus of future studies. Additionally, no changes in global CpG methylation levels or significant DMRs could be detected (Supplementary Figure 3F) between the cKO and WT brain, suggestive of mCH deposition being a major function of Dnmt3a enzymes in zebrafish. This is contrary to what has been observed in mouse neurons, where Dnmt3a KO resulted in ∼10% decrease in global mCG/CG levels (Lavery et al., 2020). However, given the incomplete KO, the possible redundancies with other zebrafish DNMTs (Goll and Halpern, 2011), the propagation by Dnmt1 once established, and the use of bulk cell DNA methylome data, we cannot completely rule out a role for these enzymes in neuronal CpG methylation. Nevertheless, the clear role in mCH deposition, and the fact that Dnmt3a enzymes can be traced back to the root of vertebrates (de Mendoza et al., 2021), suggests conserved functional importance of mCH in the vertebrate nervous system.