The elephant shark methylome reveals conservation of epigenetic regulation across jawed vertebrates

Tissue source and DNA extraction

Elephant shark tissue samples were sourced as by-product of deceased animals harvested from commercial fishing in the Otago coastal region. As such, no animal ethics permission was applicable in this circumstance. No animal experimentation or manipulation was undertaken as defined by the Animal Welfare Act (2009, New Zealand), or according to guidelines issued by the New Zealand National Animal Ethics Advisory Committee (NAEAC, Occasional Paper No 2, 2009, ISBN 978-0-478-33858-4).

DNA was purified using a modified magnetic bead approach⁴⁴. Briefly, cells were first homogenised in “GITC” lysis buffer (4 M Guanidine thiocyanate, Sigma G6639; 50 mM Tris, Thermo 15568-025; 20 mM EDTA; Thermo 15575-020; 2% Sarkosyl, Sigma L9150-50G; 0.1% Antifoam, Sigma A8311-50ML), and this lysate mixture was then combined with TE-diluted Sera-Mag Magnetic SpeedBeads (GE Healthcare, GEHE45152105050250) and isopropanol in a volumetric ratio of 2:3:4, respectively. Following capture with a neodymium magnet, beads were washed once with isopropanol, twice with 70% ethanol and resuspended in filter-sterile milliQ water.

Preparation of WGBS-seq libraries

WGBS-seq was undertaken using a post-bisulfite adapter tagging (PBAT) method adapted from Peat et al., 2014⁴⁵. Briefly, 50–100 ng of purified DNA was subjected to bisulfite conversion using the Imprint DNA modification kit (Sigma, MOD50). Converted DNA underwent first strand synthesis with a biotin-labelled adapter sequence possessing seven random nucleotides at its 3’ end (BioP5N7, biotin- ACACTCTTTCCCTACACGACGCTCTTCCGATCTNNNNNNN). The product of first strand synthesis was captured using streptavidin-coated Dynabeads (Thermo, 11205D) and magnetic immobilisation. Double-stranded DNA was created using the immobilized first-strand as a template and an additional adapter that also possesses seven random nucleotides at its 3’ end (P7N7, GTGACTGGAGTTCAGACGTGTGCTCTTCCGATCTNNNNNNN). Unique molecular barcodes and sequences necessary for binding to Illumina flow-cells were added to libraries by PCR using 1X HiFi HotStart Uracil+ Mix (KAPA, KK2801 and 10 μM indexed Truseq-type oligos), with thermal cycling as follows: 12× (94°C, 80 sec; 65°C, 30 sec; 72°C, 30 sec).

For deep sequencing, libraries were sequenced with a single-end 100bp protocol on a HiSeq 2500 instrument (Illumina) using rapid run mode. For low-coverage sequencing of additional samples, libraries were sequenced on a MiSeq instrument (Illumina) until the desired depth (at least 15,000 mapped CG calls) was attained.

Detailed sequencing results are provided in Table S1.

Bioinformatic processing of WGBS-seq dataset

Mapped CG methylated calls for mouse liver⁵ were downloaded from GEO (accession GSE42836, sample GSM1051157) and analysed directly. For zebrafish muscle⁸ (SRA study SRP020008, run SRR800081) and sea squirt muscle¹⁴ (GEO accession GSE19824, sample GSM497251), raw sequencing data was downloaded and processed along with elephant shark WGBS-seq data generated in this study as follows.

Trimming was performed to remove both poor-quality calls and adapters sequences using TrimGalore (v0.4.0, default parameters). For the elephant shark data, 10bp were also removed from the 5’ end of reads to account for sequence biases associated with PBAT library construction.

Trimmed reads were aligned using Bismark⁴⁶ (v0.14.3, default parameters) with the --pbat option specified for elephant shark data. The following genome assemblies were used for alignment: zebrafish, GRCz10; elephant shark, 6.1.3; sea squirt, KH. For sea squirt and elephant shark, alignment was only performed against scaffolds larger than 277kb to avoid gene annotation issues and assembly artefacts. The deep-sequenced elephant shark data generated in this study was additionally mapped to the mitochondrial genome.

