A dynamic histone-based chromatin regulatory toolkit underpins genome and developmental evolution in an invertebrate clade

doi:10.21203/rs.3.rs-5284446/v1

Download PDF

Research Article

A dynamic histone-based chromatin regulatory toolkit underpins genome and developmental evolution in an invertebrate clade

https://doi.org/10.21203/rs.3.rs-5284446/v1

This work is licensed under a CC BY 4.0 License

You are reading this latest preprint version

Eukaryotic histones protect and package nuclear DNA into nucleosomes. The dynamic addition and removal of posttranslational modifications on these proteins define regulatory regions that play a central role in genome and chromatin biology. However, our understanding of these regulatory mechanisms in animals is primarily based on a few model systems, preventing a general understanding of how histone-based regulation unfolds and promotes phenotypic variation during animal embryogenesis. Here, we apply a comprehensive multi-omics approach to dissect the histone-based regulatory complement in Annelida, one of the largest invertebrate phyla. Annelids exhibit a conserved histone repertoire organised in clusters of dynamically regulated, hyperaccessible chromatin. However, unlike other animals with reduced genomes, the worm Dimorphilus gyrociliatus shows a dramatically streamlined histone repertoire, revealing that genome compaction has lineage-specific effects on histone-based regulation. Notably, the annelid Owenia fusiformis has two H2A.X variants that co-occur in other animals, sometimes associate with fast cell divisions, and represent a unique case of widespread parallel evolution of a histone variant in Eukarya. Histone-modifying enzyme complements are largely conserved amongst annelids. Yet, temporal differences in the expression of a reduced set of histone modifiers correlate with distinct ontogenetic traits and variation in the adult landscapes of histone modifications, as revealed by quantitative mass spectrometry in O. fusiformis and Capitella teleta. Collectively, our unparalleled analysis of histone-based epigenetics within a non-model phylum informs the evolution of histone-based regulation, presenting a framework to explore how this fundamental genome regulatory layer contributes to developmental and morphological diversification in annelids and animals generally.

Genome packaging is a critical phenomenon in cell biology. Different mechanisms, like DNA supercoiling, aid in this process (Gilbert and Allan 2014), yet the packaging in certain viruses (Bryson et al. 2022; Irwin and Richards 2024), bacteria (Hocher et al. 2023), archaea (Pereira et al. 1997), and especially in eukaryotes (Khorasanizadeh 2004; Abouelmagd and Ageely 2013; Nelson et al. 2021), results from the interaction of genomic DNA with histones. Eukaryotic histones are highly basic, DNA-associated proteins that protect and package nuclear, naked DNA to arrange them into higher-order and densely compacted structures, all the way up to the chromosome level. By doing so, they play a central role as regulators of genome architecture and almost all DNA-based biological processes, especially the regulation of gene expression (Khorasanizadeh 2004). Histone-based regulation is, indeed, one of the most versatile and logically intricate mechanisms modulating genome functioning (Lee et al. 2010; Bannister and Kouzarides 2011; Millán-Zambrano et al. 2022) (Fig. 1A). These mechanism include histone dynamics and differential usage of histone variants (Henikoff and Smith 2015; Buschbeck and Hake 2017; Ferrand et al. 2020). Yet, arguably, the most functionally important regulatory mechanism involves the deposition and removal of small covalent posttranslational modifications, such as methyl and acetyl moieties, to key residues primarily located in the N-terminal end, or tail, of histones (Kimura 2013; Skvortsova et al. 2018). Certain histone posttranslational modifications (hPTMs) cause structural perturbations in the nucleosome structure, thereby changing the level of chromatin condensation or accessibility. Others act as binding marks for modules of reader proteins and complexes that recruit chromatin-modifying machinery, leading to different biological readouts (Workman and Kingston 1998; Bannister and Kouzarides 2011; Yun et al. 2011; Tessarz and Kouzarides 2014). In this way, hPTMs define regulatory regions in the genome and their levels of regulatory activity, constituting critical regulators of various processes like proliferation, differentiation, and most importantly, embryogenesis and development (Chi et al. 2010; Audia and Campbell 2016; Jambhekar et al. 2019; Völker-Albert et al. 2020; Millán-Zambrano et al. 2022). Even though the central histone-based regulation mechanisms are largely evolutionarily conserved, their vast regulatory activity also makes them a very promising source of evolutionary variation.

Little is known about the repertoire and role of hPTMs outside traditional animal models for which high-quality reference genomes and methods to interrogate histone regulation have been readily available. In Spiralia, one of the three largest clades of bilaterally symmetrical animals comprising morphologically diverse lineages and nearly half of the animal phyla (Martín-Durán and Marlétaz 2020), functional genomics work has been limited (Martín-Zamora et al. 2023a). Research on histones peaked during the second half of the 20th century. Those early works were primarily limited to biochemical studies on variations in histone sequence, structure and composition across a handful of major groups, namely molluscs (van Helden et al. 1979; Wouters-Tyrou et al. 1982; Sellos 1985; Ausio and Van Holde 1988), annelids (Kmiecik et al. 1983; Sellos et al. 1990; Davis et al. 1992), and nemerteans (Wang and Ausió 2001). Likewise, hPTM-targeting functional genomics studies with ChIP-seq (Johnson et al. 2007; Kasinathan et al. 2014) and more modern techniques (Schmidl et al. 2015; Skene and Henikoff 2017; Kaya-Okur et al. 2019; Meers et al. 2019; Kaya-Okur et al. 2020) have so far been restricted to Platyhelminthes, where efforts have focused on understanding stem cell maintenance and regeneration in the planarian Schmidtea mediterranea (Rink 2013; Duncan et al. 2015; Dattani et al. 2018; Mihaylova et al. 2018; Blanco et al. 2020; Li et al. 2021; Verma et al. 2021; Neiro et al. 2022; Pascual-Carreras et al. 2023) and the life cycle of the parasitic helminth Schistosoma mansoni to develop potential treatments to schistosomiasis, the disease it causes (Cabezas-Cruz et al. 2014; Liu 2016; Picard et al. 2016; Roquis et al. 2016; Anderson et al. 2017; Cosseau et al. 2017; Roquis et al. 2018; Picard et al. 2019; Ghazy et al. 2022). Therefore, our knowledge of histone-based regulation in non-model systems such as spiralians is quite limited, thus hampering a comprehensive understanding of how changes in genome regulation promotes variation in gene expression programmes, and ultimately phenotypic change, during animal evolution.

Variations in the timing or rate of development, commonly called heterochronies, are a key source of evolutionary variation (McNamara 2012; Keyte and Smith 2014; Buendía-Monreal and Gillmor 2018; Dobreva et al. 2022). In annelids, one of the largest spiralian clades, and likely in bilaterians generally (Degnan and Degnan 2023; Formery and Lowe 2023; Martín-Zamora et al. 2023b), developmental heterochronies correlate with the diversification of embryonic outcomes and life cycles. The annelids Owenia fusiformis Delle Chiaje, 1844; Capitella teleta Blake, Grassle & Eckelbarger, 2009; and Dimorphilus gyrociliatus (O. Schmidt, 1857) share a common early developmental programme—spiral cleavage—but have distinct strategies for cell fate specification, life cycles, and larval types (Fig. 1B) (Seaver et al. 2005; Kerbl et al. 2016a; Carrillo-Baltodano et al. 2021). In these annelids, genes involved in key cellular processes (e.g., autophagy) and enzymatic pathways (e.g., chitin synthesis pathway) are either pre- or post-displaced between larval types, i.e., they have an earlier or later expression onset, respectively (Martín-Zamora et al. 2023b). Moreover, the trunk patterning transcriptional programme, which includes the expression of the Hox genes, is progressively brought forward in development—i.e., pre-displaced—as larval traits are lost and indirect development transitions into direct development (Martín-Zamora et al. 2023b). Importantly, these changes in gene expression mirror chromatin accessibility dynamics, as observed in the open chromatin regions of the Hox cluster and the transcription factor binding dynamics of HOX DNA-binding motifs (Martín-Zamora et al. 2023b).Therefore, given their fundamental role in regulating gene expression, differences in the dynamics of histone-based regulation—particularly of hPTMs—may be upstream and account for these heterochronies, providing a mechanistic explanation of how changes in genome regulation trigger adaptive variation in developmental programmes.

Here, we begin to tackle the relationship between histone-based regulation and developmental transcriptional programmes by mining and characterising the histone complement and histone modification machinery in the annelids O. fusiformis, C. teleta, and D. gyrociliatus through a large-scale multi-omics profiling at the genomic, transcriptomic, epigenomic, and proteomic levels. We describe a minimal histone complement in the miniature genome of D. gyrociliatus and identify divergence in a key histone variant throughout Eukarya that may account for the evolution of early embryonic phenotypes in animals. Despite their differences in histone numbers, the repertoires of histone modifiers are largely complete and conserved between these annelids, although some families exhibit expansions and domain fusions that hint at neofunctionalization events. Concomitant to the heterochronies in other gene regulatory networks and pathways, some histone modifiers are shifted in their developmental time of expression between these annelids, correlating with their life cycle and larval type differences, as well as with the levels of hPTMs in the adults. Altogether, our study is the most comprehensive comparative characterisation of histone-based regulation in a spiralian clade, paving the way for the genome-wide profiling of hPTMs in annelids and the functional assessment of the interplay between histone modifications and development programmes in the phenotypic diversification of this major animal group.

Annelids exhibit a conserved histone repertoire

To investigate the histone gene repertoire of spiralians and annelids, we first reannotated the H2A, H2B, H3 and H4 histone gene models by combining a protein-to-genome alignment-based approach, transcriptomic developmental time series data, a careful manual curation, and an orthology assignment based on maximum likelihood and Bayesian reconstructions (Fig. 2A; Supplementary Fig. 1, 2). We recovered 81, 90, and 17 core histone genes for O. fusiformis, C. teleta, and D. gyrociliatus, respectively (Table 1). In these species, only six, three, and four genes, respectively, displayed an atypical divergent amino acid sequence, which we regard as unknown histone variants. Out of these genes with unknown orthology, all six in O. fusiformis and one in C. teleta were deemed pseudogenic due to having either null or negligible low expression levels throughout development (Supplementary Fig. 7E–G, 8C). Unlike in other protostomes, where certain histone variants have been lost, like H2A.X in Nematoda (Lemmens and Tijsterman 2011), or where features and roles from multiple variants get merged into a single gene, like H2A.X and H2A.Z in the H2A.v variant in Drosophila (Baldi and Becker 2013), all three analysed annelid taxa have conserved at least one ortholog of each of the core histone proteins—canonical and variant—presumed to be ancestral to Bilateria (Henikoff and Smith 2015; Grau-Bové et al. 2022), i.e., the canonical H2A, H2B, H3, and H4 proteins, and the H2A.X, H2A.Z (or H2Av), macroH2A, H3.3, and cenH3 histone variants, suggesting a potential functional conservation of the histone complement. Owenia fusiformis and C. teleta possess between 16 and 23 copies per canonical histone (Table 1). However, the compact genome of D. gyrociliatus only encodes two genes per canonical histone (Table 1), thus making this fully conserved histone complement one with the lowest (if not the lowest) described copy numbers for canonical histones in a metazoan lineage.

Table 1

Histone gene repertoire in Annelida.
Protein	O. fusiformis	C. teleta	D. gyrociliatus
H2A family	24 (21)	23	6
H2A canonical	17	20	2
H2A.X	2	1	1
H2A.Z	1	1	1
macro-H2A	1	1	1
H2A unknown variants	3 (0)	0	1
H2B family	21 (18)	19	3
H2B canonical	18	18	2
H2B unknown variants	3 (0)	1	1
H3 family	20	25 (24)	6
H3 canonical	17	22	2
H3.3	1	1	1
cenH3	1	1	1
H3 unknown variants	1	2 (1)	2
H4 family	16	23	2
H4 canonical	16	23	2
Total	81 (75)	90 (89)	17
Total number of canonical and variant histones encoded in the genomes of O. fusiformis, C. teleta, and D. gyrociliatus, classified by family (H2A, H2B, H3, and H4). For each family, the different canonical and variant histones are shown. Numbers in parentheses denote the number of non-pseudogenic genes.

Genomic organisation and gene structure of histones

Canonical histones tend to appear in tandemly repeated organised clusters, sometimes even associated with H1 linker histones (Nagel and Grossbach 2000; Buschbeck and Hake 2017; Amatori et al. 2021). To assess whether this also occurs in annelids, we characterised the genomic organisation of the core histone genes in O. fusiformis and D. gyrociliatus, the two focal species with more contiguous reference genomes (Fig. 2B; Supplementary Fig. 3A). In O. fusiformis, canonical core histones are arranged into two large clusters containing 35 and 29 canonical histones in chromosomes 5 and 8, respectively, and a smaller one containing four genes only in chromosome 7. With the exceptions of the chromosome 7 cluster and a single occurrence in the chromosome 5 cluster, core canonical histones are tandemly arrayed as a H4–H3–H2A–H2B repeating unit. In this unit, H4, H3, and H2B are transcribed from one strand, and H2A is transcribed from the opposite one (Fig. 2B, Supplementary Fig. 3A). Concomitant to the extremely reduced histone copy number, the genomic clusters in D. gyrociliatus are minimal in size, consisting of a single H3–H4–H2B–H2A repeating unit in scaffolds 2 and 29 (Fig. 2B). This may result from a genomic translocation of the H2A–H2B and H3–H4 units, with the transcripts still being synthesised from the same strands as in O. fusiformis. The presence of these H2A–H2B and H3–H4 units makes it plausible that previously described head-to-head gene pairs in C. teleta and Helobdella robusta (Grau-Bové et al. 2022) are also found in O. fusiformis and D. gyrociliatus and potentially in other annelids like Chaetopterus variopedatus (del Gaudio et al. 1998). Interestingly, Owenia fusiformis and D. gyrociliatus display different genomic organisations than those reported for the annelids Platynereis dumerilii (Sellos et al. 1990; Krawetz et al. 1993) and Urechis capo (Davis et al. 1992). Although these are not based on high-quality reference genome analyses and could be reconciled with a single translocation event from the putative O. fusiformis state, they might also indicate that different genomic organisations of the histone complement occurred in annelids.

