{'fraction_non_alphanumeric': 0.03554955563055462, 'fraction_numerical': 0.019827877151535604, 'mean_word_length': 5.288177557697409, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 1, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A phylogenetic and proteomic reconstruction of eukaryotic chromatin evolution

Histones and associated chromatin proteins have essential functions in eukaryotic genome organization and regulation. Despite this fundamental role in eukaryotic cell biology, we lack a phylogenetically-comprehensive understanding of chromatin evolution. Here, we combine comparative proteomics and genomics analysis of chromatin in eukaryotes and archaea. Proteomics uncovers the existence of histone post-translational modifications in Archaea. However, archaeal histone modifications are scarce, in contrast with the highly conserved and abundant marks we identify across eukaryotes. Phylogenetic analysis reveals that chromatin-associated catalytic functions (e.g., methyltransferases) have pre-eukaryotic origins, whereas histone mark readers and chaperones are eukaryotic innovations. We show that further chromatin evolution is characterized by expansion of readers, including capture by transposable elements and viruses. Overall, our study infers detailed evolutionary history of eukaryotic chromatin: from its archaeal roots, through the emergence of nucleosome-based regulation in the eukaryotic ancestor, to the diversification of chromatin regulators and their hijacking by genomic parasites.


Introduction
The access to genetic information in eukaryotes is controlled by a manifold nucleoproteic interface called chromatin. This nucleosomal chromatin environment defines a repressive ground state for transcription and other DNA-templated processes in eukaryotic genomes 1,2 . Multiple components associated with chromatin underlie elaborate eukaryotic genome regulation, allowing the differential access to genetic information in time/space and the maintenance of the resulting regulatory states [3][4][5][6] . Moreover, chromatin-based regulation is essential in repressing parasitic genomic elements, like transposons and viruses [7][8][9][10][11] .
The main protein components of eukaryotic chromatin are histones. All eukaryotes have four major types of histones (H2A, H2B, H3 and H4), which are combined as an octamer to form the basic repetitive unit of the chromatin: the nucleosome. Canonical histones are among the most highly conserved proteins across eukaryotes 12 and, in addition, unique histone variants (paralogs of one of the four major histone types) are found in many species, often associated with particular regulatory states [13][14][15][16][17] . Histone chemical modifications, including acetylations and methylations play a central role in genome regulation and transgenerational epigenetic inheritance 3,18-21 . These chemical moieties, known as histone post-translational modifications (hPTMs), are added and removed by specific enzymes ('writers', e.g., histone methyltransferases or acetylases; and 'erasers', e.g., histone demethylases and deacetylases). Some hPTMs (e.g., most acetylations) have a generic effect on nucleosome stability, while others are bound by specific proteins or protein complexes. These are often referred to as 'readers' and include proteins like HP1, which binds to H3K9me3, as well as a myriad of other proteins encoding Chromo, PHD, Tudor and Bromo structural domains, among others [22][23][24] . Finally, nucleosome remodellers (like SNF2 proteins) and histone chaperones are additional important players in chromatin regulation, by mediating chromatin opening, nucleosomal assembly, and histone variant interchanges [25][26][27][28] .
All eukaryotes studied to date possess histone-based chromatin organization, with the sole exception of dinoflagellates, which nonetheless encode for histone proteins in their genomes 29 . Beyond eukaryotes, histones have also been identified in Archaea, where they have been shown to form nucleosomal structures [30][31][32][33] . However, unlike eukaryotic histones, the few archaeal histones experimentally characterized so far (i) generally lack disordered N-terminal tails; (ii) do not have any known post-translational modifications 34 ; and (iii) do not seem to impose a widespread, genome-wide repressive transcriptional ground state 33,35 . Thus, chromatin-based elaborate genome regulation is often considered a eukaryotic innovation 36,37 .
From a phylogenetic perspective, our understanding of chromatin components and processes derives from a very small set of organisms, essentially animal, fungal and plant model species plus a few parasitic unicellular eukaryotes. Additional efforts have sampled specific aspects of chromatin regulation, such as histone modifications or their genome-wide distribution, in non-model animal species 38,39 , fungi (Neurospora crassa and Fusarium graminearum) 40,41 , and five other eukaryotes: the unicellular holozoan Capsaspora owczarzaki 42 , the dinoflagellate Hematodiunium sp. 29 , the brown alga Ectocarpus siliculosus 43 , the amoebozoan Dictyostelium discoideum 44 , and the ciliate Tetrahymena thermophila 45,46 . However, these organisms represent a tiny fraction of eukaryotic diversity.
Hence, we lack a systematic understanding of the evolution of eukaryotic chromatin modifications and components 47 .
In order to infer the origin and evolutionary diversification of eukaryotic chromatin, we performed a joint comparative analysis of histone proteomics data from 30 different eukaryotic and archaeal taxa, including new data for 23 species. In parallel, we analyzed the complement of chromatin-associated gene families in an additional 172 eukaryotic genomes and transcriptomes. This comprehensive taxon sampling includes representatives of all major eukaryotic lineages, as well as multiple free-living members of enigmatic early-branching eukaryotes (e.g., jakobids, malawimonads, Meteora sp. and ancyromonads, as well as Collodictyonida, Rigifilida and Mantamonadida (CRuMS); Fig. 1a). In addition, in order to trace the pre-eukaryotic origins of these chromatin gene families, we systematically searched for orthologs in archaeal, bacterial and viral genomes. Specifically, we reconstructed the evolutionary history of enzymes involved in chromatin modification and remodelling; as well as the conservation of the hPTMs effected by these enzymes. Our comparative genomics and proteomics suggest a concurrent and early origin of canonical histones, a core of quasi-universal hPTMs, and their corresponding enzymatic effectors. We also identify independent expansions in hPTM reader gene families across eukaryotes and document evidence of the capture of these reader domains by parasitic genomic elements. Overall, this work provides a phylogenetically-informed framework to classify and compare chromatin components across the eukaryotic tree of life, and to further investigate the evolution of hPTM-mediated genome regulation.
histone-like DNA binding domains, which are widely distributed in eukaryotes (POLE3, POLE4, DR1) and/or archaea (NFYB). Further analysis of the genomic distribution of these histone genes shows a frequent occurrence of H3-H4 and H2A-H2B pairs in head-to-head orientation (5' to 5'), strongly indicating co-regulation across eukaryotes (Extended Data Fig.  1b,c and Supplementary Data 2).
Next, we investigated the distribution and conservation of hPTMs across major eukaryotic groups and Archaea, including methylations, acetylations, crotonylations, phosphorylations, and ubiquitylations. To this end, histones from 19 different eukaryotic species were extracted, chemically derivatized 49 and analyzed by mass-spectrometry ( Fig. 1c and Supplementary Data 3), adding to previously available hPTM proteomics data for additional seven species. Our extensive taxon sampling covers all major eukaryotic groups, as well as hitherto unsampled early-diverging eukaryotic lineages-such as the malawimonad Gefionella okellyi, the discoban Naegleria gruberi, or the ancyromonad Fabomonas tropica -, thus providing a comprehensive comparative framework for evolutionary inference.
We focused first on hPTMs present in canonical histones, as defined by their highly conserved N-terminal regions, phylogenetic analyses, and sequence similarity to curated reference canonical histones ( Fig. 1d; see Methods). hPTMs are detected in all canonical histones from all species. After correcting by sequence coverage, we observe that hPTMs are particularly abundant in H3 canonical histones (median = 23.5 hPTMs per species, mean = 24.3), compared with H2A, H2B and H4 (medians between 6.5 and 9, means between 9.5 and 13.4; Extended Data Fig. 2a). Holozoan canonical H2As (Homo sapiens, Sycon ciliatum and Capsaspora owczarzaki) represent an exception to this trend and contain similar number of modifications to H3s in these species. We also examined the reproducibility of hPTM detection across replicate samples, showing that the majority of hPTMs (87.5%) can be found in more than one sample (Extended Data Fig. 2b,c). Despite this, it is worth emphasizing that our data may contain false negatives, beyond the lack of coverage for particular residues that we systematically report. For example, some marks might be globally too scarce in the nucleosomes of a particular species, while other modifications like phosphorylations and ubiquitination are difficult to detect by mass-spectrometry without dedicated peptide-enrichment protocols.
Canonical H3 and H4 N-terminal tails contain the majority of phylogenetically-conserved hPTMs, in stark contrast with the relative paucity of conserved hPTMs in canonical H2A and H2B. A striking example of paneukaryotic conservation comes from the acetylation of the H4 K5, K8, K12 and K16 residues (Fig. 1d, second panel), all of which mark gene expression-permissive chromatin environments in multiple eukaryotic species 22 . A similar conservation pattern is observed in the acetylation of a group of N-terminal H3 lysines (K9, K14, K18, K23, K27) associated with similar functions, while other H3 acetylations are only found in a few species (e.g., residues K4, K56 and K79). While acetylations are highly conserved, only seven histone H3/H4 methylations are broadly conserved across eukaryotic lineages: H3K4me1/2/3, H3K9me1/2/3, H3K27me1/2/3, H3K36me1/2/3, H3K37me1/2/3 and, more sparsely, H3K79me1/2 and H4K20me1. Many of these broadly conserved marks have conserved roles in demarcating active (e.g., H3K4me) and repressive chromatin states (e.g., H3K9me and H3K27me) 22,42,50 . The scarcity of conserved hPTMs in H2A and H2B canonical histones can partially explained by their higher degree of sequence divergence (Fig. 1e), which is reflected in many non-homologous lysine residues (Fig. 1d). But even among homologous positions, we found little evidence of conservation, with the exception of H2A K5ac (associated to active promoters 51 ) and, in fewer species, methylation of H2A K5 and H2B K5. Finally, we were also able to identify phosphorylations in serine and threonine residues and a few instances of ubiquitylation. In general, these marks show more restricted phylogenetic distributions than lysine acetylation or methylation, even in the tightly conserved H3 and H4 histones. We can identify conserved phosphorylations in H2A T120 and S122, which are shared by most opisthokonts, and the ubiquitylation of H2A K119 only in some holozoan species.
Mass-spectrometry analysis detected histone variants in all species included in our study, suggesting that they are relatively abundant in the chromatin of these eukaryotes (Fig. 1e). Most of these variants are lineage-specific, with the exception of the paneukaryotic variants H2A.Z, H3/cenH3 and H3.3; and the macroH2A variant found in holozoans and Meteora sp. (belonging to an orphan eukaryotic lineage). Interestingly, we find hPTMs in the vast majority of detected variants, both conserved and lineage-specific, particularly acetylations and methylations ( Fig. 1e and Extended Data Fig. 2d). Overall, our comparative proteomic analysis suggests the existence of a highly conserved set of canonical hPTMs of ancestral eukaryotic origin in H3 and H4, which co-exists with less conserved hPTMs in H2A, H2B, and lineage-specific modifications in variant histones.

Archaeal histone post-translational modifications
In contrast with the paneukaryotic distribution of histones, sequence searches show that only a fraction of archaeal genomes encode for histones (28.1% of the taxa here examined; Fig. 2a). Archaeal histones exhibit a patchy phylogenetic distribution, similar to other gene families shared with eukaryotes 52 . Among others, histones are present in Euryarchaeota, the TACK superphylum and Asgard archaea 12,[53][54][55][56] . Asgard are generally are considered to be the closest known archaeal relatives of eukaryotes 57,58 , although this sister-group relationship has been challenged by some studies 59 . Our extended sampling revealed that Asgard archaea histones, particularly in the Lokiarchaeota and Heimdallarchaeota clades 55 , often have lysine-rich N-terminal tails in the manner of eukaryotic histones (Fig. 2a-c).
These Asgard histones appear to be conserved across multiple taxa, albeit without direct sequence similarity compared to canonical eukaryotic histones (Extended Data Fig. 1d). When compared against eukaryotic sequences classified in HistoneDB 48 , these archaeal histones clearly cluster in a separate group and are most similar to either eukaryotic H4 or, to a lesser degree, H3 canonical histones, in line with previous findings 12,55,60 .
To identify potential archaeal hPTMs, we performed proteomics analysis of histones in three Euryarchaeota (the Methanobacteriota Methanobrevibacter cuticularis and the Halobacteriota Methanospirillum stamsii and Methanosarcina spelaei) and one Thaumarchaeota species (Nitrososphaera viennensis; Fig. 2b). Mass-spectrometry detects histone proteins in all of them: 2-4 in the euryarchaeotes (with 27-90% protein coverage) and one in the thaumarchaeote (80% protein coverage), including homologs with N-terminal tails encoded by each of the three euryarchaeotes in our survey (22-40 aa, 0.09-28 lysines per residue; Fig. 2c). Moreover, this proteomics analysis finds evidence of hPTMs in archaeal histones. However, in comparison with eukaryotic histones, hPTMs are extremely scarce in archaeal histones. Specifically, we identify no hPTMs in N. viennensis and M. spelaei (one and two histones detected, respectively), three acetylations and one methylation in M. stamsii (in three out of four histones detected), and one acetylation and two methylations in M. cuticularis (in two out of four histones; Fig. 2b, top). Interestingly, we find conserved lysine residues with shared modifications in M. stamsii and M. cuticularis (methylation in K54 and acetylation in K57; Fig. 2b, bottom). This result indicates that highly-abundant hPTMs represent a eukaryotic innovation, likely linked to dynamic nucleosomal regulation in eukaryotes but not in Archaea.

Taxonomic distribution of chromatin-associated proteins
hPTMs are deposited and removed by specific modifying enzymes ('writers' and 'erasers'), while 'reader' protein domains found in diverse proteins bind and recognize specific hPTMs. For example, Bromo and Chromo domains bind acetylated and methylated lysine residues, respectively. In addition, the control of histone loading/eviction from specific genomic loci is mediated by chromatin remodellers, like SNF2 proteins 27 , and histone chaperones 26 . To date, the classification and evolutionary analysis of this chromatin machinery has been based on biased, partial taxonomic samplings and has not employed phylogenetic methods 61 (with rare exceptions 12,27 ), often resulting in inaccurate orthologous relationships and confounded classification and naming schemes.
We sought to obtain a systematic, phylogenetics-based classification of histone remodellers, chaperones, readers, and modifiers in order to understand the evolutionary history of eukaryotic chromatin (Fig. 3a). To this end, we (i) compiled a taxa-rich dataset of 172 eukaryotic genomes and transcriptomes, covering all major eukaryotic supergroups and devoting particular attention to early-branching, non-parasitic lineages (Supplementary Data 1), as well as genomic data from 4,226 Archaea, 24,886 Bacteria and 185,579 viral taxa; (ii) defined a protein structural domain as a proxy for each gene family (Supplementary Data 4) and retrieved all genes in these genomes that contained these domains; and (iii) inferred accurate orthology groups from phylogenetic analyses of each gene class (next section).
We examined the taxonomic distribution and abundance of the major gene classes (Fig.  3b,c). Many domains with chromatin-associated functions in eukaryotes are also present in Archaea and Bacteria, albeit with scattered phylogenetic distributions ( Fig. 3b and Extended Data Fig. 3a,b). Families with prokaryotic homologs include mostly catalytic gene classes (writer, eraser and remodeller enzymes), whereas readers and histone chaperones are virtually absent from prokaryotes (Fig. 3b). Histone fold-encoding genes constitute a case in point for this patchy distribution of chromatin proteins in prokaryotes: they are present in most archaeal phyla, but are absent in about half of the sampled genomes within each (Fig.  3b). Yet, there is a qualitative difference between the phylogenetic distribution of archaeal and bacterial chromatin-associated gene classes: whereas archaeal histones tend to co-occur with chromatin-associated gene classes, the bacterial complement of writers and erasers is much less conserved and is uncorrelated with the extremely rare presence of histone-like genes (Fig. 3d).
Within eukaryotes, most gene structural classes associated with chromatin functions are ubiquitously distributed across all lineages here surveyed, supporting an early eukaryotic origin for the core chromatin machinery ( Fig. 3b and Extended Data Fig. 3d). In fact, the total number of chromatin writer, eraser and remodeller enzymes remains remarkably stable across eukaryotes (Fig. 3e). The only exception is the marked increase in genes encoding reader domains observed in lineages exhibiting complex multicellularity: animals, streptophyte plants, and, to a lesser degree, phaeophyte brown algae (Stramenopila). This occurs partially due to the addition of new gene classes (e.g., SAWADEE in the Plantae s.l. + Cryptista lineage, or ADD_DNMT3 in bilaterians and cnidarians), but also via the expansion of ancient, widely-distributed reader gene classes (e.g., Tudor, PHD, Chromo or Bromo domains). These taxonomic patterns indicate that chromatin modifying and remodelling catalytic activities originated in prokaryotes, while reader and chaperone structural domains are eukaryotic innovations.

Phylogenetics of chromatin modifiers and remodellers
To gain detailed insights into the origin and evolution of chromatin gene families, we used phylogenetic analysis to define orthology groups from paneukaryotic gene trees. We surveyed 172 eukaryotic species and defined a total of 1,713 gene families (orthogroups), 95% of which were conserved in two or more high-ranking taxonomic groups (as listed in Fig. 1a), and which included 51,426 genes in total (Supplementary Data 5). We annotated each gene family according to known members from eukaryotic model species. For simplicity, we use a human-based naming scheme throughout the present manuscript (unless otherwise stated), but we also provide a dictionary of orthologs in three additional model species (Arabidopsis thaliana, Saccharomyces cerevisiae and Drosophila melanogaster; see Supplementary Data 5). This phylogenetic classification scheme of eukaryotic chromatin gene families, as well as the sequences and associated phylogenetic trees, can be explored and retrieved in an interactive database: https://sebe-lab.shinyapps.io/chromatin_evolution We first investigated the potential pre-eukaryotic origins of these gene families/orthogroups by comparing their phylogenetic distance to prokaryotic sequences and to other eukaryotic orthogroups (Fig. 4a). Most eukaryotic gene families are more closely related to other eukaryotes than to prokaryotic sequences, supporting the idea that writers, erasers, remodellers and readers diversified within the eukaryotic lineage, as previously noted for histones 12 . This analysis also reveals a substantial fraction of eukaryotic gene families with close orthogroups in Archaea and Bacteria, which pinpoints components that were (i) inherited from a prokaryotic ancestor during eukaryogenesis; (ii) laterally transferred between eukaryotes and prokaryotes at later stages; or (iii) a combination of both phenomena. For example, we identified a well-supported sister-group relationship between the eukaryotic SIRT7 deacetylase and a clade of Asgard archaea Sirtuin enzymes (Heimdallarchaeota and Lokiarchaeota), a topology compatible with an archaeal origin or ancient transfers to/from Asgard and eukaryotes 62 ; whereas SIRT6 appears nested within other eukaryotic sequences (Fig. 4b, left). Likewise, the KAT14 acetylase is more closely related to bacterial enzymes than to other eukaryotic acetylases (Fig. 4b, right).
Next, we mapped the phylogenetic distribution of orthogroups in order to infer the origin and diversification of individual chromatin gene families ( Fig. 4c and Extended Data Fig.  4a). Using probabilistic inference of ancestral gene content, we reconstruct a rich Last Eukaryotic Common Ancestor (LECA) complement of chromatin-associated gene families: 65 acetylases (amongst which 61 were conserved in at least two of the most deeply sampled eukaryotic early-branching lineages, namely Amorphea, Diaphoretickes, and Discoba); 20 deacetylases (19 in these early-branching eukaryotic lineages); 59 methyltransferases (55); 42 demethylases (38); 33 remodellers (33); and 25 chaperones (18) (Fig. 4c and Supplementary Data 5). The subsequent evolution of these families is characterized by relative stasis, with few new orthologous families emerging in later-branching eukaryotic lineages. Notable exceptions include the origin of KAT5 deacetylases and KMT5B/C SET methyltransferases in Opisthokonta; KAT8 and SIRT7 in Holozoa; and Viridiplantae-specific deacetylases (homologs of A. thaliana HDA7 and HDA14 deacetylases) and SETs (A. thaliana PTAC14); among others.
In spite of their broad distributions across eukaryotes, many chromatin modifier families exhibit variation in their protein domain architectures, likely conferring them functional properties such as distinct binding preferences (Extended Data Fig. 4b). For example, most CREBBP/EP300 acetylases consist of a catalytic HAT_KAT11 domain and two TAZ and ZZ zinc finger domains, but different lineages have acquired different reader domains: an acetylation-reading Bromo domain in holozoans and stramenopiles, PHD in plants and some stramenopiles, and no known reader domains in other lineages (e.g., in the fungal orthologs of the S. cerevisiae protein RTT109). A similar pattern is apparent in SET methyltransferase families sharing a core catalytic domain (SET) harboring variable DNAand chromatin-interacting domains -animal SETDB1/2 homologs have MBD domains that bind CpG methylated DNA, while plants have SAD_SAR domains with the same function; and holozoan ASH1L homologs encode Bromo and BAH readers, whereas phaeophytes encode PHD domains (Extended Data Fig. 4b). Other architectures, however, are much more conserved, as exemplified by the presence of Tudor-knot and MYST zinc finger domains in most KAT5 deacetylases; or the ubiquitous co-occurrence of Helicase-C and SNF2_N domains in most remodellers (Extended Data Fig. 4b).
Specific examples of evolutionarily conserved chromatin gene families include the catalytic core and the subunits of well-studied chromatin complexes 63 like PRC1 (RING1/AB, PCGF), PRC2 (EZH1/2, SUZ12, EED, RBBP4/7) and Trithorax/MLL (MLL1/2/3/4, WRD5, ASH2L, RBBP5, DPY-30; Fig. 4d,e). However, when we compared the distribution of these complexes with the hPTMs they are related to, we found a generally poor co-occurrence ( Fig. 4f-h). For example, organisms like Dictyostelium discoideum and Creolimax fragrantissima lack EZH1/2 orthologs, but we detected H3K27me3 in these species; while Thecamonas trahens and Naegleria gruberi lack Dot1 orthologs but have H3K79me marks. A poor correlation is also observed between the occurrence of H3K9me and that of SUV39H1 orthologs. An exception to this pattern is the ubiquitous distribution of H4K16ac and the acetylase family KAT5/8 64 (Fig. 4h). These patterns suggest that the specificity between hPTMs and their writers might not be completely conserved across eukaryotes, with distinct members of the same gene classes (e.g., methyltransferases) performing similar roles. In this context, reading domains present in writing/erasing enzymes (directly in the same protein or as part of multi-protein complexes) are likely to play a major role in the re-purposing of chromatin catalytic activities.

Evolutionary expansion of chromatin readers
Multiple protein structural domains have been involved in the recognition of hPTMs, such as Bromo domains binding to acetylated lysines or Chromo, MBT and Tudor domains binding to methylated lysines 23,24 . These are generally small domains and can be found both as stand-alone proteins as well as in combination with other domains, often catalytic activities such as hPTM writers, erasers and remodellers. Thus, they are central in the establishment of functional connections between chromatin states. To understand the contribution of these reading domains to the evolutionary diversification of chromatin networks, we studied in detail the phylogeny and protein architecture of reader domains across eukaryotes.
We quantified the co-occurrence frequency of reader and catalytic domains, finding (i) that most reader domains are present in genes without writer, eraser or remodeller domains (87%, Fig. 5a); and (ii) that most cases of reader-catalytic co-occurrence involve PHD, Chromo and Bromo domains (Extended Data Fig. 5a). For example, the conserved architecture of the paneukaryotic CHD3/4/5 re-modellers includes Chromo readers in most species and PHD domains specifically in animals and plants (Extended Data Fig.  4b). Likewise, PHD domains are often present in the KMT2A/B and KMT2C/D SET methyltransfrase; and the ASH1L family has recruited Bromo and BAH domains in holozoans, and PHD in multicellular stramenopiles (Extended Data Fig. 4b). In spite of these redundancies, reader families typically have independent evolutionary histories, as illustrated by the fact that most reader domain-containing genes encode only one such domain (92%, Extended Data Fig. 5b).
We next performed phylogenetic analyses of individual reader domains and reconstructed the gains and losses of these reader gene families/orthogroups (Fig. 5a). Compared to the relative stasis of catalytic enzyme families, this reader-centric analysis revealed a strikingly different evolutionary pattern of lineage-specific bursts of innovation, particularly amongst PHD, Chromo and Bromo genes, as well as Tudor in animals ( Fig. 5a and Extended Data Fig. 5c). PHD, Chromo and Bromo families also appeared as the most abundant in the reconstructed LECA reader domain repertoire, which amounted to 89 gene families ( Fig.  5a, left). The distribution of gene family ages in extant species also corroborates that more readers have emerged at evolutionarily more recent nodes of the tree of life than catalytic gene families (Fig. 5b).

Co-option of chromatin machinery by transposable elements
Further examination of the domain co-occurrence networks of readers revealed that Chromo and PHD domains are often present together with protein domains found in transposable elements (TEs; Fig. 5c and Supplementary Data 6), including retrotransposons (e.g., retrotranscriptases and integrases; orange modules in Fig. 5c) and DNA transposons (e.g., DNA binding domains and transposases; red modules). It is known that some TEs show insertion-preferences associated to specific chromatin states 65 , often mediated by direct chromatin tethering mechanisms 66 . For example, the Chromo domain of the MAGGY gypsy retrotransposon of the fungus Magnaporthe grisea targets H3K9me regions 67 . Reciprocally, some protein domains of TE origin, often DNA-binding domains, have been co-opted into chromatin and transcriptional regulators 68 . Thus, we decided to explore in detail the occurrence of chromatin-associated domain (readers, but also catalytic domains) linked to TEs in the 172 eukaryotic genomes in our dataset (Fig. 5d). Moreover, we used available RNA-seq datasets in many of these species to validate some of these TE fusions (Fig. 5d-e). A fully validated fusion gene would (i) come from a non-discontinuous gene model in the original assembly, and (ii) have evidence of expression, with reads mapping along the entire region between the TE-associated domain and the chromatin-associated domain (Extended Data Fig. S6).
We identified 823 predicted gene models containing both chromatin-and TE-associated domains (Fig. 5d). Whilst these TE fusions were not exclusive of reader domains, most such fusions involved PHD and Chromo-encoding genes; followed by SNF2_N remodellers, SET methyltransferases, and others. An homology search against a database of eukaryotic TEs revealed that most of these candidate TE fusions could be aligned to known retrotransposons or DNA transposons. For example, by way of validation, our analysis identifies the SETMAR human gene, a previously-described fusion between a SET methyltranferase and a Mariner-class DNA transposon 69 . Overall, 31% of the candidate fusion genes were supported by valid gene models according to our stringent criteria (Fig. 5d). Interestingly, we find very few cases of hypothetical fusions between TEs and Bromo domains, which recognize K acetylations and are otherwise highly abundant across eukaryotes, and none of them is validated by RNA-seq data. This could be explained by the detrimental effect of targeting TE insertions to sites of active chromatin demarcated by histone acetylations, such as promoter and enhancer elements.
Some of these validated fusions have a broad phylogenetic distribution (Fig. 5e), such as a Gypsy-ERV retrotransposon with a C-terminal Chromo domain (Unk. Chromo 2.1 in Fig.   5e) that is widely distributed in animals and various microbial eukaryotes, and contains dozens of paralogs in vertebrate Danio rerio or the charophyte Chara braunii, many of which are expressed. Another widespread Gypsy-ERV retrotransposon with a Chromo domain is present in multiple expressed and highly similar copies in the fungus Rhizopus delemar (Fig. 5f,e), suggesting a successful colonization of this genome by this TE. By contrast, other TE fusions are taxonomically restricted to one or few related species, such as the fusion of hAT activator DNA transposons with Chromo CBX and CDY readers in the sponge Ephydatia muelleri; or multiple instances of fusions with Chromo and PHD readers in cnidarians. A common fusion in cnidarians involves different retrotransposon classes with PHD domains orthologous to the PYGO1/2 protein (Fig. 5e), which is known to recognize specifically H3K4me 70 . Globally, this analysis reveals that recruitment of chromatin reading and even modifying domains by TE has occurred in many eukaryotic species, in a way that might facilitate the evasion from suppressing mechanisms in the host genomes as suggested by the expansion of Chromo-fused TEs in the genomes of Chara braunii (Viridiplantae), Chromera velia (Alveolata) and Rhizopus delemar (Fungi).

Chromatin components in viral genomes
In addition to TEs, chromatin is also involved in the suppression of another type of genomic parasites: viruses. Some chromatin-related genes, including histones, have been found in viral genomes, especially among the nucleocytoplasmic large DNA viruses -also known as giant viruses. Eukaryotic core histones have been even hypothesized to have evolved from giant virus homologs, after the discovery that certain Marseilleviridae genomes encoded deeply-diverging orthologs of the four canonical histones 71 . These viral histones have been recently shown to form nucleosome-like particles that package viral DNA 72,73 .
We analyzed the distribution and abundance of chromatin-related protein domains among viruses, including data from 1,816 giant virus genomes. Based on structural domain searches, we identified 2,163 viral chromatin-related proteins ( Fig. 5g and Supplementary Data 6). The majority of these proteins are encoded by giant viruses (55%), followed by Caudovirales (37%). Among these two groups, only giant virus genomes encode histones -specifically, the Iridoviridae, Marseilleviridae, Mimiviridae, Pithoviridae, and Phycodnaviridae families. Concordantly with previous studies 74 , we also identify remodellers in all giant virus families; as well as less abundant components of the chromatin writer/eraser/reader toolkit (Fig. 5g).
We then investigated the phylogenetic affinities of these viral chromatin proteins, starting with histones (Fig. 5h). Our analysis recovers the phylogenetic affinity of Marseilleviridae histones with specific eukaryotic histone families 71 , and makes this pattern extensive to Mimiviridae, Iridoviridae, and Pithoviridae giant viruses (Fig. 5h), with the caveat of the ambiguous clustering of the H4-like viral histones with either H4 eukaryotic or archaeal HMfB genes. In all these lineages, we identify genes encoding two histone-fold domains orthologous to H2B + H2A (inset table in Fig. 5h), whereas the H4 + H3 histone doublet genes appears to be exclusive to Marseilleviridae. By contrast, histone homologs in Phycodnaviridae, Pandoraviridae (also giant viruses), and Polydnaviridae (incertae sedis) are never found as either doublets or as early-branching homologs of eukaryotic histones, suggesting recent acquisition from eukaryotes.
Unlike histones, most of the viral chromatin-associated genes exhibited a mixture of prokaryotic and eukaryotic phylogenetic affinities and often lack affinity to any specific eukaryotic gene family ( Fig. 5i and Extended Data Fig. 7). Viral readers, on the other hand, are often embedded within eukaryotic clades in gene trees and are similar to bona fide eukaryotic families, exhibiting topologies consistent with recent, secondary acquisitions. This is the case of BIRC2/3/XIAP readers widespread in the Baculoviridae, which encode BIR domains that are often hijacked from their hosts 75 . We also find a number of viral Chromo-encoding genes, which fall in two main taxonomic categories: (i) giant virus homologs of the eukaryotic CBX1/3/5 family (present in Mimiviridae, Iridoviridae and Phycodnaviridae); and (ii) homologs from various Adintoviridae, which are closely related to animal Chromo genes encoding rve integrase domains 76 (Fig. 5i). Finally, we also identify a handful of eukaryotic-like viral genes with deep-branching positions relative to core eukaryotic gene families, as seen in histones (Fig. 5h). This includes Mimiviridae homologs of the eukaryotic methyltransferases SMYD1-5 and DOT1 (Extended Data Fig. 7d,e), as well as SNF remodeller families with homologs in distinct giant virus clades (HLTF/TTF2 in Phycodnaviridae, Mimiviridae and Iridoviridae). These results indicate that cases of horizontal transfer from eukaryotes to viruses are common in different chromatin-related gene families, including histones. Therefore, it is likely that basally-branching giant virus histones were similarly acquired from a stem eukaryotic lineage and this would explain the observed histone tree topology with extant eukaryotic species. In any case, most of the eukaryotic chromatin machinery appears to have cellular roots.

Discussion
Our comparative proteogenomics study reconstructs in detail the origin and evolutionary diversification of eukaryotic chromatin components, from post-translational modifications to gene family domain architectures. We looked first at the pre-eukaryotic roots of chromatin. Multiple aspects of archaeal chromatin have been studied in recent years, including nucleosomal patterns 31 and the structure of the archaeal nucleosome 30 . A recent taxonomic survey of archaeal nucleoid-associated proteins revealed multiple independent diversifications of DNA-wrapping proteins and a strong association between high levels of chromatinization and growth temperature, overall suggesting a structural, non-regulatory role for archaeal chromatin 77 . Our proteomics data support this notion by showing the scarcity of hPTMs in four species belonging to two different archaeal lineages (Euryarchaeota and Thaumarchaeota). An earlier proteomics study reported the complete absence of hPTMs in the euryarchaeote Methanococcus jannaschii 34 . Here we do identify a few instances of modified lysine residues in Euryarchaeota, which is in line with the recently reported acetylations in Thermococcus gammatolerans histones 78 . It remains to be seen if hPTMs are frequently present in Asgard and other unsampled archaeal linages, where other eukaryotic-like features have been found 57,79,80 . In fact, some of these Asgard, particularly Lokiarchaeota, encode for histones with long, K-rich N-terminal tails but that bear no similarity with eukaryotic histones and are, therefore, most probably the result of convergent evolution. Interestingly, Lokiarchaeota genomes also frequently encode histone modifiers such as SET methyltransferases and MOZ_SAS acetylases. However, overall our results suggest that extensive usage of hPTMs is an eukaryotic innovation (Fig. 6a). Similarly, while we find the majority of catalytic domains of hPTM writers, hPTM erasers and chromatin remodellers in Archaea and even Bacteria, these appear only scattered in a small fraction of the examined taxa. In contrast, hPTM reader domains and histone chaperones are eukaryotic innovations, further supporting the idea that the functional readout of hPTMs and the role for histone variants in defining chromatin states are both exclusive to eukaryotes (Fig. 6a).
The origin of eukaryotes represents a major evolutionary transition in the history of life 81 . Thanks to sequencing and comparative analysis of archaeal and eukaryotic genomes, we also have a detailed reconstruction of the massive innovation in gene repertoires that occurred at the origin of eukaryotes. This gene innovation in the Last Eukaryotic Common Ancestor (LECA) includes cytoskeletal proteins and associated motors like myosins 82,83 and kinesins 84 , vesicle trafficking apparatus 85 , splicing machinery 86 , ubiquitin signalling systems 87 and a large repertoire of sequence-specific transcription factors 37 . Combining parsimony analysis and knowledge on gene function in extant lineages (mostly vertebrates, yeast and plants), our results allow us to reconstruct a complex LECA repertoire of hPTMs and associated writing, eraser and reader gene families (Fig. 6b,c). We infer 23 to 29 highly-conserved lysine acetylations in canonical histones (e.g., H3K9ac and H3K27ac) and a repertoire of 65 and 20 histone acetylase and deacetylase families, respectively. With the exception of H4K16ac 64 , most histone acetylations are thought to exert a generic, perhaps additive, effect on the opening of chromatin 22 . As such, acetylation marks like H3K27ac have been found to be enriched in promoters of active genes in diverse eukaryotes 42 . In contrast, histone methylations often have very specific readouts and they can be linked both to active and repressive chromatin states. We infer between 13 and 25 conserved methylated lysine residues in LECA histones, including marks typically associated to active promoters (H3K4me1/me2/me3), gene bodies (H3K36me3, H3K79me1/2, H4K20me1), and repressive chromatin states (H3K9me2/me3, H3K27me3, H4K20me3) 88,89 . Finally, we also infer the existence of five histone variants in the LECA (cenH3, H3.3, H2A.Z, macroH2A and H2A.X), as well 33 chromatin remodellers (e.g., EP400/SWR1 and INO80, involved in loading and removal of H2A.Z, respectively) and 25 histone chaperones (e.g., ASF1A/B and NPM1/2/3). This indicates that, in addition to an extensive repertoire of hPTMs, the regulation of nucleosomal histone composition was also an important feature in the LECA.
Chromatin evolution after the origin of eukaryotes is characterized by an expansion of lineage-specific histone variants harboring unique hPTMs and a net expansion in the number of reader gene families, as opposed to the relatively static catalytic gene families (writers, erasers and remodellers). This is particularly relevant as it suggests extensive remodelling of chromatin networks during eukaryote evolution, that is, changes in the coupling of particular hPTMs to specific functional chromatin states. An example of such changing state-definitions comes from looking at the hPTMs associated to TEs in different organisms: H3K9me3+H4K20me3 in animals, H3K27me3 in some plants 90 , H3K79me2+H4K20me3 in the brown multicellular algae Ectocarpus siliculosus 43 , and H3K9me3+H3K27me3 in the ciliate Paramecium tetraurelia 91 . In the context of the histone code hypothesis 3,20,92-94 , our findings indicate that, while there is an ancient core of conserved hPTMs across eukaryotes, evidence for a universal code/functional-readout is limited, with perhaps the exception of the highly conserved configuration of ancient hPTMs around active promoters across many eukaryotes 42 . Another interesting observation related to the evolution of chromatin networks is the capture of chromatin reader domains by TEs. We find evidence of this phenomenon in a number of species with a scattered phylogenetic distribution, suggesting that it is a recurrent process and that it often leads to the successful propagation of the TE in the host genome. We hypothesize that this process facilitates the targeting of TEs to specific chromatin states, as it has been described in the case of MBD DNA methylation readers captured by TEs 95,96 .
In the future, a broader phylogenetic understanding of the genome-wide distribution of hPTMs, as well as the direct interrogation of hPTM binders in different species [97][98][99] , will be crucial to further clarify questions such as the ancestral role of specific hPTM and the co-option of ancient hPTMs into novel functions.

Europe PMC Funders Author Manuscripts
Europe PMC Funders Author Manuscripts from Air Liquide GmbH, Schwechat, Austria. Nitrososphaera viennensis EN76 was grown in continuous culture in a bioreactor as previously described 100 .
Cells were harvested via centrifugation at 21,000xg 4°C 1h (Thermo scientific, Sorvall Lynx 4000 centrifuge), the supernatant discarded and the resulting pellet resuspended in 1ml of spent medium, followed by another round of centrifugation at 21,000xg 4°C for 1h (Eppendorf, Centrifuge 5424R). Pellets were stored at -70°C. All archaeal histones were extracted as described below.

Histone acid extraction
Starting material was a pellet of 50-100M cells (washed once with cold PBS) or a flashfrozen tissue homogenate in liquid nitrogen using a ceramic mortar grinder. Cells were washed first in 10ml of buffer I (10 mM TrisHCl pH 8, 10 mM MgCl 2 , 0.4M Sucrose). After 5min incubation, samples were centrifuged at 8.000g for 20min at 4°C and supernatant was removed. The resulting pellet was resuspended in 1.5ml of Buffer II (10 mM TrisHCl pH 8, 10 mM MgCl 2 , 0.25M Sucrose, 1% Triton X-100, 1% Igepal Ca-630) and incubated 15min on ice. In specific cases, cells at this stage were broken using a 2ml Dounce homogenizer (with Pestle B) or with a 20G syringe. Then samples were centrifuged at 15.000g for 10min at 4°C and supernatant was removed. The resulting pellet was then slowly resuspended in 300μL of Buffer III (10 mM TrisHCl pH 8, 2 mM MgCl 2 , 1.7M Sucrose, 1% Triton X-100) and then resulting resuspended nuclei were layered on top of another 300μL of Buffer III. Sample was centrifuged at 20.000g for 1h at 4°C and supernatant was removed, resulting in a nuclear pellet ready for acid histone extraction. All buffers were supplemented with spermidine (1:1000), beta-mercaptoethanol (1:1000), protease inhibitors (1x cOmplete cocktail Roche #11697498001, 1mM PMSF, 1:2000 Pepstatin), phosphatase inhibitors (1x phoSTOP cocktail Roche #4906845001) and deacetylase inhibitors (10mM Sodium butyrate).
For samples processed using a high-salt + HCl extraction protocol 101,102 , the pellet was resuspended in 500μL of High Salt Extraction Buffer (20 mM TrisHCl pH 7.4, CaCl 2 1M and protease, phosphatase and deacetylase inhibitors, same as above). Sample was incubated on ice for 30min and then pure HCl has added to a final 0.3N concentration (12.82μL to the initial 500μL). Samples were incubated for at least 2h on a rotor at 4°C and then centrifuged at 16.000g for 10min at 4°C to remove cellular/nuclear debris. The resulting supernatant containing solubilized histones was transferred to a clean 1.5ml tube and Trichloroacetic Acid (TCA) was added drop-wise to 25% final concentration (171μL TCA to an approximate initial 513μL sample) and left overnight at 4°C to precipitate histones. Samples were then centrifuged at 20.000g for 30min at 4°C and the supernatant removed. The pellet was then washed twice with 500μL of cold acetone and then dried for 20min at room temperature. Finally, clean histone pellets were resuspended in 30-50μL of ultrapure water. Protein concentration in the sample was measured using BCA and extraction was examined using an SDS-PAGE protein gel with Coomassie staining.
For samples processed using H 2 SO 4 102 , the protocol was exactly the same except that 400μL 0.4N H 2 SO 4 (freshly diluted) was used instead, with a similar incubation time of at least 2h at 4°C.

Histone chemical derivatization
Histones samples were quantified by the BCA method and 10 μg of each sample were derivatized with propionic anhydride, digested with trypsin and derivatized again with phenylisocyanate as previously described 49 . Briefly, samples were dissolved in 9 μL of H2O and 1 μL of triethyl ammonium bicarbonate was added to bring the pH to 8.5. The propionic anhydride was prepared by adding 1 μL of propionic anhydride to 99 μL of H2O and 1 μL of propionic anhydride solution was added immediately to the samples with vortexing and incubation for 2 minutes. The reaction was quenched with 1 μL of 80mM hydroxylamine and samples were incubated at room temperature for 20 minutes. Tryptic digestion was performed for 3 h with 0.1 μg trypsin (Promega Sequencing Grade; Madison, WI) per sample. A 1% v/v solution of phenyl isocyanate (PIC) in acetonitrile was freshly prepared and 3 μl added to each sample (17 mM final concentration) and incubated for 60 min at 37 °C. Samples were acidified by adding 50 μL of 5% formic acid, vacuum dried and desalted with C18 ultramicrospin columns (The Nest Group, Inc, Southborough, MA).

Liquid Chromatography-Tandem Mass Spectrometry Sample Acquisition
A 2-μg aliquot of the peptide mixture was analyzed using a LTQ-Orbitrap Fusion Lumos mass spec-trometer (Thermo Fisher Scientific, San Jose, CA) coupled to an EASY-nLC 1000 (Thermo Fisher Scientific, San Jose, CA) with both collision induced dissociation (CID) and high energy collision dissociation (HCD) fragmentation.
Peptides were loaded directly onto the analytical column and were separated by reversedphase chromatography using a 50-cm column with an inner diameter of 75 μm, packed with 2 μm C18 particles spectrometer (Thermo Scientific, San Jose, CA, USA) with a 90 min chromatographic gradient. The mass spectrometer was operated in positive ionization mode using a data dependent acquisition method. The "Top Speed" acquisition algorithm determined the number of selected precursor ions for fragmentation.

Mass-spectrometry Data Analysis
Acquired data were analyzed using the Proteome Discoverer software suite (v2.0, Thermo Fisher Scientific), and the Mascot search engine (v2.6, Matrix Science 103 ) was used for peptide identification using a double-search strategy. First, data were searched against each organism protein database plus the most common contaminants considering Propionylation on N-terminal, Propionylation on Lysines and Phenylisocyanate on N-terminal as variable modifications. Then a new database was generated with the proteins identified in the first search,, and a second search was done considering Propionylation on N-terminal, Propionylation on Lysines, Phenylisocyanate on N-terminal, Dimethyl lysine, trimethyl lysine, propionyl + methyl lysine, acetyl lysine, crotonyl lysine as variable modifications. Precursor ion mass tolerance of 7 ppm at the MS1 level was used, and up to 5 missed cleavages for trypsin were allowed. False discovery rate (FDR) in peptide identification was set to a maximum of 5%. The identified peptides were filtered by mascot ion score higher than 20 and only PTMs with a localization score ptmRS 104 higher than 45 were considered. The raw proteomics data have been deposited to the PRIDE 105 repository with the dataset identifier PXD031991.

Analysis of hPTM conservation
Identification of canonical and variant histones-We classified histone protein domains from a database of eukaryotic, prokaryotic and viral sequences (see details below) according to their similarity to known canonical (H2A, H2B, H3, H4) and variant histones (e.g., H2A.Z, macroH2A, cenH3 or H3.3), as well as other gene families with histone-like protein folds (e.g., the transcription factors DR1, DRAP1, NFYB/C, POLE3/4, SOS, TAF, or CHRAC). To that end, we used diamond to perform local alignments of each histone domain against (i) a set of curated histone variants obtained from HistoneDB 2.0 48 , and (ii) annotated each domain according to the best hit in the reference database, which allowed us to classify histone fold-containing proteins as canonical histones (H2A, H2B, H3, H4) or their main variants (H2A.Z, macroH2A and cenH3). This best-hit strategy performs well in distinguishing canonical histones from each other, as well as each canonical histone from its main variants (H3 from cenH3, and H2A from H2A.Z and macroH2A; Extended Data Fig.  1a).
Then, we built a graph of pairwise similarity between histones, with edges weighted by the alignment bitscore (discarding edges with bitscore < 20). We created visualisations of each connected component in this graph using the spring layout algorithm implemented in the networkx 2.4 Python library (100 iterations, weighted by alignment bitscore) 106 .
We selected the four connected components in the graph that matched the four canonical eukaryotic histones (H2A, H2B, H3, H4; discarding edges with bitscore < 20), retrieved the protein sequences for each of them, aligned them using mafft (E-INS-i mode, 1,000 iterations) 107 , and built phylogenetic trees with IQ-TREE 2.1.0 (-fast mode) 108 .
Identification of hPTM homology-We retrieved the protein sequences of the canonical histones identified in each of the 26 species and we used them for the proteomic analysis of hPTMs, and aligned them using mafft (G-INS-i mode, up to 10,000 refinement iterations).
For this subset of species, histone class identity was cross-referenced with the HistoneDB search tool. Then, we manually aligned the peptides mapping onto these proteins to identify the position of each hPTM along a consensus alignment. In the case of H3, H4, and macroH2A, the majority of alignment positions were conserved across most eukaryotes in our dataset, and we used a consensus numbering scheme. In the case of H2A, H2A.Z, and H2B, non-conserved insertions and deletions at the N-terminal tail precluded the use of a paneukaryotic numbering scheme. Instead, we reported hPTM positions based on the human homolog (if possible), or relative to taxonomically restricted conserved positions. In cases where position-wise homology could not be established, we grouped multiple amino-acids into stretches of unclear homology, which we report separately from conserved positions (question mark symbols in Fig. 1). The complete list of hPTMs and their position-wise coordinates relative to the consensus alignment is available in Supplementary Data 3.
Furthermore, we also reported the presence (in any position) of modifications in lessconserved histone variants, as well as the linker histone H1.

Comparative genomics analysis of chromatin-associated proteins
Data retrieval-We identified homologs of gene families associated with eukaryotic chromatin, using a database of predicted proteomes from a selection of eukaryotic species from all major supergroups (n = 172 species; see Supplementary Data 1 for their taxonomic classification and data sources), as well as archaeal and viral peptides available in the NCBI nr peptide collection (as of 25th of April, 2020) and bacterial peptides available in RefSeq Orthology identification-We aimed to identify groups of orthologs within each of the 61 chromatin-associated gene families using targeted phylogenetic analyses. We followed the following strategy for each of the 59 sets of eukaryotic genes. First, we partitioned each set into one or more homology groups based on pairwise local sequence alignments using diamond 0.9.36.137 (high sensitivity all-to-all search) 124 , followed by clustering of the resulting pairwise alignments graph with MCL 14.137 (--abc mode) 125  Then, we parsed the species composition of each gene tree in order to identify groups of orthologous proteins using the POSSVM pipeline 128 . Specifically, we used the species overlap algorithm 129 implemented in the ETE toolkit 3.1.1 130 , which identifies pairs of orthologous genes in a phylogenetic tree by examining the species composition of each subtree, and classifying internal nodes as paralogy nodes (if there is overlap in the species composition between each of its two descendant subtrees) or orthology nodes (if there is no overlap). Pairs of genes linked by an orthology node are then recorded as orthology pairs. In our analysis, we used an overlap threshold=0 (i.e. any species composition overlap between the two descendant subtrees is classified as a paralogy event). The resulting list of pairwise orthology relationships between genes was clustered into groups of orthologs (orthogroups) using MCL. We further annotated each orthogroup with a string denoting the gene symbols of the human proteins therein (if any).
Overall, we classified 51,426 proteins from 61 gene classes (defined by protein structural domains), divided into 242 gene trees and 1,713 gene families (orthogroups). The source peptide sequences and gene trees used for these analyses are available in Supplementary Data 7 and 8.
Ancestral reconstruction of gene content-We inferred the presence, gain and loss of each orthogroup along the eukaryotic tree of life, using a phylogenetic birth-and-death model 131 implemented in Count 132 . This tool takes a numeric profile of gene family presence/absence in extant species (172 in our dataset) and a phylogenetic tree defining their evolutionary relationships, and infers the probabilities of gain and loss of each family at each ancestral node along the tree.
First we trained the probabilistic model in Count. As a training set, we used a random sample of 1,000 PFAM domains annotated in the 172 species of interest (restricting the sampling to domains present in at least 5% of species). The final model consists of gain, loss and transfer rates with two Γ categories each, and a constant duplication rate (given that we only recorded gene presence/absence, duplication events are not included in our downstream analyses). This model was obtained in three sequential rounds of training, so as to sequentially add zero, one and two Γ categories to each evolutionary rate. Each round consisted of up to 100 iterations, and stopped when the relative change in the model log-likelihood fell by 0.1% in two consecutive rounds. The final evolutionary rates and the Newick-formatted species tree used in this step are available in the Supplementary Data 1 and Extended Data Fig. 3a.
Second, we calculated the posterior probability of gain, loss and presence of each orthogroup in our dataset with Count. The aggregated counts of gains and losses of the various classes of chromatin-associated proteins (acetylases, deacetylases, methyltransferases, demethylases, readers and remodellers) along the eukaryotic tree were obtained by summing the probabilities of gain, presence or loss of all orthogroups of a given class at each ancestral node. To investigate the evolutionary histories of specific orthogroups at a given node in the tree, we applied a probability threshold of 0.9 (for presence) or 0.5 (to identify the most probable gain and loss node). The Count model was not able to calculate ancestral probabilities for a few orthogroups with widespread phylogenetic distributions, due to violations of the birth-and-death model (25 out of 1,713 families). In order to be able to report presence probabilities in the LECA for these orthogroups, we inferred their presence in this ancestor using the Wagner parsimony procedure implemented in Count with a gain-to-loss penalty g = 5, and recorded their presence as binary values (0/1) accordingly.
Protein domain architecture analyses-We annotated the Pfam domains present in each protein from the gene classes listed in Supplementary Data 4, using Pfamscan 1.6-3 and the Pfam 33.0 database 122 . We visualized the networks of protein domain co-occurrence from the point of view of the core domain(s) that define each gene class, using the networkx Python library (version 2.4) 106 . Specifically, we built a graph where each node represented 'accessory' domains (i.e. domains that co-occur with the 'core' domain that defines given gene class), node size reflected number of co-occurrences with the 'core' domain, and edges reflected co-occurrences between accessory domains. We identified communities of frequently co-occurring accessory domains using the label propagation algorithm implemented in networkx (communities submodule), which we used as a basis to manually annotate groups of co-occurring domains of interest (Fig. 5c). Network visualizations were created using the NEATO spring layout algorithm from the Graphviz 2.40.1 Python library 133 .
In parallel, we also recorded the presence of Pfam domains within individual orthogroups, and their taxonomic distribution. The archaeal and bacterial protein sets were also reduced with CD-HIT (at 95% and 90% sequence similarity, respectively). Each set of sequences was then partitioned into low-granularity homology clusters using the MCL-based strategy described above (inflation I = 1.2), and a phylogenetic tree was then constructed from each homology cluster with IQ-TREE (as described above).

Prokaryotic roots of the eukaryotic chromatin machinery-We
Then, we mapped each eukaryotic gene to its orthogroup (obtained from eukaryoticonly analyses, see above) and used the distribution of phylogenetic distances from the prokaryotic+eukaryotic gene trees to classify them according to their similarity to (i) eukaryotic genes in other orthogroups, (ii) archaeal homologs, or (iii) bacterial homologs. Specifically, we used a majority-voting procedure in which we recorded the number of sequences of eukaryotic, archaeal or bacterial origin amongst the ten nearest neighbors of each gene (measuring intergenic distances as substitutions per site), and assigned the most common taxonomic group as the 'closest' homolog of that gene (minimum 50% agreement). This fraction is termed 'Phylogenetic affinity score' and reported in Supplementary Data 5. The pairwise distances were obtained from each gene tree using the cophenetic distance method in the cophenetic.phylo utility of the ape 5.4 R library 135 .
domains (version 33.0), using a list compiled from 68,136 (complete list in Supplementary Data 4), totaling 823 candidate fusions from 91 species (listed in Supplementary Data 6). We annotated these genes to their most similar known TE element by aligning them against the Dfam 3.3 database 137 using the tblastn program in BLAST 2.2.31 138 .
We validated each candidate fusion using the following criteria: (i) contiguity of the gene model on the genome assembly, i.e., recording which genes were interrupted by poly-N stretches (which might indicate an incorrect gene model); (ii) evidence of expression in at least one sample from a range of publicly available transcriptomic experiments (from the NCBI SRA repository); (iii) evidence of contiguous expression, i.e., whether an expressed transcript had mapped reads along the entire region located between the 'core' and 'TEassociated' domains; (iv) we also recorded the number of exons per gene; and (v) located near any other candidate fusion gene in the genome. Analysis of viral homologs-We investigated the homology of the viral chromatinassociated genes (which included 19 out of 61 families present in our survey) using joint phylogenetic analyses of protein domains from virus, prokaryotic and eukaryotic genes. We used the same method described above to investigate the prokaryotic roots of eukaryotic gene classes: we aligned viral domains against a database of cellular homologs (high sensitivity diamond search), followed by low-granularity MCL clustering (inflation I = 1.2) and phylogenetic tree building (IQ-TREE). Then, we used the same majority-voting procedure described above to classify viral homologs according to their similarity to eukaryotic, archaeal or bacterial gene families based on their distribution of phylogenetic distances. For viral genes that were most similar to eukaryotic genes, we used the same procedure to map them to their closest eukaryotic orthogroup.
The complete list of viral genes and their phylogenetic annotation is available in Supplementary Data 6. Out of 2,163 viral genes in our dataset, 2,144 could be annotated as similar to a particular cellular group using this procedure (99.1%), and the majority of these genes had a high agreement in the annotations of their nearest neighbors (2,096 with ≥50% agreement; 1,449 with ≥90% agreement).
In the case of viral histones, we built a separate phylogeny with a few modifications in our protocol: (i) we used additional viral genes obtained from 71 as a reference; (ii) we omitted the CD-HIT reduction and MCL partitioning steps, and jointly analyzed the entire set of homologs instead; and (iii) in the phylogenetic reconstruction step, we used the approximate Bayes posterior probabilities 141 implemented in IQ-TREE.  Fig. 1d). Alignments were plotted using the msa 1.24.0 library in R 142 . than one canonical histone or major variant (macroH2A, H2A.Z or cenH3), the scatter plots represent the relative identity between the primary (horizontal axis) and secondary alignment(s) (vertical axis). b, Aggregated counts of histone gene pairs, classified according to histone type and orientation. c, Presence of histone variants (left) and number of collinear pairs of histone-encoding genes (right) per species, classified according to their histone types and relative orientation (head-to-head, hh; head-to-tail, ht; and tail-to-tail, tt). Source data available in Supplementary Data 2. Histone variant classification is based on the highest-scoring HMM profile from HistoneDB. Asterisks colors in the macroH2A column indicate species where histone-less Macro domains orthologous to the macroH2A genes are found (see panel d). Lighter colors in the variant classification indicate ambiguously classified histones (i.e. cases in which the highest-scoring HMM profile exhibited a low bitscore, defined as a probability below 0.05 in the profile-wise distribution function of scaled bitscores; or cases in which the first-to-second ratio between high scoring profiles was below 1.01). d, Alignments of putatively conserved histone N-tails in archaea.

Identification of archaeal N-terminal histone tails-We
Conserved amino-acids are color-coded according to chemical properties. Dots next to species names are color-coded according to taxonomy (same as Fig. 2c). e, Phylogenetic analysis of the Macro motif of macroH2A histones across eukaryotes, highlighting the macroH2A ortholog group (green), and, within this group, Macro-containing genes lacking histone domains (orange), and their protein domain architectures. In each heatmap, we indicate the fraction of genes within an orthogroup (rows) that contain a specific protein domain (columns). Domains in bold are catalytic (black) or reader (purple) functions. At the right of each heat-map, we summarize the presence/absence profile of each orthogroup across eukaryotic lineages (as listed in Fig. 1a). Fig. 5. Evolution of the hPTM reader toolkit. a, Pie plot representing the number of genes classified as part of the catalytic (acetylases, deacetylases, methyltransferases, demethylases, remodellers or chaperones) or reader families, or as both. The barplot at the right shows the most common reader domains in genes classified with both reader and catalytic functions. b, Pie plot representing the number of reader domain-encoding genes classified according to whether they contain one type of reader domain (e.g., PHD) or more than one (e.g., PHD + PWWP). The barplot at the right shows the most common combinations of reader domains among genes with multiple reader domains. c, Summary of gene family gains per reader family, with example cases highlighted in selected nodes. Node size is proportional to number of gains at 90% probability.

Extended Data
Extended Data Fig. 6

Supplementary Material
Refer to Web version on PubMed Central for supplementary material.    a, Summary of the seven classes of genes with chromatin-related activity covered in our survey: histone-specific hPTM writers (acetylases and methyltransferases), erasers (deacetylases and demethylases), readers, remodellers, and chaperones. b, Percentage of surveyed taxa containing homologs from each chromatin-associated gene class, for eukaryotes (top), archaea, bacteria, and viruses (bottom). Species-level tables are available in Extended Data Fig. 3. c, Number of eukaryotic genes classified in each of the chromatinassociated modification enzymes, readers, remodellers, and chaperones. d, Overlap between the taxon-level phylogenetic distribution of histones and chromatin-associated domains in a, Summary of phylogenetic affinities of the eukaryotic homologs of gene classes that are also present in prokaryotes. For each gene family, we evaluate whether it is phylogenetically closer to a majority (≥50%) of eukaryotic sequences from a different orthogroup (indicating intra-eukaryotic diversification), or to sequences from Bacteria or Archaea. b, Left, gene tree of eukaryotic and prokaryotic Sirtuin deacetylases, showcasing an example of a eukaryotic family that diversified within eukaryotes (SIRT6) and another one with close relatives in Asgard archaea (SIRT7). Right, gene tree of KAT14 acetylase, a eukaryotic orthogroup with bacterial origins. Statistical supports (UF bootstrap) are shown at selected internal nodes of the highlighted clades. c, Evolutionary reconstruction of hPTM writer and eraser gene families, remodellers, and histone chaperones along the eukaryotic phylogeny, including the number of genes present in the last eukaryotic common ancestor (LECA). Barplots indicate the number of orthologs of each gene family present at the LECA (at 90% posterior probability; see Methods) and whether the presence of a given orthogroup at LECA is supported by its conservation in various early-branching eukaryotic lineages (Amorphea, Discoba, Diaphoretickes and others). The list of ancestral gene families below each plot is non-exhaustive. Two ancestral gene counts are provided: all families at presence probability above 90%, and, in brackets, the subset of these that is present in at least two of the main eukaryotic early-branching lineages (Amorphea, Diaphoretickes, and Discoba). Source data in Supplementary Data 5. d-e, Reconstructed evolutionary origins of the different subunits of the Polycomb repressive complexes (PRC2 and PRC1) and Trithorax-group complexes (KMT1 to 5). f-h, Side-by-side comparison of the presence of individual hPTM marks and various subunits of the Polycomb and Trithorax complexes, as well as other hPTM writers, responsible for their deposition. a, Evolutionary reconstruction of reader gene families along the eukaryotic phylogeny, highlighting the number of gains along the eukaryotic phylogeny (at 90% posterior probability). The Euler diagram at the top shows the overlap between presence of chromatinassociated catalytic domains and readers. The barplot at the left indicates the number of orthologs of each gene family present at the LECA and whether their presence is supported by its conservation in various early-branching eukaryotic lineages (Amorphea, Discoba, Diaphoretickes, and others). Pie plots at the right summarize the number of orthogroups from each gene family gained within selected lineages: Metazoa, Holomycota, Viridiplantae and SAR+Haptophyta. b, Number of reader or catalytic orthogroups gained at each node in the species tree, for selected species. Source data in Supplementary Data 5. c, Networks of protein domain co-occurrence for Chromo and PHD readers. Each node represents a protein domain that co-occurs with Chromo or PHD domains, and node size denotes the number of co-occurrences with either Chromo or PHD. Edges represent co-occurrences between domains. Groups of frequently co-occurring protein domains have been manually annotated and color-coded, which has revealed sub-sets of retrotransposon and DNA transposon-associated domains.  a, Summary of events in chromatin evolution prior to, during and after the origin of eukaryotes. b, Number of chromatin-related gene families and hPTM marks inferred to have been present at the LECA. Ancestral gene counts are indicated at >90% probability. For gene counts, numbers within bars indicate the subset of families present in at least two of the most deeply-sampled early-branching eukaryotic lineages (Amoropha, Diaphoretickes, and Discoba). For hPTMs, the ancestral counts have been inferred using Dollo parsimony assuming a Diaphoratickes -Amorphea split at the root of eukaryotes, and numbers within bars indicate the number of hPTMs whose ancestral presence is supported by more than one species at both sides of the root. c, hPTMs inferred to be present in the last eukaryotic common ancestor (LECA) based on Dollo parsimony. Only amino-acid positions conserved in all eukaryotes in our dataset are shown. Asterisks indicate modifications whose presence at the LECA is supported by just one species at either side of the root. The inferred LECA presence of known writing/erasing enzymes associated to these hPTM is indicated.