Bismark mapping reports were used to determine global methylation levels for low-coverage elephant shark data. All other datasets were deduplicated and CG methylation calls extracted using Bismark (--comprehensive and --merge_non_CG options specified).

The number of mapped cytosine calls for sequencing performed in this study are provided in Table S1. The frequency of non-CG methylation indicates the maximum rate of non-conversion during the bisulfite treatment step; by this measure, all libraries had a bisulfite conversion efficiency of at least 98.9%.

Bootstrap sampling to determine margin of error in low-coverage WGBS-seq

In order to determine the number of CG methylation calls required to accurately predict genome-wide methylation levels, bootstrap sampling of reads from the deep-sequenced male elephant shark dataset was performed to generate regular intervals of CG calls from approximately 100 to 30,000. These reads were trimmed, mapped and methylation quantified as described above, and following 1000 iterations, the proportion of data falling within the 0.5-99.5 percentiles was calculated to generate a 99% confidence interval. An asymptotic model described by the equation $y=2.208/\sqrt{x}$ was used to fit a curve to the data. At our minimum sequencing depth of 15,000 CG calls, bootstrap sampling predicts a margin of error (99% confidence interval) of approximately ±1.8 methylation percentage points.

Analysis of deep-sequenced WGBS-seq datasets

CG methylation calls were imported into the SeqMonk program (v1.37.1) for analysis. For elephant shark and sea squirt, custom SeqMonk genomes were built using GFF annotation files downloaded from NCBI and Ensembl, respectively.

To analyse methylation at the level of individual CG dinucleotides, we generated an annotation track of each CG site using Bowtie v1.1.2⁴⁷. A minimum of five methylation calls was required for inclusion of a CG site in analyses.

For mouse, zebrafish and elephant shark, precompiled annotation tracks of repetitive elements generated using the RepeatMasker program were downloaded from UCSC. For sea squirt, we generated these annotations by running the RepeatMasker program (v4.0.6) on the KH assembly with the -s option and specifying Ciona intestinalis as the species. The various classes of transposable elements were extracted from these annotation files and where indicated, merged for analysis. A minimum of five calls was applied as a threshold for inclusion when quantifying individual elements.

To examine methylation profiles across genes or TEs and neighbouring sequences, methylation was quantified at individual CGs and the mean plotted across a size-standardised gene or TE as well as 10kb upstream and downstream regions, using the quantitation trend plot function. Figures were produced using Prism (GraphPad, v7), with smoothing applied to flanking regions by averaging 100 neighbours.

Transcription start sites were defined as 200bp centred on the first nucleotide of an annotated mRNA, and a minimum of five methylation calls was applied as a threshold for inclusion in analyses. For analysis of gene bodies, 2kb running windows were quantified (with a minimum of 50 methylation calls applied for inclusion) within annotated mRNAs, excluding 1kb at the 5’ end, and the mean was reported for each mRNA.

Violin plots and histograms were drawn using the ggplot2 package⁴⁸ in R.

Bioinformatic processing of RNA-seq datasets

We downloaded raw sequencing data from previous studies as follows; sea squirt muscle¹⁴, GEO accession GSE19824, sample GSM497252; elephant shark liver⁴⁹, SRA study SRP013772, run SRR513760; zebrafish muscle⁸, SRA study SRP020008, run SRR800045; mouse liver (ENCODE Consortium^50,51), GEO accession GSE78583, sample GSM2072415.

Trimming was performed to remove both poor-quality calls and adapters sequences using TrimGalore (v0.4.0, default parameters). In addition, 12bp were removed from the 5’ end of sea squirt reads and 10bp from the 5’ end of both elephant shark and mouse reads to avoid sequence biases.