All clustered histone genes lack introns in both species, whereas histone variants are often multi-exonic (Fig. 2B, Supplementary Fig. 3A). Interestingly, differences in the overall gene length between species are not statistically significant (Fig. 2C; Supplementary Fig. 3B, C). However, the overall intron length of histones, as well as the total and average intron length per gene of intron-containing genes, but not the number of introns per gene, are significantly lower in both C. teleta and D. gyrociliatus compared to O. fusiformis (Fig. 2D–G; Supplementary Fig. 3D–F). This reduction of average intron length can even be traced gene by gene in the single-copy histone variant genes h2ax, h2az, macroh2a, cenh3, and h3–3 (Supplementary Fig. 3G–I). Therefore, annelids organise their histone genes in tandemly repeated organised clusters. Still, significant changes have occurred in both the cluster structure and number, as well as the gene structure, throughout annelid diversification.

Histone genes are in dynamic hyperaccessible chromatin regions

To better understand chromatin regulation around histone genes and their genomic clusters, we leveraged our publicly available ATAC-seq developmental time series for O. fusiformis and C. teleta (Martín-Zamora et al. 2023b), as well as data for D. gyrociliatus female adults (Martín-Durán et al. 2021). Histone genes are amongst those with the highest ATAC-seq enrichment at every developmental point in O. fusiformis and C. teleta, albeit more starkly in the former (Fig. 3A, B, Supplementary Fig. 4A–D). In these species, open chromatin is most prevalent at the transcription start sites and promoter regions, progressively decreasing along the gene body towards the transcription end sites. During the embryogenesis of O. fusiformis and C. teleta, chromatin around histone genes becomes more accessible after the blastula stage during gastrulation, reaching its peak openness and remaining at very elevated levels during larval development (Fig. 3A, B, Supplementary Fig. 4A–D). At the end of the developmental time course, however, chromatin around the histone clusters becomes more compact. This above-average enrichment is present in the clusters that contain canonical histones—especially in the large clusters of O. fusiformis, where it is very significant—and less so in the loci of the variant histone genes (Fig. 3C, D, Supplementary Fig. 4C, D, 5A–C, 6A, B). In D. gyrociliatus, histone genes are also in loci of highly open chromatin (Supplementary Fig. 4E, F). However, their chromatin accessibility landscape differs from those of O. fusiformis and C. teleta, with chromatin being open almost exclusively at the TSS in defined peaks (Supplementary Fig. 5D, 6C). Therefore, histones are in dynamically regulated, hyperaccessible chromatin regions in annelids, suggesting their expression might be more variable than expected for genes essential for DNA compaction and regulation.

Histone gene expression dynamics in Annelida

To investigate histone gene expression dynamics across annelid development, we mapped the developmental RNA-seq time courses for all three species (Burns and Pechenik 2017; Martín-Durán et al. 2021; Martín-Zamora et al. 2023b) against the new gene models containing the curated histone genes. Canonical histone gene expression levels were calculated as the sum of the expression of all the genes encoding the same protein (e.g., all 17 genes for caH2A in O. fusiformis). All four canonical histones (caH2A, caH2B, caH3, and caH4) present similar expression patterns within species (Supplementary Fig. 7A–D, 8A). Moreover, in both species with a larval stage, i.e., O. fusiformis and C. teleta, histones are expressed in the zygote and cleavage stages—particularly in O. fusiformis—but reach maximum expression levels after the zygotic genome activation, at the gastrula stage and immediate post-gastrula stages. Thereafter, they progressively decrease until the competent larva and juvenile stages (Supplementary Fig. 7A–C, 8A). This mirrors and correlates with the developmental dynamics of the accessible chromatin regions in which they are located. In the direct developer D. gyrociliatus, however, maximum expression of canonical histones occurs at the early embryogenesis time point and plummets immediately after (Supplementary Fig. 7A, D, 8A). More interesting are the expression patterns of histone variants, whose expression dynamics across the development of all three studied annelids are conserved, thereby hinting possible conservation of their roles in development (Supplementary Fig. 7A, E–G, Supplementary Fig. 8B). The h2az histone, which is known to maintain the pluripotency state of stem cells in development (Giaimo et al. 2019; Dijkwel and Tremethick 2022), is expressed at high levels in cleavage and early embryogenesis, and it decreases after gastrulation when cell fates are acquired and terminal differentiation into multiple cell lineages starts. Very similar is the pattern of h3–3, which incorporates into chromatin to sustain several differentiation programmes across metazoans (Sakai et al. 2009; Bush et al. 2013; Jang et al. 2015). The centromeric cenh3, which is critical in kinetochore positioning and assembly during cell division (Quénet and Dalal 2012; Kixmoeller et al. 2020), is also highly expressed at cleavage stages when cells are dividing quickly, but its levels dilute earlier and decrease even before the blastula forms, unlike h3–3 and h2az. The only case in which expression increases dramatically throughout development to reach a maximum in the adult is the macroh2a gene, which might be biologically relevant given its known role in transcriptional repression in maintaining the identity of terminally differentiated cells (Buschbeck et al. 2009; Gaspar-Maia et al. 2013). We also found distinct expression patterns of the unknown/unidentified histone variants across all species (Supplementary Fig. 7A, H–J, 8C), particularly that of h3uv1 in O. fusiformis and h2buv1 in C. teleta, which appear to be extremely abundant genes. We systematically explored the conservation of known hPTMs sites in these genes and all other mined histone genes (see Data Availability section), but we still we have no hints of their functional relevance.

Histone H2A.X variants have evolved in parallel in Metazoa

The histone variant H2A.X, critical in DNA repair (Rogakou et al. 1999; Shroff et al. 2004) and embryonic stem cell development (Andäng et al. 2008; Wu et al. 2014), among other non-canonical functions (Turinetto and Giachino 2015), is strikingly divergent across the studied annelids. While other variants have a single ortholog in all three species, C. teleta and D. gyrociliatus display a single h2ax gene, but O. fusiformis encodes two different H2A.X proteins that share an 88.1% sequence identity (Fig. 2A, B; Supplementary Fig. 1, 2, 3A; Table 1). Both paralogs are maternally deposited in the egg and highly abundant in early embryonic stages, and progressively lose expression after the blastula and gastrula stages up to their minima in the juvenile adult stage (Fig. 4A; Supplementary Fig. 8B). The same expression dynamics are evident for C. teleta and D. gyrociliatus (Fig. 4A, Supplementary Fig. 8B). Nevertheless, while developmental dynamics of both paralogs in O. fusiformis are similar, the expression of h2ax2 is up to 150-fold greater than that of h2ax1, confirming that h2ax2 is the dominant H2A.X gene (Fig. 4A, Supplementary Fig. 7E). Yet, despite the lower levels of expression, we could confidently prove h2ax1 is expressed throughout the life cycle of O. fusiformis by specifically amplifying both genes from a cDNA pool of mRNA coming from all developmental time points of O. fusiformis (Supplementary Fig. 9). Remarkably, some of the differences in protein sequence between h2ax1 and h2ax2 lie in critical regions, most strikingly in the C-terminal motif of the protein (Fig. 4B, C; Supplementary Fig. 10A). In mammals and most chordates, Tyr142 is a conserved position that undergoes biologically relevant post-translational modifications (Cook et al. 2009; Xiao et al. 2009). Where h2ax1 displays this phosphorylatable Tyr142 (from here on also H2A.X-Y), h2ax2 has a sterically homologous yet unmodifiable Phe142 (thus being called H2A.X-F) (Fig. 4B). Indeed, this residue change has been observed before in the H2A.X-F gene of Xenopus, where it is expressed in eggs and early embryos and where it has been hypothesised to facilitate rapid early-embryo cell divisions (Shechter et al. 2009; Wang et al. 2014), but also in some fish (Shechter et al. 2009), a mollusc (González-Romero et al. 2012), and in plants (Shechter et al. 2009; Lorković et al. 2017; Lei and Berger 2020). Indeed, we found that some of the few curated spiralian H2A.X sequences—and also some arthropod and other eukaryote ones—display this change (Fig. 4D; Supplementary Fig. 10B–F). Therefore, the distribution of H2A.X-Y and H2A.X-F variants is more widespread than previously foreseen. Yet, it is unknown whether these originated in the last common eukaryotic ancestor or are the result of parallel evolution.

To disentangle the evolutionary history of H2A.X, we performed PHI-BLAST searches for H2A.X-Y and H2A.X-F orthologs across Eukarya combined with phylogenetic reconstruction. We recovered 1,656 H2A.X sequences (660 H2A.X-Y and 996 H2A.X-F), mostly belonging to the highly sequenced Metazoa and Plantae clades, across 29 different eukaryotic phyla (Supplementary Fig. 11B–G). Few clear trends could be elucidated beyond almost all plants displaying an H2A.X-F variant, with the few plant H2A.X-Y variants (nine out of 635) being almost entirely restricted to chlorophytes. In arthropods, the ancestral and most common scenario seems to be an H2A.X-Y from which divergent H2A.X-F proteins evolved. Within Chordata, on the other hand, two distinct clades correspond to each of the two H2A.X variants. In each clade, some conversions to the opposite scenario (Y-to-F and F-to-Y transitions) occur, yet generally, there is strong conservation of the amino acid in position 142. This robustly suggests that in chordates these are different genes with distinct evolutionary origins (Fig. 4E; Supplementary Fig. 11A) but that other evolutionary dynamics are also present in animal and other eukaryote lineages (Fig. 4F). Previous research suggested possessing both a H2A.X-Y and a H2A.X-F variant could be an adaptation for quickly developing aquatic species, yet 36.5% of the species with both variants do not display aquatic lifestyles (Fig. 4G). Here we show, however, that species with this two-H2A.X variant setting are almost exclusively found within Metazoa and are more common than previously thought, as shown in species analysed from the Chordata (23.40%), Cnidaria (16.67%), and Mollusca (30%) (Fig. 4H), likely reflecting H2A.X-Y and H2A.X-F animal variants to be two distinct genes. Altogether, these results demonstrate that histone variants can display intricate phylogenetic histories, highlighting the need for further functional comparative analyses of the genome regulatory implications of the two H2A.X variants in eukaryotes and during animal embryogenesis.

Histone modifiers show a complex evolution throughout annelid diversification

Besides the histone complement, diversity in the repertoire of histone-modifying enzymes (HME) is the other major potential source of evolutionary variation within histone-based regulation. To address this, we investigated the writers and erasers of histone acetylation and methylation: histone deacetylases (HDAC), histone demethylases (HDM), histone methyltransferases (HMT) of both the lysine-specific (KMT) and arginine-specific (PRMT) types, and histone acetyltransferases (HAT), of both type A and type B. We used a mutual best BLAST hit approach to find candidate sequences in the annelid genomes, which were then subjected to orthology assignment using both maximum likelihood and Bayesian reconstructions in six different phylogenetic analyses (Fig. 5; Supplementary Fig. 12–23). We identified a total of 78 clades corresponding to 77 different genes, of which only three (3.9%) were not supported by phylogenetic methods, and only five (6.5%) had bootstrap values or posterior probabilities below 70 or 0.7, respectively, thus making our gene identification highly robust and reliable.

Histone deacetylases

Annelids display at least a hdac1/2, hdac3, hdac4/5/7/9, hdac8, and hdac11 ortholog, as well as one ortholog of each of the known NAD-dependent sirtuins, that is, sirt1, sirt2, sirt3, sirt4, sirt5, sirt6, and sirt7 (Fig. 5A, B; Supplementary Fig. 12, 13). There is an additional well-assigned sirtuin clade, which we have here termed sirtuin novel (sirtn), with a single ortholog in all three species, as well as a clade with high similarity to both sirt6 and sirt7 genes (sirt6/7-like, or sirt6/7l), with two copies present in both O. fusiformis and C. teleta, but not in D. gyrociliatus. The latter presents a duplication of the hdac1/2 gene, whereas the genome of C. teleta encodes for two different hdac6/10 genes. Intriguingly, C. teleta also has an additional HDAC gene that clusters independently, and for which the most similar gene is the poorly characterised hdac12 gene in zebrafish (Fig. 5A, B; Supplementary Fig. 12, 13).

Histone demethylases

In terms of histone demethylases, annelids encode orthologs of kdm1a (lsd1), kdm1b (lsd2), kdm2a/b (fbxl10/11), kdm3 (jmjd1x), kdm4 (jmjd2x), kdm5 (jarid1x), kdm6 (utx/uty), kdm7 (phf2/8), kdm8 (jmjd5), riox1/2 (no66/mina53), jmjd6, and jmjd7 (Fig. 5A, C; Supplementary Fig. 14, 15). All three species encode a single ortholog of these genes, except for kdm5, which is only present in C. teleta, the secondary losses of kdm1b (lsd2) and kdm7 (phf2/8) in D. gyrociliatus, and a lineage-specific duplication of kdm2a/b (fbxl10/11), also in the direct developer D. gyrociliatus. We identified a clade with high similarity to jmjd6, which we have denoted as jmjd6-like (jmjd6l), that is nevertheless distinct from jmjd6 and is present as a single gene in each of the annelids’ genomes. Two other clades phylogenetically closer to jmjd7 yet automatically assigned as kdm8 genes were identified. The analysed genomes contain one copy of each, except in C. teleta, which harbours two genes within the kdm8-like 1 (kdm8l1) clade. Similarly to all other invertebrates, the protein arginine deiminase padi4 is absent in our annelids (Fig. 5A, C; Supplementary Fig. 14, 15).

