dataesr/scientific-paragraphs-categorization huggingface dataset:
Results Comparative proteomics of eukaryotic histone modifications. We analysed the phylogenetic distribution and evolutionary history of histone proteins. To this end, we surveyed the presence of histone-fold proteins across 172 eukaryotic and 4,226 archaeal taxa, using HMM searches (Fig. 1a ,b and Supplementary Data 1). Histone proteins are found in all eukaryotic genomes. We clustered the identified 8,576 histone-encoding proteins using pairwise local alignments and then classified individual sequences in these clusters on the basis of pairwise alignments to a reference database48 (Fig. 1a  and Extended Data Fig. 1a ). This reveals four broad clusters corresponding to the four main eukaryotic histones (H2A, H2B, H3 and H4) and their variants (H2A. Z, macroH2A and cenH3), as well as a fifth cluster composed of archaeal HMfB homologues.
We sought to obtain a systematic, phylogenetics-based classification of histone remodellers, chaperones, readers and modifiers to understand the evolutionary history of eukaryotic chromatin (Fig. 3a ). To this end, we (1) 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; (2) 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 (3) 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.
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 histones12. This analysis also reveals a substantial fraction of eukaryotic gene families with close orthogroups in archaea and bacteria, which pinpoints components that were (1) inherited from a prokaryotic ancestor during eukaryogenesis; (2) laterally transferred between eukaryotes and prokaryotes at later stages; or (3) a combination of both phenomena.
Thus, we decided to explore in detail the occurrence of chromatin-associated 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 (1) come from a nondiscontinuous gene model in the original assembly and (2) have evidence of expression, with reads mapping along the entire region between the TE-associated domain and the chromatinassociated domain (Extended Data Fig. 6 ). 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 methyltransferase and a Mariner-class DNA transposon69. 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 lysine (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.
Peptides were loaded directly onto the analytical column and were separated by reversed-phase chromatography using a 50-cm column with an inner diameter of 75 μm, packed with 2 μm C18 particles spectrometer (Thermo Scientific) 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 analysed using the Proteome Discoverer software suite (v. 2. 0, Thermo Fisher Scientific) and the Mascot search engine (v. 2. 6, Matrix Science103) 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 phenyl isocyanate on Nterminal 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 Nterminal, propionylation on lysines, phenyl isocyanate on N-terminal, dimethyl lysine, trimethyl lysine, propionyl + methyl lysine, acetyl lysine and crotonyl lysine as variable modifications. Precursor ion mass tolerance of 7 ppm at the MS1 level was used and up to five missed cleavages for trypsin were allowed. False discovery rate in peptide identification was set to a maximum of 5%. The identified peptides were filtered by mascot ion score >20 and only PTMs with a localization score ptmRS104 >45 were considered. The raw proteomics data have been deposited to the PRIDE105 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 and H4) and variant histones (for example, H2A. Z, macroH2A, cenH3 or H3. 3), as well as other gene families with histone-like protein folds (for example, 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: (1) a set of curated histone variants obtained from HistoneDB 2. 0 (ref. 48 ) and (2) 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 and 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 visualizations 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 and 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.
siliculosus43, the diatom Phaeodactylum tricornutum109, the ciliate T. thermophila46,110-112, the ascomycete N. crassa113, S. cerevisiae and Schizosaccharomyces pombe46 and the plant A. thaliana114-116. When available in public repositories, we reanalysed these datasets using the strategy described above. Finally, we also complemented our own proteomics data using previously published hPTM data from H. sapiens46,117-120 and C. owczarzaki42. Comparative genomics analysis of chromatin-associated proteins. Data retrieval. We identified homologues 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 non-redundant peptide collection (as of 25 April 2020) and bacterial peptides available in RefSeq (release 99, 11 May 2020).
The database of viral sequences was complemented with peptides from 501 genomes of nucleocytoplasmic large DNA viruses121. Gene family searches. We defined 61 gene classes associated with eukaryotic chromatin, based on HMM models obtained from the Pfam database (release 33. 0)122. This list included canonical and linker histones (n = 2 families), chromatin-specific lysine acetylases (n = 5), deacetylases (n = 2), methyltransferases (n = 2), demethylases (n = 2), chromatin readers (n = 16), remodellers (n = 1) and chaperones (n = 13), as well as multiple families associated with the Polycomb complexes (n = 18). The complete list of gene families, including the associated HMM models, is available in Supplementary Data 4. For each gene family, we retrieved all homologues from the eukaryotic, archaeal, bacterial and viral databases using the hmmsearch tool from the HMMER 3.
3 toolkit123 and the gathering threshold defined in each Pfam HMM model. We recorded the taxonomic profile of each homologue. Orthology identification. We aimed to identify groups of orthologues 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 on the basis of 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, using low inflation values (Supplementary Data 4) to favour inclusive groupings. Second, we performed multiple sequence alignments of each homology group with mafft 7.
471 (ref. 107)  under the E-INS-i mode (optimized for multiple conserved regions), running up to 10,000 refinement iterations. Third, we trimmed the resulting multiple sequence alignments using clip-kit 0. 1 (kpic-gappy mode)126. Fourth, we built phylogenetic trees for each trimmed alignment using IQ-TREE 2. 1. 0 (ref. 108) , selecting the best-fitting evolutionary model using its ModelTest module (according to the Bayesian Information Criterion) and using 1,000 UFBS bootstrap supports127. Each tree was run for up to 10,000 iterations until convergence was attained (at the 0. 999 correlation coefficient threshold and for at least 200 iterations). Then, we parsed the species composition of each gene tree to identify groups of orthologous proteins using the POSSVM pipeline128. Specifically, we used the species overlap algorithm129 implemented in the ETE toolkit 3.
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 model131 implemented in Count132. 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 chromatinassociated 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).
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 database122. 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 (v. 2. 4)106. Specifically, we built a graph where each node represented 'accessory' domains (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 library133. In parallel, we also recorded the presence of Pfam domains within individual orthogroups and their taxonomic distribution. Prokaryotic roots of the eukaryotic chromatin machinery. We retrieved all eukaryotic domains from gene class shared with prokaryotes (Histones, Acetyltransf_1, GNAT_acetyltr_2, MOZ_SAS, Hist_deacetyl, SIR2, DOT1, SET, CupinJmjC, ING, MBT, PWWP and SNF2_N), collapsing identical sequences at 100% similarity with CD-HIT 4.
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 neighbours of each gene (measuring intergenic distances as substitutions per site) and assigned the most common taxonomic group as the 'closest' homologue 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 library135. Characterization of fusions with transposon-associated domains. We retrieved all classified genes from our eukaryotic dataset that contained transposon-associated Pfam domains (v. 33. 0), using a list compiled from refs.
68,136 (complete list in Supplementary Data 4), totalling 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 database137 using the tblastn program in  BLAST 2. 2. 31 (ref. 138). We validated each candidate fusion using the following criteria: (1) contiguity of the gene model on the genome assembly, that is, recording which genes were interrupted by poly-N stretches (which might indicate an incorrect gene model); (2) evidence of expression in at least one sample from a range of publicly available transcriptomic experiments (from the NCBI SRA repository); (3) evidence of contiguous expression, that is, whether an expressed transcript had mapped reads along the entire region located between the 'core' and 'TE-associated' domains; (4) we also recorded the number of exons per gene; and (5) located near any other candidate fusion gene in the genome.
The list of SRA experiments used for these validation steps is available in Supplementary Data 1. This list includes 64 out of 91 species for which transcriptomics datasets are publicly available and covers 768 out of the 822 TE fusion candidates (93%). RNA-seq read mapping was performed with bwa mem 0. 7. 17-r1188 (ref. 139 ) using the complete set of spliced transcripts of each species as the reference database. We used bedtools 2. 29. 2 (ref. 140)  to identify poly-N stretches in the genome assembly (assembly contiguity criterion). We identified regions of low coverage along the transcript sequence (expression contiguity criterion) using the bedtools genomecov utility, requiring that the coverage along both domains involved in each fusion and their intermediate regions be higher or equal to two reads.
Analysis of viral homologues. We investigated the homology of the viral chromatin-associated 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 homologues (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 homologues according to their similarity to eukaryotic, archaeal or bacterial gene families on the basis of 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 most of these genes had a high agreement in the annotations of their nearest neighbours (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: (1) we used additional viral genes obtained from ref. 71 as a reference; (2) we omitted the CD-HIT reduction and MCL partitioning steps and jointly analysed the entire set of homologues instead; and (3) in the phylogenetic reconstruction step, we used the approximate Bayes posterior probabilities141 implemented in IQ-TREE.
Acknowledgements We thank  A. de Mendoza  for critical input on the analysis of TE fusions. We also thank  J. Casacuberta  for P. patens samples,  H. J. G. Meijer  for P. infestans samples,  M. Adamska  for S. ciliatum samples and  A. Simpson  for access to the G. okellyi culture (made possible by his funding from  NSERC, Canada). Research in the A. S. -P. group was supported by the  European Research Council (ERC)  under the European Union's  Horizon 2020 Research and Innovation Programme  (grant agreement no. 851647 ) and the  Spanish Ministry of Science and Innovation  ( PGC2018-098210-A-I00 ). We also acknowledge support of the  Spanish Ministry of Science and Innovation  to the EMBL partnership, the  Centro de Excelencia Severo Ochoa  and the  CERCA Programme (Generalitat de Catalunya ).
Petzold et al (2017) , showed that Lp isolates, sharing a unique ST, that fell into two separate clusters, were only different because of recombination events in specific regions of the genome. In Bresso, Italy, two outbreaks associated with ST23 were identified in 2014 (6 cases) and in 2018 (52 cases). The only suspected, yet unlikely, source of infection was identified as a public fountain, where Lp1 Philadelphia and Lp1 France/Allentown ST23 matching the clinical isolates was found  (Faccini et al, 2020). The SBT database was queried and ST23 was found to be widespread but in a restricted number of countries including Italy, where it caused all up to now described outbreaks and also several sporadic cases  (Rota et al, 2005 (Rota et al, , 2011;; Scaturro et al, 2015). It also represents the most frequently community-and travel-associated ST  (Fontana et al, 2014).
The highest number of cg differences was found in 299C from Verona. The ST23 isolated from Danish patients (orange nodes in Fig 3A ) differed for 1-28 loci from ST23 isolated in Italy, France, South Europe, and Switzerland. In particular, the genomes having DK02 and DK03, from a Danish patient, showed only 12 cg loci differences encoding Icm proteins of the Dot/Icm type IV secretion system, but as determined by Geneious analysis the variations were all conserved. ST23 isolated in Scotland from patients who had travelled to Southern Europe differed by between 1 and 33 cg loci from ST23 isolated in Italy and France. Generally, the loci differences found in ST23 genomes mainly corresponded to genes encoding proteins related to the carbohydrate, amino acids, and lipid metabolism, and several for transport and binding (Table  S9 ).