Trimmed reads were aligned to the reference genomes described above with HISAT2⁵² (v2.0.5) using single-end or paired-end mode, as appropriate. Known splice sites were specified from a file built from GTF annotation files downloaded from Ensembl (release 87) using the HISAT2 python script. No GTF file was available for elephant shark, so a GFF annotation file downloaded from NCBI was first converted to GTF format using the gffread program (https://github.com/gpertea/gffread).

Analysis of RNA-seq datasets

Alignments from HISAT2 were imported into the SeqMonk program, specifying a minimum mapping quality of 60 to select only uniquely aligned reads.

The RNA-seq quantitation pipeline was used to generate raw read counts across the exons of nuclear protein-coding genes with a correction for any DNA contamination. Counts were corrected by transcript length and genes were divided into quintiles according to expression level.

Genome-wide methylation profiles of the elephant shark, Callorhinchus milii

To generate genome-wide methylation profiles, we extracted DNA from the liver tissue of one female and one male adult elephant shark and performed whole-genome bisulfite sequencing (WGBS-seq). Detailed sequencing results are provided in Table S1.

As described in the somatic tissues of other vertebrates, we found that methylation is much more prevalent in nuclear DNA at CG dinucleotides (69 – 71.6%) than in non-CG context (0.8 – 1%) or mitochondrial DNA (1.6 – 2.5%; Figure 2A). Low-coverage WGBS-seq demonstrated similar global methylation levels in three additional individuals for liver, and in spleen and pancreas samples (Figure 2B). While we observed a small trend for lower methylation in female samples (Figure 2B; female mean 66.4%, male mean 68.6%), this was not significant according to a t-test (p=0.2308) and within the margin of error expected at this sequencing depth (Figure S1).

Figure 2. Global methylation levels in elephant shark somatic tissues.

A: Global methylation levels of deep-sequenced liver samples in different contexts. ‘CG’ refers to symmetrical CG dinucleotides; ‘Non-CG’ indicates all other sequence contexts. B: Global CG methylation levels in elephant shark tissues examined by low-coverage sequencing. The horizontal bar indicates the mean; gold dots, female samples; blue dots, male samples. The difference between female and male liver samples is not significant according to a t-test, and within the technical margin of error expected at the threshold sequencing depth used (±1.8 methylation percentage points; Figure S1).

We proceeded with further analysis of CG methylation in deep-sequenced liver datasets as an example of the elephant shark somatic methylome, and combined male and female samples to enhance sequencing coverage.

The elephant shark genome is pervasively methylated

Existing data indicate that methylation patterns differ markedly between vertebrates and invertebrates. In order to delineate the characteristics of these disparate systems and establish their relationship to the elephant shark methylome, we reanalysed published WGBS-seq data from two vertebrates, mouse (Mus musculus)⁵ and zebrafish (Danio rerio)⁸, as well as an invertebrate from the closest chordate outgroup, the sea squirt Ciona intestinalis¹⁴ (Table 1A).

As expected from analysis of global levels, examination of methylation at individual CG dinucleotides in the elephant shark showed that the majority of sites are highly methylated (≥ 80%), and fewer than one tenth are unmethylated (Figure 3A). Both this pattern and the global methylation level are comparable to mouse and zebrafish (Figure 3A–B). In contrast, mean methylation in the invertebrate sea squirt is only 22.9%, and over two thirds of CG sites are unmethylated.

Figure 3. Global structure of the elephant shark methylome.

A: Distribution of methylation at individual CG dinucleotides. ‘M’ denotes percentage CG methylation. B: Mean methylation of CG dinucleotides. C: Distribution of methylation within 2kb running windows covering the entire genome. Black dots denote the median. D: Genome screenshots of methylation quantified in 2kb running windows over the first 3Mb of chromosome 1 in sea squirt, zebrafish and mouse, and of the largest scaffold (NW_006890054.1) in elephant shark. These regions were arbitrarily chosen as an unbiased section of each genome.

A further striking distinction is evident when the genome is profiled in 2kb running windows. High methylation levels are pervasive in the elephant shark genome (Figure 3C–D), resembling the structure of other vertebrate methylomes. In contrast, the sea squirt methylome is characterised by a bimodal but largely unmethylated distribution (Figure 3C), resulting from a mosaic pattern in which background hypomethylation is punctuated by shorter stretches of methylated sequences (Figure 3D). Interestingly, running windows show a broader distribution of methylation in elephant shark than in mouse or zebrafish (Figure 3C). Whether this is a feature of basal vertebrates generally or of elephant shark specifically will require analysis of methylation patterns in additional cartilaginous fish.

Transposable elements are hypermethylated in the elephant shark

Having established that the global structure of the elephant shark methylome is characteristic of vertebrates, we sought to determine the profile and impact of methylation at specific functional elements.

Transposable elements (TEs) are highly methylated in vertebrate genomes, a feature which is linked to the necessity of repressing their transcription to prevent destabilising transposase activity^{8,9,16,22–24}. The generally low levels of methylation at TEs in invertebrates such as the sea squirt do not appear to regulate their activity^14,38,39.

Examination of methylation patterns at TEs and flanking sequences showed that the elephant shark exhibits hypermethylation at the large majority of TEs and a slight increase in mean methylation relative to adjacent regions (Figure 4A–B), conforming to the pattern of other vertebrates. While mean methylation levels of TEs in sea squirt are moderately elevated compared to flanking sequences, the large majority of TEs are hypomethylated. Little variation in methylation was observed between the two predominant TE classes in the elephant shark genome, long interspersed nuclear elements (LINEs) and short interspersed nuclear elements (SINEs; Figure 4C–D), indicating that – as in other vertebrates^8,9,16,24 – hypermethylation of TEs is ubiquitous.

Figure 4. Methylation patterns at transposable elements.

A: Distribution of methylation at transposable elements. Mean methylation values are divided into 10 bins. B: Mean CG methylation across transposable elements and 10kb flanking regions. Quantification was performed at the level of individual CG dinucleotides. Flanking regions were smoothed by averaging 100 neighbours. C: Distribution of methylation at long interspersed nuclear elements (LINEs) and short interspersed nuclear elements (SINEs) in the elephant shark genome. Mean methylation values are divided into 10 bins. D: Methylation at long interspersed nuclear elements (LINEs) and short interspersed nuclear elements (SINEs) in the elephant shark genome, plotted as in (B).

Methylation at elephant shark transcription start sites is associated with gene silencing

Silencing of gene expression through the deposition of methylation at transcription start sites (TSSs) constitutes an important regulatory mechanism in vertebrates, but appears to be absent from invertebrates^{5–10,12–16,35,40,41}. To compare the relationship of methylation and transcription in elephant shark with higher vertebrates and the sea squirt, we made use of tissue-matched published RNA-seq datasets^8,14,49,50 (Table 1B) to classify protein-coding genes into expression quintiles.

Hypomethylation at the TSS of expressed genes constitutes a conspicuous exception to the otherwise pervasively methylated elephant shark genome, matching the higher vertebrates examined (Figure 5A–C). Significantly, we document an inverse relationship between TSS methylation and expression level in the elephant shark (Figure 5A, Figure 5E). A bimodal distribution in which a large proportion of sequences are methylated at low expression levels contrasts with negligible methylation at most TSSs of intermediate and highly expressed genes. The association of TSS methylation with transcriptional silencing is a distinguishing feature of higher vertebrate methylomes^5,8,12–16 that is recapitulated here for zebrafish and mouse (Figure 5B–C, Figure 5E), and its presence in the elephant shark indicates that methylation at the TSS induces repression in a similar manner. Consistent with reports showing that invertebrates lack this wide variation in TSS methylation as a function of expression level^14,35,41, the large majority of sea squirt TSSs are hypomethylated at all expression levels and methylation levels at the TSS are comparable to intergenic sequences (Figure 5D–E).

Figure 5. Relationship between methylation and gene expression.

A – D: Mean CG methylation across genes and 10kb flanking regions, classified into quintiles according to expression level in RNA-seq datasets (5 = highest). Quintile 4 is omitted for clarity. Quantification was performed at the level of individual CG dinucleotides. Flanking regions were smoothed by averaging 100 neighbours. E – F: Distribution of methylation at the transcription start site (E) and within the body (F) of genes classified into quintiles according expression level (5 = highest). Each violin is scaled to the same maximum width (total area is not constant between violins) to demonstrate distributions for each quintile. Black dots denote the median.

Interestingly, a larger number of TSSs at highly expressed genes remain methylated in elephant shark compared to mouse and zebrafish. This may suggest that the association of methylation with repression is less absolute than in higher vertebrates, but could also be attributed to poorer TSS annotation in the less intensively-studied and incompletely assembled elephant shark genome.

The methylomes of higher vertebrates and invertebrates also differ within gene bodies. While intragenic methylation in sea squirt forms the bimodal distribution reported in invertebrates^35,37, and most silenced genes lack methylation, vertebrate gene bodies are generally hypermethylated at all expression levels (Figure 5F). Intragenic methylation in the elephant shark is characteristic of this vertebrate pattern. In addition, higher expression levels are associated with moderately elevated gene body methylation in elephant shark liver, but not in zebrafish muscle or mouse liver. Given the limited understanding of the role played by intragenic methylation in the regulation of vertebrate gene expression, the functional relevance of this relationship is unclear.

Evolutionary history of the vertebrate methylation system

The observation that methylation patterns in a cartilaginous fish are characteristic of higher vertebrates implies the conservation of a complex methylation system across jawed vertebrates separated by 465 million years of evolution (Figure 1). Of particular note, they support the common presence of a regulatory architecture that links methylation at the TSS to transcriptional repression.

Preprint methylome data from the sea lamprey Petromyzon marinus, a basal jawless vertebrate, indicate that this species lacks the genome-wide hypermethylation and functional relationships of higher vertebrates (https://doi.org/10.1101/033233). While the data from this study has not yet been released, the authors state that methylation patterns in sea lamprey more closely resemble those of the sea squirt and appear to represent a transitional intermediate. In the context of our findings, this implies that the evolution of the higher vertebrate methylation system was achieved after the emergence of jawed vertebrates (~600 Mya⁴³), but before the divergence of bony and cartilaginous fish (~465 Mya⁴³; Figure 1). These data further identify cartilaginous fish as the most divergent class to possess a DNA modification system similar to our own, and position the elephant shark as a valuable model to examine the function and evolution of the vertebrate methylation system. As the slowest evolving vertebrate documented⁴⁹, the elephant shark bears the closest resemblance to the most recent common ancestor of all jawed vertebrates, enhancing its appeal in this respect. Moreover, the extensive orthology of its small genome to those of tetrapods⁴⁹ facilitates comparative studies.

Transposon aggressiveness correlates with the degree of sexual outcrossing in the host, and repression of this destabilising activity has been proposed as a major reason for genome-wide hypermethylation in sexually-reproducing organisms such as plants and vertebrates^14,38,54. This control mechanism appears to have been discarded as unnecessary in early asexual metazoans, and alternative suppression systems such as the piwi-piRNA pathway were developed in their sexually-reproducing invertebrate descendants^54,55. The reason for the apparent reinvention of methylation-based silencing in vertebrates is unclear. Comparison of TE dynamics in the cells of elephant shark and basal chordates offers the opportunity to determine whether the need for additional control mechanisms was a primary driver for genome-wide hypermethylation in jawed vertebrates.

We note that in addition to substantial physiological changes, the emergence of jawed vertebrates was accompanied by major innovations in gene regulatory networks, notably non-coding RNA elements⁴⁹. These advances may have facilitated, or conversely been enabled by, the development of a complex methylation system during the same time period. The role of the whole-genome duplications that occurred in vertebrate progenitors⁵⁶ in the acquisition of components that act downstream of the methylation signal, or as a stimulus for new mechanisms of regulating gene dosage, also merits further investigation.

Methylation of elements that modulate gene expression forms an epigenetic memory that plays an important role in defining and stabilising cell identity in higher vertebrates^5,17–21. The reprogramming of this specification in the germline to regenerate full developmental competence after fertilisation, and the pathways employed to achieve this – such as active demethylation by ten-eleven-translocase (TET) enzymes, vary considerably across vertebrates⁵⁷. Examination of these phenomena in the elephant shark will provide insight into the evolutionary history of epigenetic control in the life cycle and its consequences for vertebrate development.

Our findings provide fresh perspective on an important epigenetic modification. The elephant shark methylome delineates the evolutionary extent of the complex methylation system found in higher vertebrates, and sets the scene for comparative studies that will address longstanding questions about the primary purpose of this system and how these functions evolved from the mosaic pattern of invertebrates. It will be particularly pertinent to understand the development of the mechanism that links TSS methylation to transcriptional repression. Epigenetic studies in the elephant shark also open promising avenues to explore the ways in which methylation is put to use during development and the specification of cell fate, and the conservation of these strategies amongst vertebrates.