Histone acetyltransferases

A complete repertoire of nuclear/type A and cytoplasmic/type B occurs in all three annelids (Fig. 5A, D; Supplementary Fig. 16–19). The only exception is kat12 (gtf3c4), which despite being part of a general transcription factor (TFIIIC90) of RNA Pol III (Male et al. 2015; Seifert-Davila et al. 2023), seems to be absent in D. gyrociliatus. At least one ortholog was found for all three lineages for kat2a/b (gcn5/pcaf), kat3a/b (crebbp/ep300), kat4 (taf1), kat5 (tip60), kat6a/b (myst3/4, moz/morf), kat7 (myst2, hbo1), kat8 (myst1, mof), kat9 (elp3), kat13a/b/c (ncoa1/3/2), kat13d (clock), kat14 (csrp2bp), hat1, naa60, and naa40. Unexpectedly, the compact genome of D. gyrociliatus showed the largest expansions, with kat4 (taf1), kat5 (tip60), and kat8 (myst1, mof) duplicated in a lineage-specific manner. But we could also find a likely ancestral duplication common to O. fusiformis and C. teleta of the kat13d (clock) gene (Fig. 5A, D; Supplementary Fig. 16–19).

Lysine-specific histone methyltransferases

The KMT family showed the largest variations between annelids and model organisms. Annelids do not have an ortholog of the chordate-specific kmt3c (smyd2) and kmt3d (smyd1) genes (Calpena et al. 2015), or kmt7 (setd7). They do display orthologs for kmt1c/d (ehmt2/1), kmt1e/f (setdb1/2), kmt2a/b (mll1a/b), kmt2c/d (mll3/2), kmt2e (mll5), kmt2f/g (setd1a/b), kmt2h (ash1l), kmt3a (setd2), kmt3b/f/g (nsd1/3/2), kmt3e (smyd3), kmt4 (dot1l), kmt5a (setd8), kmt5b/c (suv420h1/2), kmt6a/b (ezh2/1), setd6, setmar, and smyd5 (Fig. 5A, E; Supplementary Fig. 20, 21). Strikingly, we could only assign a kmt1a/b (suv39h1/2) ortholog in the case of C. teleta. Dimorphilus gyrociliatus has suffered losses of kmt2h (ash1l), kmt5b/c (suv420h1/2), and setd6, but carries two different genes of the kmt1c/d (ehmt2/1) and kmt2f/g (setd1a/b) classes. Owenia fusiformis has a duplication of the kmt2h (ash1l) gene, whereas C. teleta has one for kmt3e (smyd3). Moreover, we uncovered an expansion of kmt5a (setd8) orthologs exclusive of C. teleta, which has three different paralogs of this gene. Two sequences of the SET and MYND domain-containing family (SMYD)—one belonging to O. fusiformis and another one to D. gyrociliatus—did not cluster with any identified KMT clade and were classified broadly as smyd-like (smydl) genes. Within the kmt8 (prdmx/mecom) clade of genes, both D. gyrociliatus and C. teleta contained three different genes that we could not subclassify but that may represent different genes. Surprisingly, O. fusiformis displayed a large expansion in this family, with up to six different genes that could not be resolved further (Fig. 5A, E; Supplementary Fig. 20, 21).

Arginine-specific histone methyltransferases

Lastly, we also studied the complement of PRMT enzymes (Fig. 5A, F; Supplementary Fig. 22, 23). Excluding prmt2, which could not be found in any polychaete, all others are present, namely, prmt1/8, prmt3, prmt4, prmt5, prmt6, prmt7, and prmt9. A further PRMT gene was identified in all three species that only clustered with the fruit fly Art8 gene and that we deemed art8-like (art8l). The expansions of prmt1/8 and prmt6 in the genome of D. gyrociliatus are thus the only changes within the PRMT family in the studied annelids.

PRMT6 expansions led to domain fusions and likely catalytically dead enzymes

Considering that D. gyrociliatus is a case of extreme genome compaction, frequently associated with gene losses, it is striking that it appeared to possess four different prmt6 genes (hereon referred to as prmt6-a through prmt6-d) that subcluster independently of the rest of the sequences of the other species within the prmt6 clade (Fig. 5A, F; Supplementary Fig. 22, 23). PRMT6 enzymes catalyse type II arginine methylation reactions and yield monomethyl arginine (Rme1/MMA) and asymmetric dimethyl arginine (Rme2a/ADMA) (Fig. 6A) in both histones and non-histone proteins and have been shown to regulate cell proliferation and senescence (Blanc and Richard 2017; Lorton and Shechter 2019). Both the PRMT6-A and PRMT6-C proteins are significantly longer—with a full length of 871 and 1090 amino acids long, respectively—than what is expected of a PRMT6 ortholog (340–380 residues) (Fig. 6B). We analysed their domain and region composition and predicted their protein structure (Fig. 6B–D, Supplementary Fig. 24, 25) and uncovered that while O. fusiformis and C. teleta have a very conserved protein structure and domain length, this is only true for prmt6-b in D. gyrociliatus. Meanwhile, both prmt6-c and prmt6-d have a shorter S-adenosyl methionine (SAM)-dependent methyltransferase class I domain, and prmt6-a and prmt6-c contain additional regions and domains in the N-terminal region of the protein, namely three consecutive galactose/rhamnose-binding lectin domains and a progesterone-induced blocking factor 1 family region, respectively (Fig. 6B, D, Supplementary Fig. 25A). Transcriptomic validation of the fusion genes showed that the one present in prmt6-c is, however, likely a false positive and a result of a wrongful annotation. Unlike in prmt6-a, where both presupposed fused parts of the gene show continuous transcription at similar levels, the prm6-c gene showed large discrepancies in read density across the model and displayed a large number of unexplained antisense reads (Supplementary Fig. 25B, C). No PRMT6 ortholog contains the functionally critical tyrosine dyad (Y47 and Y51) in the PRMT characteristic motif (Fig. 6C, Supplementary Fig. 24). There are also a very high number of amino acid changes in highly conserved positions known to interact with either the SAM cofactor (e.g., R66, E141, and S169) and/or with the arginine residue (e.g., E164 and E317) (Fig. 6C, Supplementary Fig. 24). In fact, structural changes affecting the characteristic PRMT domain are evident, most notably in the lack of tertiary structures of the alpha helices where the PRMT characteristic motif and the tyrosine dyad are located (Fig. 6D, Supplementary Fig. 25A). These sequence and structure alterations are likely disrupting the potentially PRMT6 enzymatic function, thus perhaps relaxing selection pressures and allowing for various gene expansions and divergence events to occur. Altogether, our data indicate that in an annelid with a highly compact genome and a reduced and simplified histone complement, like D. gyrociliatus, the plasticity in histone-based regulation required during embryogenesis most likely resides in the conservation and, in some cases, diversification of the HME repertoire.

Life cycle correlates with heterochronies of histone modifiers

Histone-modifying enzymes (HMEs) constitute a highly plastic source of regulative variation that underpins changes in gene expression patterns (Bannister and Kouzarides 2011; Zhou et al. 2011; Butler et al. 2012; Zhao and Shilatifard 2019; Morgan and Shilatifard 2020). To investigate how temporal shifts in the deployment of orthologous developmental programmes could correlate with large changes in HME expression, we leveraged the entire developmental RNA-seq time series, from the zygote to competent larva/adult stages, for all three annelid species. All expressed transcripts were clustered using an unbiased soft k-means approach into an optimal number of 15 clusters (O. fusiformis and C. teleta) or eight clusters (D. gyrociliatus) (see Methods, Supplementary Fig. 26A–C). Next, we assigned the clusters to their corresponding developmental stages according to their maxima and classified them as cleavage or post-cleavage clusters (O. fusiformis and C. teleta) or early development and late development clusters (D. gyrociliatus). We further subdivided post-cleavage clusters in the indirect developers into pre-larval and post-larval clusters (Supplementary Fig. 26D–F). We profiled the expression dynamics of all five superfamilies of histone modifiers described above across the development of the three species. Given the broad array of biological readouts different modifiers from the same superfamily can have, as expected no clear superfamily-specific trends were common to all three taxa (Supplementary Fig. 27, 28). We then decided to focus on gene-wise expression dynamics. To do this, we compared the timing of deployment of every histone modifier with single-copy orthology either between both indirect-developing species or common to all three species (Fig. 7; Supplementary Fig. 29, 30). When comparing O. fusiformis and C. teleta, we unravelled how there are more histone modifiers under heterochronic shifts between cleavage and post-cleavage clusters (29 genes, 44.6% of the total) than sharing the same expression pattern (27 genes, 41.5%). This is mostly due to 21 enzymes (32.3%) being expressed early during cleavage stages in O. fusiformis but late in post-cleavage time points in C. teleta; yet the opposite is also true for eight of the enzymes (12.3%) (Fig. 7A, top). Heterochronic genes belong to all five superfamilies of modifiers. In the first group of genes (group I), or those delayed in the lecithotrophic larva, we can find, among others, the H3K27ac acetyltransferase kat3a/b (crebbp/ep300) and one of the circadian clock gene orthologs kat13d-a (clock-a); four different sirtuin genes, namely sirt3, sirt4, sirt5, and sirt6; both lysine-specific demethylases kdm1a (lsd1) and kdm1b (lsd2); the H3K79-specific kmt4 (dot1l) and H3K27-specific kat6a/b (ezh2/1); as well as prmt4 and prmt5. Displaced towards a late expression in the planktotrophic larva (group II), examples include the ortholog to the steroid receptor co-activators kat13a/b/c (ncoa1/3/2), hdac11 and sirt1, the bifunctional histone demethylase/ribosome hydroxylase riox1/2 (no66/mina53), the inactive kmt2e (mll5), and the previously discussed prmt6 gene (Fig. 7A, top, D; Supplementary Fig. 29). Contrary to our null expectation of finding heterochronic shifts between pre-larval and post-larval expression, all genes with a post-cleavage expression were consistently deployed at the same life cycle time point—either before or after larval formation—between C. teleta and O. fusiformis (Fig. 7A, bottom), suggesting that histone-based regulatory mechanisms in later developmental stages might be broadly conserved between indirect developing annelids.

We thus focused on the comparison between direct and indirect developers, i.e., between D. gyrociliatus and either C. teleta or O. fusiformis. As observed with transcription factors (Martín-Zamora et al. 2023b), a significant number of HMEs are pre-displaced from late/post-cleavage expression in C. teleta (group III, 14 genes, 29.8%) and O. fusiformis (group IV, 11 genes, 23.4%) to early development in D. gyrociliatus (Fig. 7B, C). Genes involved in these shifts are varied, with all superfamilies except the PRMT enzymes being represented in both groups of heterochronic genes (Fig. 7E, F; Supplementary Fig. 30). A critical advantage of these pair-wise comparisons is that we can intersect them and find those changes that are consistent and specific to direct and indirect developers (Fig. 7G). By comparing groups I. and III., those where genes are shifted from early/cleavage expression in either O. fusiformis or D. gyrociliatus to late/post-cleavage expression in C. teleta, we can robustly describe modifiers specific to the development of the lecithotrophic larva of C. teleta (Fig. 7D, top, E, G, left). These are kat3a/b (crebbp/ep300), kat14 (csrp2bp), hdac8, jmjd6l, kdm1a (lsd1), kdm3 (jmjd1x), kmt4 (dot1l), kmt6a/b (ezh2/1), and smyd5. The equivalent can be produced for O. fusiformis by comparing groups II and IV, which yields the modifiers specific to the development of the planktotrophic larva of O. fusiformis (Fig. 7D, bottom, F, G, centre), namely kat6a/b (myst3/4, moz/morf), kat7 (myst2, hbo1), hdac11, sirt1, and kmt2e (mll5). Lastly, we compared groups III and IV, i.e., those genes that are consistently pre-displaced to early expression in D. gyrociliatus from late/post-cleavage expression in both indirect developers (Fig. 7E, F, G, right), which included kat2a/b (gcn5/pcaf), three different demethylases, jmjd6, kdm4 (jmjd2x), and kdm6 (utx/uty), and kmt2c/d (ehmt2/1). Therefore, a few orthologous HMEs exhibit different temporal expression dynamics between the studied annelids. These might represent key epigenetic regulators of specific developmental programmes that differ between direct and indirect annelids, thus emerging as a tractable gene set for functional manipulations in the future.

Adult annelids harbour distinct histone modification enrichments

Finally, given the differences in the HME repertoires between species and their expression timing across development, we investigated whether adults may consequently display different levels of key histone modifications. To assess this, we used quantitative mass spectrometry-based detection of hPTMs in adult specimens of O. fusiformis and C. teleta and focused on the methylation and acetylation of residues in H3 and H4 histones (Fig. 8A, B; Supplementary Fig. 31A–C). The three biological replicates with two technical replicates strongly clustered by species (Fig. 8C and Supplementary Fig. 31D). Significant differences in key hPTMs were immediately obvious in several peptides. Most notable were the large differences in H3K56 methylation and H3K79 methylation. H3K56 is mostly unmethylated (H3K56me0) or di- and trimethylated (H3K56me2 and H3K56me3) in C. teleta but more abundant in its monomethylated form (H3K56me1) in O. fusiformis (Fig. 8D, left). Notably, these differences are consistent with the differences in the presence/absence of HMEs between species and their developmental heterochronies. H3K56me3 is known to be introduced in humans by the KMT1A (SUV39H1) and KMT1B (SUV39H2) enzymes. The almost exclusive presence of H3K56me3 in C. teleta matches our predictions, given that the annelid KMT1A/B (SUV39H1/2) ortholog is present in C. teleta but absent in O. fusiformis. On the contrary, the ortholog to the human H3K56me1-depositing enzyme KMT1D (EHMT2), which is KMT1C/D (EHMT2/1) in annelids, is present in O. fusiformis, likely explaining the accumulation of H3K56me1 in this species. These different hPTM abundances may have deep biological consequences, as it has been described that while H3K56me1 is involved in regulating DNA replication (Yu et al. 2012), H3K56me3 is a heterochromatic mark that largely overlaps with the constitutive heterochromatin mark H3K9me3 (Jack et al. 2013; Colmenares et al. 2017). We also found large differences in unmodified and monomethylated H3K79, with the former being the dominant form in C. teleta and the latter in O. fusiformis (Fig. 8D, right). Some well-studied and conserved methylation hPTMs like H3K4me2 and H3K4me3—among many others—appeared to be more enriched in one species from EpiProfile analyses, though statistical quantifications showed no significant differences between species (Fig. 8; Supplementary Fig. 31E). Other methyl marks like H3K36me1/2/3 and H3K27me2/3 were significantly more abundant in one of the annelid lineages. We also found this to be true for two acetylations, namely H3K23ac and H3K27ac, both significantly more abundant in C. teleta (Fig. 8F; Supplementary Fig. 31F–H), which could be the result of the heterochronic shift between early/cleavage expression in O. fusiformis and late/post-cleavage expression in C. teleta of five different HAT genes: kat5 (tip60), kat9 (elp3), kat13d-a (clock-a), kat14 (csrp2bp), but most interestingly, kat3a/b (crebbp/ep300). H3K9ac, however, escapes this pattern and appears more abundant in O. fusiformis (Supplementary Fig. 31F). Likewise, this may be associated with chromatin and transcription differences in the annelids, as histone acetylations like H3K23ac and H3K27ac are known to be tightly associated with increased transcription across species (Creyghton et al. 2010; Martin et al. 2021; Zhang et al. 2023). Altogether, our data describe a core repertoire of annelid histone methylation and acetylation of histones H3 and H4, revealing some key differences in the adult levels in the abundance of key hPTM, likely to be the result of differences in the repertoire and expression dynamics of HMEs.

In this study, we profiled the histone, HME, and hPTM repertoires to trace the evolution and developmental dynamics of histone-based regulation in three annelid lineages with contrasting genomic features and ontogenetic trajectories (Fig. 1B). Through a multi-omics approach combining comparative genomics, public bulk RNA-seq and ATAC-seq datasets, and the generation of mass spectrometry proteomic data, we identified multiple layers of histone-based regulation common and specific to these annelids that, in some cases, correlate with their life cycles and developmental modes, raising testable and tractable hypotheses to further explore how this genome regulatory mechanism underpins animal embryonic and life cycle evolution.

The in-depth characterisation of histones, histone modifiers, and histone modifications revealed largely stable and conserved repertoires in our three focal annelids, vastly differing from the large differences we recently described in the DNA methylation levels and associated machinery for these same species (Guynes et al. 2024). The genomic organisation of histone genes of O. fusiformis and D. gyrociliatus determined in this study revealed that H3 and H4 genes are contiguous and always transcribed from the same strand. Yet, when compared with other annelids and molluscs (del Gaudio et al. 1998; Eirín-López et al. 2002; Albig et al. 2003), at least two or more chromosomal rearrangements are required to explain the genomic organization differences, suggesting that either Mollusca or Annelida (or both) have experienced multiple changes in the histone loci from their last common ancestor. Remarkably, D. gyrociliatus displays the lowest characterised histone copy number in a free-living animal but a largely complete (and expanded in some instances) HME repertoire. While we cannot exclude the latter to be explained by the multitude of non-histone targets from HMEs, paradoxically, though, this setting is opposite to Oikopleura dioica, another animal lineage with an extremely compact genome that has, however, experienced rapid innovation in the histone repertoire (Moosmann et al. 2011). Therefore, genome reduction impacts histone evolution in different ways, putting D. gyrociliatus forth as an ideal animal to study the functional impact of hPTMs in animals, as it would simplify the more complex functional investigations implemented in other organisms with more redundant histone repertoires, such as D. melanogaster (Günesdogan et al. 2010; McKay et al. 2015) and, most recently, mice (Sankar et al. 2022).

While previous work recapitulating molecular heterochrony invoked several epigenetic regulators, most notably microRNAs (miRNAs) (Slack and Ruvkun 1997; Moss 2007; Zou et al. 2013; Ouchi et al. 2014; Armakola and Ruvkun 2019), less is known about the role of the histone-based regulation machinery in this fundamental developmental and evolutionary process. Through comparative gene expression analyses of HMEs, we discovered that the expression dynamics of HMEs are not always conserved between annelids. Instead, many of them—in some cases more than 50% of the studied genes—are pre- or post-displaced between larval phenotypes and between species possessing or lacking a larval stage. Very interestingly, HDMs were the most popular category of HMEs consistently shifted between late expression in both indirect developers and early expression in D. gyrociliatus, the same heterochrony type previously described for the trunk patterning programme (Martín-Zamora et al. 2023b). These included kdm4 (jmjd2x), kdm6 (utx/uty) and jmjd6, whose known roles in mammals (Greenfield et al. 1998; Whetstine et al. 2006; Chang et al. 2007; Shpargel et al. 2012; Liu et al. 2013) might indicate a potential removal of methylation marks of both permissive (H3K36me3) and repressive nature (H3K27me3 and H3K9me3) earlier in direct development to accelerate the expression and subduing of pre-displaced transcriptional programmes. Altogether, our work sets the foundations for the study of histone-based regulation in annelids and other non-model systems, defining a set of tractable HMEs and hPTMs, whose detailed study, from spatial expression and regulatory dynamics to, most importantly, loss-of-function studies, will help understand how variation in chromatin regulation contribute to the fascinating phenotypic evolution found in Spiralia and animals generally.

Adult culture

Sexually mature Owenia fusiformis Delle Chiaje, 1844 adults were obtained from subtidal waters near the Station Biologique de Roscoff (Centre National de la Recherche Scientifique (CNRS)-Sorbonne University, Roscoff, France) and cultured in-house as described before (Carrillo-Baltodano et al. 2021). Capitella teleta Blake, Grassle & Eckelbarger, 2009 and Dimorphilus gyrociliatus (O. Schmidt, 1857) were cultured following previously described protocols (Seaver et al. 2005; Kerbl et al. 2016b).

Histone genes mining

Human core histone protein sequences, i.e., H2A, H2B, H3, and H4 proteins, were retrieved from NCBI following the gene definitions provided in the HUGO Gene Nomenclature Committee (HGNC). Curated histone variants sequences from model organisms corresponding to the H2A.X, H2A.Z/H2Av, H3.3, cenH3/CENP-A, and macroH2A subfamilies were downloaded from the NCBI Histone Variants Database 2.0 (HistoneDB 2.0) (https://www.ncbi.nlm.nih.gov/research/histonedb) (Draizen et al. 2016) (Supplementary Table 1). The protein2genome model of the exonerate tool (EMBL-EBI) was employed to align all fetched sequences to either the unmasked or softmasked scaffold-level genome assemblies of O. fusiformis (Martín-Zamora et al. 2023b), C. teleta (Simakov et al. 2013)d gyrociliatus (Martín-Durán et al. 2021). The resulting gff files containing putative orthologs were formatted to gff3 format using custom code. The AGAT v.0.6.0. suite of scripts (Dainat et al. 2021) and the gffutils v.0.8.4 package (https://github.com/daler/gffutils) were used to merge transcripts with overlapping codifying sequences, remove identical transcripts, and keep the longest isoform of gene models containing monoexonic transcripts only. Resulting putative sequences were aligned by histone family using MAFFT v.7. (Katoh et al. 2019) with the G-INS-i iterative refinement method, BLOSUM62 as a scoring matrix and a 1.53 gap penalty. Multiple sequence alignments (MSAs) were used as support, together with visual inspections of the corresponding genomic loci in IGV v.2.8.13. (Robinson et al. 2011), to manually curate histone sequences. Transcripts around or in poorly sequenced regions of the genomes and those with non-canonical splicing sites were discarded. Sequences with frameshifts or in-frame termination codons were also discarded. The sequence of multiexonic transcripts was inferred from full-length RNA-seq transcripts from previously published data for all three species (Burns and Pechenik 2017; Martín-Durán et al. 2021; Martín-Zamora et al. 2023b) found through a tblastn search (Camacho et al. 2009) of partial alignments. Only a single transcript per gene and locus was kept, as alternative splicing leading to different protein sequences is rare in histone genes (Bönisch et al. 2012; Martire and Banaszynski 2020) and could not be inferred with certainty from a genome alignment-based approach. Gene models were then lifted off from the scaffold-level assembly of O. fusiformis (GenBank: GCA_903813345.1) to the updated chromosome-level assembly (GenBank: GCA_903813345.2) using Liftoff v.1.6.1 (Shumate and Salzberg 2021). Curated protein sequences were then extracted using gffread v.0.12.1 (Pertea and Pertea 2020). Alignments and the corresponding alignment files in .fasta and .nexus formats, updated filtered gene models, transcript and protein files, and genome annotation reports for all three annelid taxa are available in the GitHub repository (see Data Availability section).

Histone genes orthology assignment

Curated core histone and histone variant sequences from model organisms were retrieved from the NCBI (Supplementary Table 2) and aligned with the curated annelid histones. Sequences were aligned using MAFFT v.7. (Katoh et al. 2019) as described above. To validate a preliminary sequence-based orthology assignment, a phylogenetic analysis was performed on the resulting MSA with a rtREV amino acid substitution matrix (Dimmic et al. 2002) to account for transition rates. We allowed for a proportion of invariable sites (+ I) and used the discrete gamma model (Yang 1994) with four rate categories (G4) to describe site evolution rates, together with an optimisation of amino acid frequencies using maximum likelihood (ML) in IQ-TREE v.2.0.3 (Nguyen et al. 2015). 1,000 ultrafast bootstraps (BS) (Hoang et al. 2018) were used to extract branch support values, but bootstrapping convergence was reached earlier. Further support in the form of posterior probabilities was obtained from Bayesian reconstruction in MrBayes v.3.2.7a (Ronquist and Huelsenbeck 2003), using the same rtREV + F + I + G4 model. Two runs with four chains of Markov chain Monte Carlo (MCMC) analyses were run for 50,000,000 generations. Tree visualisation and editing were done using FigTree v.1.4.4 (http://tree.bio.ed.ac.uk/software/figtree/). The resulting histone gene orthology assignment is summarised in Supplementary Table 3 and described in detail for each species in Supplementary Tables 4–6. Alignments and the corresponding alignment files in .fasta and .nexus formats are available in the GitHub repository (see Data Availability section).

Histone gene structure analysis

Histone gene structure was analysed using the GenomicFeatures v.1.46.5 and GenomicRanges v.1.46.1 packages (Lawrence et al. 2013). Gene length, aggregated intron length, total and average intron length/gene, and total and average intron number/gene were analysed. For gene structure comparisons between the three species of the histone variants with known role and single-copy orthology, i.e., H2A.Z, H2A.X, macroH2A, cenH3, and H3.3, repeated measures one-way ANOVAs were performed, followed by two-tailed paired Student’s t-test, when applicable. For O. fusiformis, we considered the h2ax2 gene (H2A.X-F variant) the H2A.X ortholog, as it is the paralog with the highest expression levels in every developmental stage. For all other comparisons, we employed one-way ANOVAs, followed by two-tailed post-hoc Tukey tests for pair-wise comparisons, when applicable. All P values derived from pair-wise comparisons were adjusted using the stringent Bonferroni method for multiple testing correction (Supplementary Table 7).

Chromatin accessibility profiling

Chromatin accessibility dynamics were analysed using publicly available coverage .bw files corresponding to ATAC-seq experiments for all three species (Martín-Durán et al. 2021; Martín-Zamora et al. 2023b) (Supplementary Table 8). Metagene enrichment analysis was computed using computeMatrix and plotHeatmap commands in deepTools v.3.4.3 (Ramírez et al. 2016). Visualization of peak tracks and gene structures was conducted using pyGenomeTracks v.2.1 (Lopez-Delisle et al. 2020) and deepTools v.3.4.3 (Ramírez et al. 2016).

Gene expression profiling

Gene expression dynamics were re-profiled using publicly available developmental transcriptomes for all three species (Burns and Pechenik 2017; Martín-Durán et al. 2021; Martín-Zamora et al. 2023b) (Supplementary Table 9) to include the previously incomplete histone gene models, from zygote to juvenile stages for O. fusiformis, from zygote to competent larva stages for C. teleta, and from early embryogenesis to female adult for D. gyrociliatus, as described before. Trimmomatic v.0.3980 (Bolger et al. 2014) was used to remove sequencing adaptors. Clean reads were pseudo-aligned to the updated filtered gene models using kallisto v.0.46.2114 (Bray et al. 2016), and DESeq2 v.1.30.1 (Love et al. 2014) was used to normalise counts between samples within each species (Supplementary Tables 10–15). Unidentified variant genes with DESeq2 normalised expression values below ten for all developmental time points and biological replicates were flagged as putative pseudogenes (Supplementary Tables 4–6). Gene expression matrices in TPM and DESeq2 normalised values for all three species are available in the GitHub repository (see Data availability section).

Amplification of H2A.X variants in O. fusiformis

To validate the RNA-seq data and further confirm the expression of the h2ax1 (h2ax-y) and h2ax2 (h2ax-f) genes, we used four different combinations of two forward (Fw) and two reverse (Rv) primers (Supplementary Table 16) per gene to selectively amplify the transcript from a cDNA pool obtained from a mixture of developmental stages from O. fusiformis ranging from the zygote to the juvenile adult. For each gene, 1 µl cDNA was mixed with 1 µl of 10 µM Fw primer (1 or 2), 1 µl of 10 µM Rv primer (1 or 2), 0.5 µl of 10 mM dNTPs mix, 2.5 µl of 10× ThermoPol reaction buffer, and 0.15 µl of 5,000 U·ml^− 1 Taq DNA polymerase, and molecular biology-grade water for a final reaction volume of 25 µl. Samples were amplified by PCR in a thermocycler with a heated lid under the following conditions: initial denaturation at 95 ºC for 2 min; 30 cycles of 95 ºC for 15 s, 57 ºC for 15 s, and 68 ºC for 1 min; final extension at 68 ºC for 5 min and hold at 12 ºC. Length of resulting amplicons (Supplementary Table 17) was analysed via agarose electrophoresis.

Evolutionary analysis of histone H2A.X variants

An MSA with representative H2A.X sequences across Metazoa (Supplementary Table 18) was obtained using MAFFT v.7. (Katoh et al. 2019) as described above. To mine histone H2A.X variants across Eukarya, we performed a protein PHI-BLAST (Zhang et al. 1998) against the non-redundant protein sequences (nr) public database using the H2A.X-Y (search 1) and the H2A.X-F (search 2) orthologs from O. fusiformis, using SQ[DE]F and SQ[DE]Y as PHI patterns, for search 1, and search 2, respectively, and a maximum of 1,000 target sequences. This ensured the presence of the SQ[DE][FY] carboxyterminal (C-terminal) motif of H2A.X and simplified downstream analysis. Non-eukaryotic sequences were discarded (four for search 1, one for search 2) (Supplementary Tables 19, 20). Taxonomical information for sequences was obtained using the ETE 3 library (Huerta-Cepas et al. 2016: 3), and the eukaryotic supergroup was assigned manually based on the assigned phylum and according to current eukaryotic phylogeny (Burki et al. 2020) (Supplementary Tables 21–24). Sequences were aligned using MAFFT v.7. as described above. ML phylogenetic inference was performed in IQ-TREE v.2.0.3110 (Nguyen et al. 2015) with automatic model selection in ModelFinder Plus. The optimal model was JTTDCMut + R10 (Kosiol and Goldman 2005; Soubrier et al. 2012). Branch support values were extracted from 1,000 ultrafast BS (Hoang et al. 2018). GraPhlAn v.1.1.3 (Asnicar et al. 2015) was used to generate circular representations of the trees. To generate sequence logos for the C-terminal motif, H2A.X variants scored against pre-built Hidden Markov Models (HMM) were retrieved from the NCBI Histone Variants Database 2.0 (HistoneDB 2.0) (Draizen et al. 2016) for the Chordata, Streptophyta, Mollusca, Platyhelminthes, and Brachiopoda clades (these last three belonging to the Spiralia clade), aligned using MAFFT v.7. as described above. Alignments were trimmed to the last 12 positions of the alignment, corresponding to the positions 131 to 142 of the histone standardised nomenclature for H2A.X, using Jalview v.2.11.2.6 (Waterhouse et al. 2009). Trimmed alignments were turned into information content sequence logos using WebLogo v.2.8.2. (Crooks et al. 2004). Alignments and the corresponding alignment files in .fasta and .aln formats are available in the GitHub repository (see Data Availability section).

Histone modifier and reader genes mining and orthology assignment

To mine histone modifier orthologs, we split the search into five different analyses corresponding to the five families of interest: histone deacetylases (HDAC), histone acetyltransferases (HAT), type B HAT, histone demethylases (HDM), arginine-specific histone methyltransferases (PRMT), and lysine-specific histone methyltransferases (KMT) (Supplementary Table 25). Gene models corresponding to seven model organisms, namely Homo sapiens, Mus musculus, Danio rerio, Xenopus tropicalis, Caenorhabditis elegans, Drosophila melanogaster, and Saccharomyces cerevisiae were downloaded from NCBI and used to create BLAST local nucleotide databases (Camacho et al. 2009), alongside with the ones corresponding to the three annelid taxa of interest (Supplementary Table 26). Histone modifier protein sequences from these seven model organisms (Supplementary Table 27) were used as queries to find annelid orthologs following the mutual best BLAST hit approach (e value ≤ 10^− 5) (Supplementary Tables 28–32), obtaining 67, 97, 120, 40, and 155 unique annelid ortholog candidates for the HDAC, HAT, HDM, PRMT, and KMT families, respectively. HAT type A (nuclear) and type B (cytoplasmic) sequences were split into two groups. Appropriate outgroup sequences were chosen for each of the six alignments (Supplementary Table 33). MAFFT v.7 (Katoh et al. 2019) in the L-INS-I iterative refinement method and default scoring parameters were used to generate all six distinct multiple sequence alignments. For orthology assignment, six phylogenetic analyses were performed on selected candidate sequences, which included the longest isoform for each species-gene combination, given that it contained a properly aligned fragment within the family-specific common domain (see below) and that it was not located in a genomic locus of poor sequencing quality. Sequences were trimmed using Jalview v.2.11.2.6 (Waterhouse et al. 2009) to the family-specific domain(s), i.e. the histone deacetylase (DAC; HDAC), N-acetyltransferase (NAT; HAT type A and HAT type B), Jumonji C (JmjC; HDM), S-adenosylmethionine methyltransferase protein arginine methyltransferase-type (SAM MTase PRMT-type; PRMT), and the Su(var)3–9, Enhancer-of-zeste and Trithorax (SET; KMT) domains using the domain boundaries defined by ProSITE domain annotation for human HDAC11 (UniProt: Q96DB2; HDAC), KAT2A (UniProt: Q92830; HAT type A), NAA40 (UniProt: Q86UY6; HAT type B), KDM2A (UniProt: Q9Y2K7; HDM), PRMT1 (UniProt: Q99873; PRMT), and SUV39H1 (UniProt: O43463; KMT) proteins, respectively. Trimmed alignments were used for ML phylogenetic inference with automatic model selection using ModelFinder Plus in IQ-TREE v.2.0.3110 (Nguyen et al. 2015). The optimal models were LG + R7 (HDAC) (Le and Gascuel 2008; Soubrier et al. 2012), Q.insect + R5 (HAT type A) (Minh et al. 2021), LG + G4 (HAT type B) (Yang 1994), LG + R6 (HDM), LG + F + R5 (PRMT), and LG + R6 (KMT) (Kosiol and Goldman 2005). Bayesian inference in MrBayes v.3.2.7a (Ronquist and Huelsenbeck 2003) was also performed with an LG replacement matrix. Branch support values were extracted from 1,000 ultrafast BS (Hoang et al. 2018), and posterior probabilities were estimated from two runs with four chains of MCMC analyses run for 50,000,000 generations (80,000,000 for KMT) (Supplementary Table 34). All trees were composed and edited in FigTree v.1.4.4 (http://tree.bio.ed.ac.uk/software/figtree/). The resulting histone modifier orthology assignment is summarised in Supplementary Table 35 and described in detail for each species in Supplementary Table 36. Alignments and the corresponding alignment files in .fasta and .nexus formats are available in the GitHub repository (see Data Availability section).

PRMT6 sequence and architecture analysis

Candidate PRMT6 sequences were manually selected from a BLAST search of O. fusiformis’ PRMT6 protein sequence against the non-redundant (nr) protein database. All ten selected metazoan sequences, together with the annelid PRMT6 sequences (Supplementary Table 37), were aligned using MAFFT v.7 (Katoh et al. 2019) in the L-INS-I iterative refinement method and default scoring parameters. Alignment was trimmed to the conserved region of the protein between residues R43 and Y359, using Jalview v.2.11.2.6 (Waterhouse et al. 2009), using the nomenclature of human PRMT6 (UniProt: Q96LA8). Residues contained in the PRMT characteristic motif, the SAM-binding pocket, the arginine-binding pocket, as well as the SAM-interacting and arginine-interacting residues were derived from previous work (Schapira and Freitas 2014; Gupta et al. 2021). Domains and regions with functional relevance were obtained for 11 sequences using InterProScan (Blum et al. 2021). Transcriptomic validation of key genes was performed via direct visualisation of .bam files in Seqmonk v.1.48.1. Alignment and the corresponding alignment files in .fasta and .aln format are available in the GitHub repository (see Data Availability section).

Protein structure prediction

3D structural models of O. fusiformis H2A.X-F and H2A.X-Y proteins, as well as of O. fusiformis PRMT6, C. teleta PRMT6, and D. gyrociliatus PRMT6-A, PRMT6-B, PRMT6-C and PRMT6-D, were created using the AlphaFold Server (https://alphafoldserver.com/about) implementation of AlphaFold 3 (Abramson et al. 2024). Resulting predictions were rendered in PyMol v.3.0.3 (https://pymol.org/), where the relevant phosphate groups were added to Ser139 and Tyr142 of H2A.X-F and H2A.X-Y, where relevant. Protein structure predictions are available in the GitHub repository (see Data Availability section).

Gene clustering and comparative gene expression analyses

Full transcriptomes were clustered according to their normalised DESeq2 expression dynamics using soft k-means clustering (or soft clustering) in the mfuzz v.2.52 package (Futschik and Carlisle 2005) (Supplementary Tables 38–40). Transcripts that were not expressed at any developmental stage were discarded (154 out of 31,979 for O. fusiformis, 1,242 out of 41,238 for C. teleta, and 200 out of 17,387 for D. gyrociliatus). The Calinski-Harabasz index (Caliński and Harabasz 1974) was computed to determine an optimal number of 15 (O. fusiformis and C. teleta) and eight temporally co-regulated gene clusters (D. gyrociliatus) using the NbClust v.3.0.1 package (Charrad et al. 2014). For interspecies comparisons of histone modifiers expression dynamics (Supplementary Tables 41–43), clusters were classed as early/cleavage (O. fusiformis: 1–5; C. teleta: 1–5; D. gyrociliatus: 1–2) or late/post-cleavage (O. fusiformis: 7–15; C. teleta: 7–15; D. gyrociliatus: 4–6). Post-cleavage clusters in O. fusiformis and C. teleta were further subclassified as pre-larval (O. fusiformis: 7–9; C. teleta: 7–8) and post-larval (O. fusiformis: 11–15; C. teleta: 11–15). This way, four different quadrants were rendered for each species pairwise comparisons: early_{species A}–early_{species B}, early_{species A}–late_{species B}, late_{species A}–early_{species B}, and late_{species A}–late_{species B}. Clusters with peak expression in the female adult of D. gyrociliatus (7 and 8) were discarded for these purposes. To construct the gene sets of genes under heterochronic shifts, we only considered the histone modifiers with a single copy ortholog in both O. fusiformis and C. teleta (for comparisons between O. fusiformis and C. teleta, Supplementary Table 44) or in all three species (for comparisons between D. gyrociliatus and O. fusiformis or C. teleta, Supplementary Table 45).

Histones acid extraction

Histones were acid-extracted following an adaptation of a previous protocol (Dickman et al. 2013). 3× O. fusiformis adults and 15× C. teleta adults were food-deprived for 24 h in artificial seawater (ASW), washed in 1× PBS and spun down for 1.5 min at 3,000 × g. Animals were homogenised with a pellet pestle motor in 200 µl Tissue Extraction Buffer (TEB) (0.5% Triton X-100, 5 mM sodium butyrate, 2 mM phenylmethylsulphonyl fluoride (PMSF), 0.02% sodium azide, supplemented with 1× cOmplete™ EDTA-free Protease Inhibitor Cocktail Tablets (Roche)). Homogenates were layered over 2.5 ml of a 1.8 M sucrose solution in TEB and centrifuged at 49,000× g and 4°C for 1 h in a rate-zonal centrifugation. Supernatants were removed, and the pellets were resuspended in 0.5 ml TEB, transferred to a clean tube, and pelleted at 21,000 × g and 4°C for 2 min. Histones were extracted from the nuclear pellets in 1.2 ml 0.5 M HCl for 48 h at 4°C. Crude extracts were centrifuged at 6,500 × g and 4°C for 10 min to remove insoluble debris. 4× cycles of ultrafiltration in Amicon Ultra-0.5 Centrifugal Filter Units with Ultracel-10 regenerated cellulose membranes (10 kDa nominal molecular weight limit, NWML) (Millipore) using mQ water as exchange buffer were performed following the manufacturer’s recommendations, to remove acid and concentrate histones to 20–25 µl. Aliquots were obtained at various steps to assess protein size and purity via sodium dodecyl sulphate-polyacrylamide gel electrophoresis (SDS-PAGE). For O. fusiformis, an additional cleaning step was introduced before the ultrafiltration to remove high molecular weight co-extracted proteins using Amicon Ultra-0.5 Centrifugal Filter Units with Ultracel-30 regenerated cellulose membranes (30 kDa, NWML) (Millipore).

Histone derivatisation and preparation

Histones were prepared for mass spectrometry as previously described (Cole et al. 2021). Briefly, histones underwent chemical derivatisation by adding 10 µl of 100 mM ammonium bicarbonate pH 8.0 and 4 µl of ammonium hydroxide to 10 µg of histone sample then adding 10 µl of propionic anhydride in isopropanol and ammonium hydroxide were then used to maintain a pH higher than 8.0. After a 15-min incubation at 37°C, samples were dried down in a vacuum centrifuge (Concentrator plus, Eppendorf), and the whole process was repeated. Samples were then re-suspended in 40 µl of 100 mM ammonium bicarbonate and digested with trypsin overnight. Digestion was stopped by adding glacial acetic acid and freezing at − 80 ºC for 5 min. Samples were vacuum-centrifuged and then underwent two more rounds of proprionylation. Lastly, we used a Hypersep™ Hypercarb™ tip to desalt the samples following the manufacturer’s protocol (Thermo Fisher Scientific) (Minshull et al. 2016).

Liquid-chromatography tandem-mass spectrometry (LC-MS/MS)

Histone samples were in 0.1% trifluoroacetic acid (TFA) and analysed on an Ultimate 3000 online nano-LC system with a PepMap300 C18 trapping column (Thermo Fisher Scientific) coupled to a Q Exactive HF Orbitrap (Thermo Fisher Scientific). Peptides were eluted onto a 50 cm × 75 µm Easy‐spray PepMap C18 analytical column at 35°C at a flow rate of 300 nl·min^− 1 using a gradient of 3–25% over 55 min, and then 25–60% until 81 min. Loading solvent was 0.1% TFA and 3% acetonitrile (ACN), solvent A comprised 0.1% formic acid (FA) and 3% ACN, and solvent B was 0.1% FA and 80% ACN. Samples were run in data-independent acquisition mode. Histone posttranslational modifications (hPTM) were identified, and their relative abundance was quantified in EpiProfile 2.0, with manual verification of key hPTMs in Skyline, as previously described (Cole et al. 2019) (Supplementary Tables 46, 47).

Competing interests

Francisco M. Martín-Zamora is an employee of Altos Labs.

Author Contribution

JMM-D, PJH, and FMM-Z conceived and designed the study. JC and MJD performed LC-MS/MS analyses. KG, AMC-B, and RDD contributed to animal culture and in vitro fertilisations. FMM-Z performed all other lab work and analyses described here. JMM-D and FMM-Z drafted the manuscript. All authors critically read and commented on the manuscript.

Acknowledgement

We thank Ben Hadley and Natasha Sandikyan for their help with the initial characterisation of histone modifiers, as well as the technical support staff at the Department of Biology at Queen Mary University of London. This study employed computational resources from the High-performance computing facility Apocrita at Queen Mary University of London. This work was funded by grants to José M. Martín-Durán from the European Union Horizon 2020 Framework Programme under the European Research Council (ERC Starting Grant agreement number 801669), to Francisco M. Martín-Zamora (ASSEMBLE Plus Grant agreement number 730984), and to Paul J. Hurd from the Biotechnology and Biological Sciences Research Council (BBSRC Grant agreement numbers BB/L023164/1 and BB/V009311/1).

Data Availability

Supplementary Figures and Supplementary Tables are available in the Supplementary Information section of this manuscript and in a public GitHub repository (see below). Accession codes and unique identifiers to previously publicly available datasets we used for this study are listed in Supplementary Table 1 (sequence identifiers used in histone genes mining), Supplementary Table 2 (sequence identifiers used in histone genes orthology assignment), Supplementary Table 8 (ATAC-seq datasets used in chromatin accessibility profiling), Supplementary Table 9 (RNA-seq datasets used in gene expression profiling), Supplementary Table 18 (sequence identifiers used in evolutionary analysis of histone H2A.X variants), Supplementary Table 26 (genome files used in histone modifier genes mining and orthology assignment), Supplementary Table 27 (sequence identifiers used in histone modifier genes mining and orthology assignment) and Supplementary Table 37 (sequence identifiers used in PRMT6 sequence and architecture analysis). The histone variants database, that is, the NCBI Histone Variants Database 2.0 (HistoneDB 2.0), can be found at https://www.ncbi.nlm.nih.gov/research/histonedb/. The AlphaFold Server can be accessed at https://alphafoldserver.com/. Updated filtered gene models, transcripts and protein files, genome annotation reports, alignment files, gene expression files, protein structure predictions, and any other relevant additional files are all publicly available in the GitHub repository associated with this study: https://github.com/ChemaMD/HistonesAnnelida. All custom code used is available in the GitHub repository associated with this study: https://github.com/ChemaMD/HistonesAnnelida.

Code availability

All custom code used is available in the GitHub repository associated with this study: https://github.com/ChemaMD/HistonesAnnelida.

Abouelmagd A, Ageely HM. 2013. Basic Genetics. Universal.
Abramson J, Adler J, Dunger J, Evans R, Green T, Pritzel A, Ronneberger O, Willmore L, Ballard AJ, Bambrick J, et al. Accurate structure prediction of biomolecular interactions with AlphaFold 3. Nature. 2024;630:493–500.
Albig W, Warthorst U, Drabent B, Prats E, Cornudella L, Doenecke D. Mytilus edulis Core Histone Genes Are Organized in Two Clusters Devoid of Linker Histone Genes. J Mol Evol. 2003;56:597–606.
Amatori S, Tavolaro S, Gambardella S, Fanelli M. The dark side of histones: genomic organization and role of oncohistones in cancer. Clin Epigenetics. 2021;13:71.
Andäng M, Hjerling-Leffler J, Moliner A, Lundgren TK, Castelo-Branco G, Nanou E, Pozas E, Bryja V, Halliez S, Nishimaru H, et al. Histone H2AX-dependent GABAA receptor regulation of stem cell proliferation. Nature. 2008;451:460–4.
Anderson L, Gomes MR, daSilva LF, Pereira A, da SA, Mourão MM, Romier C, Pierce R, Verjovski-Almeida S. Histone deacetylase inhibition modulates histone acetylation at gene promoter regions and affects genome-wide gene transcription in Schistosoma mansoni. PLoS Negl Trop Dis. 2017;11:e0005539.
Armakola M, Ruvkun G. 2019. Regulation of Caenorhabditis elegans neuronal polarity by heterochronic genes. Proceedings of the National Academy of Sciences 116:12327–12336.
Asnicar F, Weingart G, Tickle TL, Huttenhower C, Segata N. Compact graphical representation of phylogenetic data and metadata with GraPhlAn. PeerJ. 2015;3:e1029.
Audia JE, Campbell RM. Histone Modifications and Cancer. Cold Spring Harb Perspect Biol. 2016;8:a019521.
Ausio J, Van Holde KE. The histones of the sperm of Spisula solidissima include a novel, cysteine-containing H-1 histone. Cell Differ. 1988;23:175–89.
Baldi S, Becker PB. The variant histone H2A.V of Drosophila—three roles, two guises. Chromosoma. 2013;122:245–58.
Bannister AJ, Kouzarides T. Regulation of chromatin by histone modifications. Cell Res. 2011;21:381–95.
Blanc RS, Richard S. Arginine Methylation: The Coming of Age. Mol Cell. 2017;65:8–24.
Blanco E, González-Ramírez M, Alcaine-Colet A, Aranda S, Croce LD. The Bivalent Genome: Characterization, Structure, and Regulation. Trends Genet. 2020;36:118–31.
Blum M, Chang H-Y, Chuguransky S, Grego T, Kandasaamy S, Mitchell A, Nuka G, Paysan-Lafosse T, Qureshi M, Raj S, et al. The InterPro protein families and domains database: 20 years on. Nucleic Acids Res. 2021;49:D344–54.
Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30:2114–20.
Bönisch C, Schneider K, Pünzeler S, Wiedemann SM, Bielmeier C, Bocola M, Eberl HC, Kuegel W, Neumann J, Kremmer E, et al. H2A.Z.2.2 is an alternatively spliced histone H2A.Z variant that causes severe nucleosome destabilization. Nucleic Acids Res. 2012;40:5951–64.
Bray NL, Pimentel H, Melsted P, Pachter L. Near-optimal probabilistic RNA-seq quantification. Nat Biotechnol. 2016;34:525–7.
Bryson TD, De Ioannes P, Valencia-Sánchez MI, Henikoff JG, Talbert PB, Lee R, La Scola B, Armache K-J, Henikoff S. A giant virus genome is densely packaged by stable nucleosomes within virions. Mol Cell. 2022;82:4458–e44705.
Buendía-Monreal M, Gillmor CS. The Times They Are A-Changin’: Heterochrony in Plant Development and Evolution. Front Plant Sci. 2018;9:1349.
Burki F, Roger AJ, Brown MW, Simpson AGB. The New Tree of Eukaryotes. Trends Ecol Evol. 2020;35:43–55.
Burns R, Pechenik J. Transcriptomic Basis of Metamorphic Competence in the Salt-Marsh-Dwelling Polychaete Capitella teleta. Biol Bull. 2017;232:158–70.
Buschbeck M, Hake SB. Variants of core histones and their roles in cell fate decisions, development and cancer. Nat Rev Mol Cell Biol. 2017;18:299–314.
Buschbeck M, Uribesalgo I, Wibowo I, Rué P, Martin D, Gutierrez A, Morey L, Guigó R, López-Schier H, Di Croce L. The histone variant macroH2A is an epigenetic regulator of key developmental genes. Nat Struct Mol Biol. 2009;16:1074–9.
Bush KM, Yuen BT, Barrilleaux BL, Riggs JW, O’Geen H, Cotterman RF, Knoepfler PS. Endogenous mammalian histone H3.3 exhibits chromatin-related functions during development. Epigenetics Chromatin. 2013;6:7.
Butler JS, Koutelou E, Schibler AC, Dent SYR. Histone-modifying enzymes: Regulators of developmental decisions and drivers of human disease. Epigenomics. 2012;4:163–77.
Cabezas-Cruz A, Lancelot J, Caby S, Oliveira G, Pierce RJ. Epigenetic control of gene function in schistosomes: a source of therapeutic targets? Front Genet. 2014;5:317.
Caliński T, Harabasz J. A dendrite method for cluster analysis. Commun Stat. 1974;3:1–27.
Calpena E, Palau F, Espinós C, Galindo MI. Evolutionary History of the Smyd Gene Family in Metazoans: A Framework to Identify the Orthologs of Human Smyd Genes in Drosophila and Other Animal Species. PLoS ONE. 2015;10:e0134106.
Camacho C, Coulouris G, Avagyan V, Ma N, Papadopoulos J, Bealer K, Madden TL. BLAST+: architecture and applications. BMC Bioinformatics. 2009;10:421.
Carrillo-Baltodano AM, Seudre O, Guynes K, Martín-Durán JM. 2021. Early embryogenesis and organogenesis in the annelid Owenia fusiformis. EvoDevo 12:5.
Chang B, Chen Y, Zhao Y, Bruick RK. JMJD6 Is a Histone Arginine Demethylase. Science. 2007;318:444–7.
Charrad M, Ghazzali N, Boiteau V, Niknafs A. NbClust: An R Package for Determining the Relevant Number of Clusters in a Data Set. J Stat Softw. 2014;61:1–36.
Chi P, Allis CD, Wang GG. Covalent histone modifications — miswritten, misinterpreted and mis-erased in human cancers. Nat Rev Cancer. 2010;10:457–69.
Cole J, Angyal A, Emes RD, Mitchell TJ, Dickman MJ, Dockrell DH. Pneumolysin Is Responsible for Differential Gene Expression and Modifications in the Epigenetic Landscape of Primary Monocyte Derived Macrophages. Front Immunol. 2021;12:573266.
Cole J, Hanson EJ, James DC, Dockrell DH, Dickman MJ. Comparison of data-acquisition methods for the identification and quantification of histone post-translational modifications on a Q Exactive HF hybrid quadrupole Orbitrap mass spectrometer. Rapid Commun Mass Spectrom. 2019;33:897–906.
Colmenares SU, Swenson JM, Langley SA, Kennedy C, Costes SV, Karpen GH. Drosophila Histone Demethylase KDM4A Has Enzymatic and Non-enzymatic Roles in Controlling Heterochromatin Integrity. Dev Cell. 2017;42:156–e1695.
Cook PJ, Ju BG, Telese F, Wang X, Glass CK, Rosenfeld MG. Tyrosine dephosphorylation of H2AX modulates apoptosis and survival decisions. Nature. 2009;458:591–6.
Cosseau C, Wolkenhauer O, Padalino G, Geyer KK, Hoffmann KF, Grunau C. (Epi)genetic Inheritance in Schistosoma mansoni: A Systems Approach. Trends Parasitol. 2017;33:285–94.
Creyghton MP, Cheng AW, Welstead GG, Kooistra T, Carey BW, Steine EJ, Hanna J, Lodato MA, Frampton GM, Sharp PA et al. 2010. Histone H3K27ac separates active from poised enhancers and predicts developmental state. Proceedings of the National Academy of Sciences 107:21931–21936.
Crooks GE, Hon G, Chandonia J-M, Brenner SE. WebLogo: A Sequence Logo Generator. Genome Res. 2004;14:1188–90.
Dainat J, Hereñú D, Pucholt P. 2021. AGAT: Another Gff Analysis Toolkit to handle annotations in any GTF/GFF format. Zenodo [Internet]. https://zenodo.org/record/4637977
Dattani A, Kao D, Mihaylova Y, Abnave P, Hughes S, Lai A, Sahu S, Aboobaker AA. Epigenetic analyses of planarian stem cells demonstrate conservation of bivalent histone modifications in animal stem cells. Genome Res. 2018;28:1543–54.
Davis FC, Shelton JC, Ingham LD. Nucleotide sequence of the Urechis caupo core histone gene tandem repeat. DNA Seq. 1992;2:247–56.
Degnan BM, Degnan SM. How larvae and life cycles evolve. Trends Genet. 2023;39:528–30.
Dickman MJ, Kucharski R, Maleszka R, Hurd PJ. Extensive histone post-translational modification in honey bees. Insect Biochem Mol Biol. 2013;43:125–37.
Dijkwel Y, Tremethick DJ. The Role of the Histone Variant H2A.Z in Metazoan Development. J Dev Biology. 2022;10:28.
Dimmic MW, Rest JS, Mindell DP, Goldstein RA. rtREV: An Amino Acid Substitution Matrix for Inference of Retrovirus and Reverse Transcriptase Phylogeny. J Mol Evol. 2002;55:65–73.
Dobreva MP, Camacho J, Abzhanov A. Time to synchronize our clocks: Connecting developmental mechanisms and evolutionary consequences of heterochrony. J Experimental Zool Part B: Mol Dev Evol. 2022;338:87–106.
Draizen EJ, Shaytan AK, Mariño-Ramírez L, Talbert PB, Landsman D, Panchenko AR. 2016. HistoneDB 2.0: a histone database with variants—an integrated resource to explore histones and their variants. Database 2016:baw014.
Duncan EM, Chitsazan AD, Seidel CW, Sánchez Alvarado A. Set1 and MLL1/2 Target Distinct Sets of Functionally Different Genomic Loci In Vivo. Cell Rep. 2015;13:2741–55.
Eirín-López JM, González-Tizón AM, Martinez A, Méndez J. Molecular and Evolutionary Analysis of Mussel Histone Genes (Mytilus spp.): Possible Evidence of an Orphon Origin for H1 Histone Genes. J Mol Evol. 2002;55:272–83.
Ferrand J, Rondinelli B, Polo SE. Histone Variants: Guardians of Genome Integrity. Cells. 2020;9:2424.
Formery L, Lowe CJ. 2023. Integrating Complex Life Cycles in Comparative Developmental Biology. Annu Rev Genet 57.
Futschik ME, Carlisle B. Noise-robust soft clustering of gene expression time-course data. J Bioinform Comput Biol. 2005;3:965–88.
Gaspar-Maia A, Qadeer ZA, Hasson D, Ratnakumar K, Adrian Leu N, Leroy G, Liu S, Costanzi C, Valle-Garcia D, Schaniel C, et al. MacroH2A histone variants act as a barrier upon reprogramming towards pluripotency. Nat Commun. 2013;4:1565.
del Gaudio R, Potenza N, Stefanoni P, Chiusano ML, Geraci G. Organization and Nucleotide Sequence of the Cluster of Five Histone Genes in the Polichaete Worm Chaetopterus variopedatus: First Record of a H1 Histone Gene in the Phylum Annelida. J Mol Evol. 1998;46:64–73.
Ghazy E, Abdelsalam M, Robaa D, Pierce RJ, Sippl W. Histone Deacetylase (HDAC) Inhibitors for the Treatment of Schistosomiasis. Pharmaceuticals. 2022;15:80.
Giaimo BD, Ferrante F, Herchenröther A, Hake SB, Borggrefe T. The histone variant H2A.Z in gene regulation. Epigenetics Chromatin. 2019;12:1–22.
Gilbert N, Allan J. Supercoiling in DNA and chromatin. Curr Opin Genet Dev. 2014;25:15–21.
González-Romero R, Rivera-Casas C, Frehlick LJ, Méndez J, Ausió J, Eirín-López JM. Histone H2A (H2A.X and H2A.Z) Variants in Molluscs: Molecular Characterization and Potential Implications For Chromatin Dynamics. PLoS ONE. 2012;7:e30006.
Grau-Bové X, Navarrete C, Chiva C, Pribasnig T, Antó M, Torruella G, Galindo LJ, Lang BF, Moreira D, López-Garcia P, et al. A phylogenetic and proteomic reconstruction of eukaryotic chromatin evolution. Nat Ecol Evol. 2022;6:1007–23.
Greenfield A, Carrel L, Pennisi D, Philippe C, Quaderi N, Siggers P, Steiner K, Tam PPL, Monaco AP, Willard HF, et al. The UTX Gene Escapes X Inactivation in Mice and Humans. Hum Mol Genet. 1998;7:737–42.
Günesdogan U, Jäckle H, Herzig A. A genetic system to assess in vivo the functions of histones and histone modifications in higher eukaryotes. EMBO Rep. 2010;11:772–6.
Gupta S, Kadumuri RV, Singh AK, Chavali S, Dhayalan A. 2021. Structure, Activity and Function of the Protein Arginine Methyltransferase 6. Life 11:951.
Guynes K, Sarre LA, Carrillo-Baltodano AM, Davies BE, Xu L, Liang Y, Martín-Zamora FM, Hurd PJ, de Mendoza A, Martín-Durán JM. Annelid methylomes reveal ancestral developmental and aging-associated epigenetic erosion across Bilateria. Genome Biol. 2024;25:204.
van Helden PD, Strickland WN, von Brandt WF HOLT C. The Complete Amino-Acid Sequence of Histone H2B from the Mollusc Patella granatina. Eur J Biochem. 1979;93:71–8.
Henikoff S, Smith MM. Histone variants and epigenetics. Cold Spring Harb Perspect Biol. 2015;7:a019364–019364.
Hoang DT, Chernomor O, von Haeseler A, Minh BQ, Vinh LS. UFBoot2: Improving the Ultrafast Bootstrap Approximation. Mol Biol Evol. 2018;35:518–22.
Hocher A, Laursen SP, Radford P, Tyson J, Lambert C, Stevens KM, Montoya A, Shliaha PV, Picardeau M, Sockett RE, et al. Histones with an unconventional DNA-binding mode in vitro are major chromatin constituents in the bacterium Bdellovibrio bacteriovorus. Nat Microbiol. 2023;8:2006–19.
Huerta-Cepas J, Serra F, Bork P. ETE 3: Reconstruction, Analysis, and Visualization of Phylogenomic Data. Mol Biol Evol. 2016;33:1635–8.
Irwin NAT, Richards TA. Self-assembling viral histones are evolutionary intermediates between archaeal and eukaryotic nucleosomes. Nat Microbiol. 2024;9:1713–24.
Jack APM, Bussemer S, Hahn M, Pünzeler S, Snyder M, Wells M, Csankovszki G, Solovei I, Schotta G, Hake SB. H3K56me3 Is a Novel, Conserved Heterochromatic Mark That Largely but Not Completely Overlaps with H3K9me3 in Both Regulation and Localization. PLoS ONE. 2013;8:e51765.
Jambhekar A, Dhall A, Shi Y. Roles and regulation of histone methylation in animal development. Nat Rev Mol Cell Biol. 2019;20:625–41.
Jang C-W, Shibata Y, Starmer J, Yee D, Magnuson T. Histone H3.3 maintains genome integrity during mammalian development. Genes Dev. 2015;29:1377–92.
Johnson DS, Mortazavi A, Myers RM, Wold B. Genome-Wide Mapping of in Vivo Protein-DNA Interactions. Science. 2007;316:1497–502.
Kasinathan S, Orsi GA, Zentner GE, Ahmad K, Henikoff S. High-resolution mapping of transcription factor binding sites on native chromatin. Nat Methods. 2014;11:203–9.
Katoh K, Rozewicki J, Yamada KD. MAFFT online service: multiple sequence alignment, interactive sequence choice and visualization. Brief Bioinform. 2019;20:1160–6.
Kaya-Okur HS, Janssens DH, Henikoff JG, Ahmad K, Henikoff S. Efficient low-cost chromatin profiling with CUT&Tag. Nat Protoc. 2020;15:3264–83.
Kaya-Okur HS, Wu SJ, Codomo CA, Pledger ES, Bryson TD, Henikoff JG, Ahmad K, Henikoff S. CUT&Tag for efficient epigenomic profiling of small samples and single cells. Nat Commun. 2019;10:1930.
Kerbl A, Fofanova EG, Mayorova TD, Voronezhskaya EE, Worsaae K. Comparison of neuromuscular development in two dinophilid species (Annelida) suggests progenetic origin of Dinophilus gyrociliatus. Front Zool. 2016a;13:49.
Kerbl A, Martín-Durán JM, Worsaae K, Hejnol A. Molecular regionalization in the compact brain of the meiofaunal annelid Dinophilus gyrociliatus (Dinophilidae). EvoDevo. 2016b;7:20.
Keyte AL, Smith KK. Heterochrony and developmental timing mechanisms: Changing ontogenies in evolution. Semin Cell Dev Biol. 2014;34:99–107.
Khorasanizadeh S. The Nucleosome: From Genomic Organization to Genomic Regulation. Cell. 2004;116:259–72.
Kimura H. Histone modifications for human epigenome analysis. J Hum Genet. 2013;58:439–45.
Kixmoeller K, Allu PK, Black BE. The centromere comes into focus: from CENP-A nucleosomes to kinetochore connections with the spindle. Open Biology. 2020;10:200051.
Kmiecik D, Couppez M, Belaiche D, Sautiere P. Primary structure of histone H2A from nucleated erythrocyte of the marine worm Sipunculus nudus Presence of two forms of H2A in the sipunculid chromatin. Eur J Biochem. 1983;135:113–21.
Kosiol C, Goldman N. Different Versions of the Dayhoff Rate Matrix. Mol Biol Evol. 2005;22:193–9.
Krawetz SA, Sellos D, Dixon CH. Analysis of the core histone gene cluster of the annelid Platynereis dumerilii. DNA Seq. 1993;4:29–35.
Lawrence M, Huber W, Pagès H, Aboyoun P, Carlson M, Gentleman R, Morgan MT, Carey VJ. Software for Computing and Annotating Genomic Ranges. PLoS Comput Biol. 2013;9:e1003118.
Le SQ, Gascuel O. An improved general amino acid replacement matrix. Mol Biol Evol. 2008;25:1307–20.
Lee JS, Smith E, Shilatifard A. The Language of Histone Crosstalk. Cell. 2010;142:682–5.
Lei B, Berger F. H2A Variants in Arabidopsis: Versatile Regulators of Genome Activity. Plant Commun. 2020;1:100015.
Lemmens BBLG, Tijsterman M. DNA double-strand break repair in Caenorhabditis elegans. Chromosoma. 2011;120:1–21.
Li D, Taylor DH, van Wolfswinkel JC. PIWI-mediated control of tissue-specific transposons is essential for somatic cell differentiation. Cell Rep. 2021;37:109776.
Liu W. Epigenetics in Schistosomes: What We Know and What We Need Know. Front Cell Infect Microbiol. 2016;6:149.
Liu W, Ma Q, Wong K, Li W, Ohgi K, Zhang J, Aggarwal AK, Rosenfeld MG. Brd4 and JMJD6-Associated Anti-Pause Enhancers in Regulation of Transcriptional Pause Release. Cell. 2013;155:1581–95.
Lopez-Delisle L, Rabbani L, Wolff J, Bhardwaj V, Backofen R, Grüning B, Ramírez F, Manke T. pyGenomeTracks: reproducible plots for multivariate genomic datasets. Bioinformatics. 2020;37:422–3.
Lorković ZJ, Park C, Goiser M, Jiang D, Kurzbauer M-T, Schlögelhofer P, Berger F. Compartmentalization of DNA Damage Response between Heterochromatin and Euchromatin Is Mediated by Distinct H2A Histone Variants. Curr Biol. 2017;27:1192–9.
Lorton BM, Shechter D. Cellular consequences of arginine methylation. Cell Mol Life Sci. 2019;76:2933–56.
Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15:550.
Male G, von Appen A, Glatt S, Taylor NMI, Cristovao M, Groetsch H, Beck M, Müller CW. Architecture of TFIIIC and its role in RNA polymerase III pre-initiation complex assembly. Nat Commun. 2015;6:7387.
Martin BJE, Brind’Amour J, Kuzmin A, Jensen KN, Liu ZC, Lorincz M, Howe LJ. Transcription shapes genome-wide histone acetylation patterns. Nat Commun. 2021;12:210.
Martín-Durán JM, Marlétaz F. Unravelling spiral cleavage. Development. 2020;147:dev181081.
Martín-Durán JM, Vellutini BC, Marlétaz F, Cetrangolo V, Cvetesic N, Thiel D, Henriet S, Grau-Bové X, Carrillo-Baltodano AM, Gu W, et al. Conservative route to genome compaction in a miniature annelid. Nat Ecol Evol. 2021;5:231–42.
Martín-Zamora FM, Davies BE, Donnellan RD, Guynes K, Martín-Durán JM. Functional genomics in Spiralia. Brief Funct Genomics. 2023a;22:487–97.
Martín-Zamora FM, Liang Y, Guynes K, Carrillo-Baltodano AM, Davies BE, Donnellan RD, Tan Y, Moggioli G, Seudre O, Tran M, et al. Annelid functional genomics reveal the origins of bilaterian life cycles. Nature. 2023b;615:105–10.
Martire S, Banaszynski LA. The roles of histone variants in fine-tuning chromatin organization and function. Nat Rev Mol Cell Biol. 2020;21:522–41.
McKay DJ, Klusza S, Penke TJR, Meers MP, Curry KP, McDaniel SL, Malek PY, Cooper SW, Tatomer DC, Lieb JD, et al. Interrogating the Function of Metazoan Histones using Engineered Gene Clusters. Dev Cell. 2015;32:373–86.
McNamara KJ. Heterochrony: the Evolution of Development. Evo Edu Outreach. 2012;5:203–18.
Meers MP, Bryson TD, Henikoff JG, Henikoff S. 2019. Improved CUT&RUN chromatin profiling tools.Parker S, Weigel D, editors. eLife 8:e46314.
Mihaylova Y, Abnave P, Kao D, Hughes S, Lai A, Jaber-Hijazi F, Kosaka N, Aboobaker AA. Conservation of epigenetic regulation by the MLL3/4 tumour suppressor in planarian pluripotent stem cells. Nat Commun. 2018;9:3633.
Millán-Zambrano G, Burton A, Bannister AJ, Schneider R. Histone post-translational modifications — cause and consequence of genome function. Nat Rev Genet. 2022;23:563–80.
Minh BQ, Dang CC, Vinh LS, Lanfear R. QMaker: Fast and Accurate Method to Estimate Empirical Models of Protein Evolution. Syst Biol. 2021;70:1046–60.
Minshull TC, Cole J, Dockrell DH, Read RC, Dickman MJ. Analysis of histone post translational modifications in primary monocyte derived macrophages using reverse phase×reverse phase chromatography in conjunction with porous graphitic carbon stationary phase. J Chromatogr A. 2016;1453:43–53.
Moosmann A, Campsteijn C, Jansen PW, Nasrallah C, Raasholm M, Stunnenberg HG, Thompson EM. Histone variant innovation in a rapidly evolving chordate lineage. BMC Evol Biol. 2011;11:208.
Morgan MAJ, Shilatifard A. Reevaluating the roles of histone-modifying enzymes and their associated chromatin modifications in transcriptional regulation. Nat Genet. 2020;52:1271–81.
Moss EG. Heterochronic Genes and the Nature of Developmental Time. Curr Biol. 2007;17:R425–34.
Nagel S, Grossbach U. Histone H1 Genes and Histone Gene Clusters in the Genus Drosophila. J Mol Evol. 2000;51:286–98.
Neiro J, Sridhar D, Dattani A, Aboobaker A. Identification of putative enhancer-like elements predicts regulatory networks active in planarian adult stem cells.Bronner ME. editor eLife. 2022;11:e79675.
Nelson DL, Cox MM, Hoskins AA. Lehninger Principles of Biochemistry. Macmillan Learning; 2021.
Nguyen LT, Schmidt HA, Von Haeseler A, Minh BQ. IQ-TREE: A fast and effective stochastic algorithm for estimating maximum-likelihood phylogenies. Mol Biol Evol. 2015;32:268–74.
Ouchi Y, Yamamoto J, Iwamoto T. The Heterochronic Genes lin-28a and lin-28b Play an Essential and Evolutionarily Conserved Role in Early Zebrafish Development. PLoS ONE. 2014;9:e88086.
Pascual-Carreras E, Marín-Barba M, Castillo-Lara S, Coronel-Córdoba P, Magri MS, Wheeler GN, Gómez-Skarmeta JL, Abril JF, Saló E, Adell T. Wnt/β-catenin signalling is required for pole-specific chromatin remodeling during planarian regeneration. Nat Commun. 2023;14:1–16.
Pereira SL, Grayling RA, Lurz R, Reeve JN. Archaeal nucleosomes. Proc Natl Acad Sci USA. 1997;94:12633–7.
Pertea G, Pertea M. 2020. GFF Utilities: GffRead and GffCompare. F1000Research 9:304.
Picard MAL, Boissier J, Roquis D, Grunau C, Allienne J-F, Duval D, Toulza E, Arancibia N, Caffrey CR, Long T, et al. Sex-Biased Transcriptome of Schistosoma mansoni: Host-Parasite Interaction, Genetic Determinants and Epigenetic Regulators Are Associated with Sexual Differentiation. PLoS Negl Trop Dis. 2016;10:e0004930.
Picard MAL, Vicoso B, Roquis D, Bulla I, Augusto RC, Arancibia N, Grunau C, Boissier J, Cosseau C. Dosage Compensation throughout the Schistosoma mansoni Lifecycle: Specific Chromatin Landscape of the Z Chromosome. Genome Biol Evol. 2019;11:1909–22.
Quénet D, Dalal Y. The CENP-A nucleosome: a dynamic structure and role at the centromere. Chromosome Res. 2012;20:465–79.
Ramírez F, Ryan DP, Grüning B, Bhardwaj V, Kilpert F, Richter AS, Heyne S, Dündar F, Manke T. deepTools2: a next generation web server for deep-sequencing data analysis. Nucleic Acids Res. 2016;44:W160–5.
Rink JC. Stem cell systems and regeneration in planaria. Dev Genes Evol. 2013;223:67–84.
Robinson JT, Thorvaldsdóttir H, Winckler W, Guttman M, Lander ES, Getz G, Mesirov JP. Integrative genomics viewer. Nat Biotechnol. 2011;29:24–6.
Rogakou EP, Boon C, Redon C, Bonner WM. Megabase Chromatin Domains Involved in DNA Double-Strand Breaks in Vivo. J Cell Biol. 1999;146:905–16.
Ronquist F, Huelsenbeck JP. MrBayes 3: Bayesian phylogenetic inference under mixed models. Bioinformatics. 2003;19:1572–4.
Roquis D, Rognon A, Chaparro C, Boissier J, Arancibia N, Cosseau C, Parrinello H, Grunau C. Frequency and mitotic heritability of epimutations in Schistosoma mansoni. Mol Ecol. 2016;25:1741–58.
Roquis D, Taudt A, Geyer KK, Padalino G, Hoffmann KF, Holroyd N, Berriman M, Aliaga B, Chaparro C, Grunau C, et al. Histone methylation changes are required for life cycle progression in the human parasite Schistosoma mansoni. PLoS Pathog. 2018;14:e1007066.
Sakai A, Schwartz BE, Goldstein S, Ahmad K. Transcriptional and Developmental Functions of the H3.3 Histone Variant in Drosophila. Curr Biol. 2009;19:1816–20.
Sankar A, Mohammad F, Sundaramurthy AK, Wang H, Lerdrup M, Tatar T, Helin K. Histone editing elucidates the functional roles of H3K27 methylation and acetylation in mammals. Nat Genet. 2022;54:754–60.
Schapira M, Freitas RFde. Structural biology and chemistry of protein arginine methyltransferases. Med Chem Commun. 2014;5:1779–88.
Schmidl C, Rendeiro AF, Sheffield NC, Bock C. ChIPmentation: fast, robust, low-input ChIP-seq for histones and transcription factors. Nat Methods. 2015;12:963–5.
Seaver EC, Thamm K, Hill SD. Growth patterns during segmentation in the two polychaete annelids, Capitella sp. I and Hydroides elegans: comparisons at distinct life history stages. Evol Dev. 2005;7:312–26.
Seifert-Davila W, Girbig M, Hauptmann L, Hoffmann T, Eustermann S, Müller CW. Structural insights into human TFIIIC promoter recognition. Sci Adv. 2023;9:eadh2019.
Sellos D. The histones isolated from the sperm of the oyster Crassostrea gigas. Cell Differ. 1985;17:183–92.
Sellos D, Krawetz SA, Dixon GH. Organization and complete nucleotide sequence of the core-histone-gene cluster of the annelid Platynereis dumerilii. Eur J Biochem. 1990;190:21–9.
Shechter D, Chitta RK, Xiao A, Shabanowitz J, Hunt DF, Allis CD. A distinct H2A.X isoform is enriched in Xenopus laevis eggs and early embryos and is phosphorylated in the absence of a checkpoint. Proc Natl Acad Sci USA. 2009;106:749–54.
Shpargel KB, Sengoku T, Yokoyama S, Magnuson T. UTX and UTY Demonstrate Histone Demethylase-Independent Function in Mouse Embryonic Development. PLoS Genet. 2012;8:e1002964.
Shroff R, Arbel-Eden A, Pilch D, Ira G, Bonner WM, Petrini JH, Haber JE, Lichten M. Distribution and Dynamics of Chromatin Modification Induced by a Defined DNA Double-Strand Break. Curr Biol. 2004;14:1703–11.
Shumate A, Salzberg SL. Liftoff: accurate mapping of gene annotations. Bioinformatics. 2021;37:1639–43.
Simakov O, Marletaz F, Cho SJ, Edsinger-Gonzales E, Havlak P, Hellsten U, Kuo DH, Larsson T, Lv J, Arendt D, et al. Insights into bilaterian evolution from three spiralian genomes. Nature. 2013;493:526–31.
Skene PJ, Henikoff S. An efficient targeted nuclease strategy for high-resolution mapping of DNA binding sites.Reinberg D. editor eLife. 2017;6:e21856.
Skvortsova K, Iovino N, Bogdanović O. Functions and mechanisms of epigenetic inheritance in animals. Nat Rev Mol Cell Biol. 2018;19:774–90.
Slack F, Ruvkun G. Temporal Pattern Formation by Heterochronic Genes. Annu Rev Genet. 1997;31:611–34.
Soubrier J, Steel M, Lee MSY, Der Sarkissian C, Guindon S, Ho SYW, Cooper A. The Influence of Rate Heterogeneity among Sites on the Time Dependence of Molecular Rates. Mol Biol Evol. 2012;29:3345–58.
Tessarz P, Kouzarides T. Histone core modifications regulating nucleosome structure and dynamics. Nat Rev Mol Cell Biol. 2014;15:703–8.
Turinetto V, Giachino C. Multiple facets of histone variant H2AX: a DNA double-strand-break marker with several biological functions. Nucleic Acids Res. 2015;43:2489–98.
Verma P, Waterbury CKM, Duncan EM. Set1 Targets Genes with Essential Identity and Tumor-Suppressing Functions in Planarian Stem Cells. Genes. 2021;12:1182.
Völker-Albert M, Bronkhorst A, Holdenrieder S, Imhof A. Histone Modifications in Stem Cell Development and Their Clinical Implications. Stem Cell Rep. 2020;15:1196–205.
Wang W-L, Anderson LC, Nicklay JJ, Chen H, Gamble MJ, Shabanowitz J, Hunt DF, Shechter D. Phosphorylation and arginine methylation mark histone H2A prior to deposition during Xenopus laevis development. Epigenetics Chromatin. 2014;7:22.
Wang X, Ausió J. Histones are the major chromosomal protein components of the sperm of the nemerteans Cerebratulus californiensis and Cerebratulus lacteus. J Exp Zool. 2001;290:431–6.
Waterhouse AM, Procter JB, Martin DMA, Clamp M, Barton GJ. Jalview Version 2–a multiple sequence alignment editor and analysis workbench. Bioinformatics. 2009;25:1189–91.
Whetstine JR, Nottke A, Lan F, Huarte M, Smolikov S, Chen Z, Spooner E, Li E, Zhang G, Colaiacovo M, et al. Reversal of Histone Lysine Trimethylation by the JMJD2 Family of Histone Demethylases. Cell. 2006;125:467–81.
Workman JL, Kingston RE. Alteration of Nucleosome Structure as a Mechanism of Transcriptional Regulation. Annu Rev Biochem. 1998;67:545–79.
Wouters-Tyrou D, Martin-Ponthieu A, Briand G, Sautiere P, Biserte G. The Amino-Acid Sequence of Histone H2A from Cuttlefish Sepia officinalis. Eur J Biochem. 1982;124:489–98.
Wu T, Liu Y, Wen D, Tseng Z, Tahmasian M, Zhong M, Rafii S, Stadtfeld M, Hochedlinger K, Xiao A. Histone Variant H2A.X Deposition Pattern Serves as a Functional Epigenetic Mark for Distinguishing the Developmental Potentials of iPSCs. Cell Stem Cell. 2014;15:281–94.
Xiao A, Li H, Shechter D, Ahn SH, Fabrizio LA, Erdjument-Bromage H, Ishibe-Murakami S, Wang B, Tempst P, Hofmann K, et al. WSTF regulates the H2A.X DNA damage response via a novel tyrosine kinase activity. Nature. 2009;457:57–62.
Yang Z. Maximum likelihood phylogenetic estimation from DNA sequences with variable rates over sites: Approximate methods. J Mol Evol. 1994;39:306–14.
Yu Y, Song C, Zhang Q, DiMaggio PA, Garcia BA, York A, Carey MF, Grunstein M. Histone H3 Lysine 56 Methylation Regulates DNA Replication through Its Interaction with PCNA. Mol Cell. 2012;46:7–17.
Yun M, Wu J, Workman JL, Li B. Readers of histone modifications. Cell Res. 2011;21:564–78.
Zhang C, Lin H, Zhang Y, Xing Q, Zhang J, Zhang D, Liu Y, Chen Q, Zhou T, Wang J, et al. BRPF1 bridges H3K4me3 and H3K23ac in human embryonic stem cells and is essential to pluripotency. iScience. 2023;26:105939.
Zhang Z, Miller W, Schäffer AA, Madden TL, Lipman DJ, Koonin EV, Altschul SF. Protein sequence similarity searches using patterns as seeds. Nucleic Acids Res. 1998;26:3986–90.
Zhao Z, Shilatifard A. Epigenetic modifications of histones in cancer. Genome Biol. 2019;20:245.
Zhou VW, Goren A, Bernstein BE. Charting histone modifications and the functional organization of mammalian genomes. Nat Rev Genet. 2011;12:7–18.
Zou Y, Chiu H, Zinovyeva A, Ambros V, Chuang C-F, Chang C. Developmental Decline in Neuronal Regeneration by the Progressive Change of Two Intrinsic Timers. Science. 2013;340:372–6.

No competing interests reported.

Download PDF

Editor assigned by journal
24 Oct, 2024
Submission checks completed at journal
18 Oct, 2024
First submitted to journal
17 Oct, 2024

You are reading this latest preprint version

A dynamic histone-based chromatin regulatory toolkit underpins genome and developmental evolution in an invertebrate clade

Status:

Version 1

Abstract

Figures

Introduction

Results and Discussion

Annelids exhibit a conserved histone repertoire

Genomic organisation and gene structure of histones

Histone genes are in dynamic hyperaccessible chromatin regions

Histone gene expression dynamics in Annelida

Histone H2A.X variants have evolved in parallel in Metazoa

Histone modifiers show a complex evolution throughout annelid diversification

Histone deacetylases

Histone demethylases

Histone acetyltransferases

Lysine-specific histone methyltransferases

Arginine-specific histone methyltransferases

PRMT6 expansions led to domain fusions and likely catalytically dead enzymes

Life cycle correlates with heterochronies of histone modifiers

Adult annelids harbour distinct histone modification enrichments

Conclusions

Materials and Methods

Adult culture

Histone genes mining

Histone genes orthology assignment

Histone gene structure analysis

Chromatin accessibility profiling

Gene expression profiling

Amplification of H2A.X variants in O. fusiformis

Evolutionary analysis of histone H2A.X variants

Histone modifier and reader genes mining and orthology assignment

PRMT6 sequence and architecture analysis

Protein structure prediction

Gene clustering and comparative gene expression analyses

Histones acid extraction

Histone derivatisation and preparation

Liquid-chromatography tandem-mass spectrometry (LC-MS/MS)

Declarations

Competing interests

Author Contribution

Acknowledgement

Data Availability

Code availability

References

Additional Declarations

Supplementary Files

Status:

Version 1