Benoît Aliaga

Ingo Bulla

Gabriel Mouahid

David Duval

Christoph Grunau
email: christoph.grunau@univ-perp.fr

Universality of the DNA methylation codes in Eucaryotes

Genetics and epigenetics are tightly linked heritable information classes. Question arises if epigenetics provides just a set of environment dependent instructions, or whether it is integral part of an inheritance system. We argued that in the latter case the epigenetic code should share the universality quality of the genetic code. We focused on DNA methylation. Since availability of DNA methylation data is biased towards model organisms we developed a method that uses kernel density estimations of CpG observed/expected ratios to infer DNA methylation types in any genome. We show here that our method allows for robust prediction of mosaic and full gene body methylation with a PPV of 1 and 0.87, respectively. We used this prediction to complement experimental data, and applied hierarchical clustering to identify methylation types in ~150 eucaryotic species covering different body plans, reproduction types and living conditions. Our analysis indicates that there are only four gene body methylation types. These types do not follow phylogeny (i.e. phylogenetically distant clades can have identical methylation types) but they are consistent within clades. We conclude that the gene body DNA methylation codes have universality similar to the universality of the genetic code and should consequently be considered as part of the inheritance system.

Introduction

Living organisms are biological systems in which the complex interaction between different elements such as the nuclear genotype and epigenotype factors and the environment brings about a phenotype that develops and evolves over time [START_REF] Levine | The Future of Systems Biology[END_REF][START_REF] Cosseau | Epi)genetic Inheritance in Schistosoma mansoni: A Systems Approach[END_REF] . For a complete understanding and potential control of biological processes such as development and evolution, it is therefore necessary to understand as many elements of biological systems as possible. In the present work, we focus on the epigenotype unit that we operationally define as any modification of the chromatin-DNA complex that has an impact on the expression and function of genes [START_REF] Nicoglou | Epigenetics: A way to bridge the gap between biological fields[END_REF] . Epigenetic information can be stored in a multitude of bearers such as histone modifications, non-coding RNA, the topology of the nucleus, and methylation of DNA. DNA methylation has been one of the most studied epigenetic marks since its discovery in 1948 [START_REF] Hotchkiss | The quantitative separation of purines, pyrimidines and nucleosides by paper chromatography[END_REF] . Methylation occurs at positions 4 and 5 of the pyrimidine ring of cytosine forming either 4-methyl-cytosine (4mC) or 5-methyl-cytosine (5mC), or at position 6 of the purine ring in 6-methyl-adenine (6mA). 6mA and 4mC were believed to occur only in bacteria but recent advances in sequencing technology made it possible to detect them also in eukaryotic species. A specific database (MethSMRT) was dedicated to these modifications [START_REF] Ye | MethSMRT: an integrative database for DNA N6-methyladenine and N4-methylcytosine generated by single-molecular real-time sequencing[END_REF] , and the available experimental data were used to train an algorithm to predict the occurrence of 4mC [START_REF] Chen | iDNA4mC: identifying DNA N4-methylcytosine sites based on nucleotide chemical properties[END_REF] in DNA based on sequence features. We will focus here on 5mC and to facilitate the readability use the term DNA methylation for this purpose.

In most eukaryotes, 5mC is overrepresented or even restricted to the dinucleotide CpG context, where 'p' stands for the phosphodiester linkage between the cytosine (C) and the guanine (G). In plants, the 5mC can occur in other contexts such as CpHpG or CpHpH, where 'H' stands for A, C or T (reviewed in Vanyushin 7 ). In contrast, in certain molds, methylation occurs preferentially (>60%) in CpAs [START_REF] Cambareri | Repeat-induced G-C to A-T mutations in Neurospora[END_REF] . DNA methylation is catalyzed by a family of enzymes called DNA methyltransferase (DNMT) composed of 3 canonicals members (DNMT 1, 2 and 3) [START_REF] Lyko | The DNA methyltransferase family: a versatile toolkit for epigenetic regulation[END_REF] . After replication, 5mC will be maintained by the activity of DNMT1, which has a high affinity to hemimethylated DNA, and that methylates immediately after replication the newly synthesized strand, reproducing methylation patterns in CpG dinucleotide with a fidelity of roughly 99.9% [START_REF] Riggs | Methylation dynamics, epigenetic fidelity and X chromosome structure[END_REF] thus allowing for mitotic heritability of DNA methylation patterns. The role of DNMT2 is controversial because it has little DNA methylation activity [START_REF] Hermann | The human Dnmt2 has residual DNA-(Cytosine-C5) methyltransferase activity[END_REF] but is able to methylate cytosine 38 in the anticodon loop of aspartic acid transfer RNA [START_REF] Goll | Methylation of tRNAAsp by the DNA Methyltransferase Homolog Dnmt2[END_REF] and some authors propose therefore to replace DNMT2 by tRNA (Cytosine(38)-C(5))-Methyltransferase TRDMT1 [START_REF] Albalat | Evolution of DNA-methylation machinery: DNA methyltransferases and methyl-DNA binding proteins in the amphioxus Branchiostoma floridae[END_REF] . There are species, such as the model organism Drosophila melanogaster, that have only DNMT2 and do not possess 5-methyl-cytosine in their genome, or DNA methylations is so low that it is very difficult to detect [START_REF] Hermann | The human Dnmt2 has residual DNA-(Cytosine-C5) methyltransferase activity[END_REF][START_REF] Schaefer | Solving the Dnmt2 enigma[END_REF][START_REF] Raddatz | Dnmt2-dependent methylomes lack defined DNA methylation patterns[END_REF] . These enzymes have distinct roles due to the presence of different domain structures. DNA methylation is established by DNMT3 that can methylate the two strands of the DNA i.e. has a de novo methylation function [START_REF] Okano | Cloning and characterization of a family of novel mammalian DNA (cytosine-5) methyltransferases[END_REF] . Various DNA methylation contexts are found across the plant and animal kingdoms. There are species where 5mC is present all over the genome (global methylation) while others can be entirely devoid of methylation. In species with global DNA methylation, only small regions, among them promoters and other regulatory elements, are methylation-free [START_REF] Suzuki | DNA methylation landscapes: provocative insights from epigenomics[END_REF] . If 5mC occurs in the promoters of vertebrates, it has a repressive action on the gene transcription [START_REF] Rivenbark | Epigenetic reprogramming of cancer cells via targeted DNA methylation[END_REF] .

Invertebrates often have a mosaic-type methylation pattern with high methylation in almost all CpG in large blocks of genomic DNA interspersed with almost entirely unmethylated blocks. Changes in DNA methylation occur during organ regeneration [START_REF] Casimir | Evidence for dedifferentiation and metaplasia in amphibian limb regeneration from inheritance of DNA methylation[END_REF] , aging [START_REF] Mugatroyd | The janus face of DNA methylation in aging[END_REF][START_REF] Zampieri | Reconfiguration of DNA methylation in aging[END_REF] , in response to bacterial infection [START_REF] Dowen | Widespread dynamic DNA methylation in response to biotic stress[END_REF] , as well as flowering time and root length in Arabidopsis thaliana [START_REF] Cortijo | Mapping the Epigenetic Basis of Complex Traits[END_REF] , just to name a few examples.

DNA methylation was therefore proposed as a language in which environment and genome talk to each other. Other authors have seen DNA methylation primarily as a genomic defense system against parasitic genomic elements 17 .

Given the apparent heterogeneity of DNA methylation patterns and the multitude of biological processes involved it was suggested that DNA methylation evolved in every phylogenetic clade towards a specific role in controlling gene expression. Essentially, the question is whether or not there is universality in the DNA methylation epigenetic code, conceptually similar to the universality of the genetic code, or if DNA methylation is non-universal and specific to every evolutionary unit. We wished to address this question through the analysis of the evolution of DNA methylation. Here, an obstacle is that a comprehensive analysis of DNA methylation patterns in a wide range of different species is missing. There are currently methods available that allow, in principle, for determining genome-wide DNA methylation patterns ("methylomes") at a single base resolution. Since DNA methylation patterns can be different in different organisms of a species or even tissues of an individual, for a given species several methylomes can exist. A review of the available data in different databases and in the literature showed that there is a strong bias towards model organisms: there are at least 300 methylomes available for human, mouse and the model plant to derive the general rules (if any) that govern DNA methylation in the different branches of the "tree of life". A potential solution to the caveat that experimental "wet bench" data is missing is to infer DNA methylation indirectly with computational method [START_REF] Yi | Computational approaches for understanding the evolution of DNA methylation in animals[END_REF][START_REF] Bulla | Notos -a galaxy tool to analyze CpN observed expected ratios for inferring DNA methylation types[END_REF] . The basis for this is that methylated CpG sites mutate relatively frequently compared to the other dinucleotides over evolutionary time [START_REF] Bird | DNA methylation and the frequency of CpG in animal DNA[END_REF] . If a cytosine is deamined, a deoxy-uracil will form, which is not stable in DNA: it will rapidly be excised by uracil glycosylase and replaced by cytosine. In contrast, 5mC deamination generates thymine, which is less efficiently processed by the DNA repair machinery.

Despite the existence of a specific repair mechanism that restores G/C mismatch, the mutation rate from 5mC to T is therefore 10-fold to 50-fold higher than other transitions, depending on local GC content [START_REF] Fryxell | CpG mutation rates in the human genome are highly dependent on local GC content[END_REF] . For Humans, it was estimated that within 20 years, 0.17% of 5mC in the genome were converted into thymine [START_REF] Cooper | Cytosine methylation and the fate of CpG dinucleotides in vertebrate genomes[END_REF] . If 5mC occurs predominantly in CpG pairs, the above-mentioned mechanism will increase the mutation rate from CpG to TpG or CpA and induce an underrepresentation of CpG [START_REF] Jabbari | Cytosine methylation and CpG, TpG (CpA) and TpA frequencies[END_REF] . Therefore, CpG observed/expected ratio (CpG o/e) in gene bodies can be used to predict if a species' DNA is methylated in gene bodies or not [START_REF] Bird | DNA methylation and the frequency of CpG in animal DNA[END_REF][START_REF] Razin | Distribution of 5-methylcytosine in chromatin[END_REF][START_REF] Suzuki | CpG methylation is targeted to transcription units in an invertebrate genome[END_REF] . In other words, gene bodies of a species are not methylated when the CpG o/e ratio is in average close to 1, and methylated for an average ratio far below 1. Low CpG o/e is not a condition for methylation but a consequence of it. It is important to note that this is a species-level prediction that uses methylation signatures that pass through the germ-line and need several generations of mutation accumulation to be detectable. It cannot be used to predict methylation changes of individual genes within shorter periods.

These predictions were tested in at least 13 studies comparing CpGo/e to methylation levels obtained with various methods (Table 1 and Supplementary file 1). All came to the conclusion that CpGo/e ratios correlate well (inversely) with methylation levels when species were compared.

Nevertheless, there remain technical challenges. For instant, in the past, prediction of in silico DNA methylation based on Gaussian distributions, that are relatively straightforward to implement, were used to describe the frequency distribution of CpGo/e ratios [START_REF] Gavery | DNA methylation patterns provide insight into epigenetic regulation in the Pacific oyster (Crassostrea gigas)[END_REF][START_REF] Park | Comparative analyses of DNA methylation and sequence evolution using Nasonia genomes[END_REF][START_REF] Dixon | Bimodal signatures of germline methylation are linked with gene expression plasticity in the coral Acropora millepora[END_REF][START_REF] Walsh | A functional DNA methylation system in the pea aphid, Acyrthosiphon pisum[END_REF] . But in many species, frequency distributions of CpGo/e ratios are complex or skewed and Gaussian distribution is not suitable. In our hands, only for 58 out of 83 cases (65%) Gaussian mixtures allowed for description of the distribution [START_REF] Bulla | Notos -a Galaxy tool to analyze CpN observed expected ratios for inferring DNA methylation types[END_REF] . These values are comparable to what was found by Bewick and colleagues who used CpG o/e ratios in transcriptomes of 124 species of which only 50 (40.32%) were described correctly with Gaussian mixtures [START_REF] Bewick | Evolution of DNA Methylation across Insects[END_REF] . We have also tested non-Gaussian distributions and the results were even less conclusive than Gaussian distributions: out of 83 only 41% delivered an exploitable result.

Therefore, we have developed a new tool, called Notos, to identify DNA methylation signatures within CpGo/e ratios based on kernel density estimations [START_REF] Bulla | Notos -a Galaxy tool to analyze CpN observed expected ratios for inferring DNA methylation types[END_REF] . This novel algorithm delivers robust descriptions of frequency distributions of CpGo/e ratios for up to 172,000 input sequences.

Here, we have applied this software to predict DNA methylation with CpGo/e ratios in a total of 634 species and to use the results in combination with publicly available experimental data to infer evolution of DNA methylation over the eukaryotic tree of life. We applied Notos on coding sequences coming from three databases (dbEST, CleanEST, and CDS/cDNA). Our results show clearly (i) that DNA methylation prediction by CpGo/e ratio is robust, (ii) that only four types of DNA methylation can be identified in all species despite their wide range of genome sizes, environments, body plans, reproduction types etc., and (iii) that DNA methylation types does not follow phylogeny but is consistent within clades suggesting evolutionary constraints. Taken together our analysis delivers arguments for the idea of the universality of the role of DNA methylation that is preserved through evolution.

Results

CDS and cDNAs are the less biased data and thus the best choice for a pan-species study

We focused in this study on gene body DNA methylation. Annotated genomes are now available for many species, but messenger RNA sequence data is even more abundant and mRNA is representative of gene body DNA sequences. They could therefore be used instead of DNA sequence, but mRNA data is for historical reasons stored in different forms and in different databases. We reasoned that data quality will be critical for providing unbiased estimation of DNA methylation in gene bodies and therefore conducted a comparative pilot study to identify the best possible data source for the subsequent pan-species study. We used coding/transcript sequences from full genome annotations (CDS), dbEST, and cleanEST (details in Supplementary file 2). A total of 127 species are in common between CDS and dbEST, and 92 species were in common between dbEST and cleanEST. Only 29 species were common to all three databases (Supplementary file 3). We produced Notos CpGo/e profiles for all intersecting datasets and proceeded to visual inspection. In 11 out of 29 cases (38%) we identified discrepancies in at least one out of the three profiles and decided to clarify their origins by a detailed analysis of the sequences under the differential peaks. An in-depth analysis revealed that these discrepancies were either due to contaminations during the sequencing process, reflect co-occurrence of other species, or are due to bias in data acquisition. For instance, for Trichoplax adhaerens, Anolis carolinensis (green anole lizard) and Cordyceps militaris one or two additional shoulder peaks in dbEST and CleanEST datasets. We isolated the sequences contained in these peaks (dbEST and Clean EST)

and performed a Blast2GO analysis with the aim to know their origins and functions. For the anole lizard (Supplementary figure 1), two peaks were isolated (peak 1: 0.92-1.08 and peak 2: 1. T. adherens is known to contain intracellular bacteria [START_REF] Driscoll | Bacterial DNA sifted from the Trichoplax adhaerens (Animalia: Placozoa) genome project reveals a putative rickettsial endosymbiont[END_REF] we believe that these sequences originate from them (Supplementary figure 2). For the mold C. militaris, two peaks were isolated (1.14-1.22, and 1.26-1.32). For the sequences under these peaks, homologies with other fungi sequences were found. We conclude that, in all three species, the additional modes occurred due to presence of sequences from other species, either through contamination during RNA extraction and library preparation, or as co-purification from naturally occurring symbionts or parasites.

In other species, we identified other sources of bias in the transcript data. For instance, in Bombyx mori an ovarian library cleanEST showed an additional weak shoulder peak. We isolated the 769 sequences under this peak (0.40-0.60). The gene ontology showed that all sequences coded for the ribosomal protein SA (RPSA). In human, RPSA genes are indeed highly expressed in the ovary but no data are available for other species. Nevertheless, it seems unlikely that the high abundance of RPSA ESTs reflects an expression bias. We speculate that the research interest of the submitters focused on this particular gene and that therefore many individual EST were submitted (Supplementary figure 3). Also in the duck (Anas platyrhynchos), we found a shoulder peak at 0.57 and 0.59 in data from dbEST and Clean EST. Sequences under this peak corresponded to an EST library exclusively composed of immunoglobulins (146 sequences), reflecting probably a bias introduced by specific research interests (Supplementary figure 4). Interestingly, these sequences had a CpGo/e ratio between 0.5 and 0.8 suggesting hypomethylation, and in human, immunoglobulin genes in lymphoid cells are indeed undermethylated during differentiation [START_REF] Storb | Methylation patterns of immunoglobulin genes in lymphoid cells: correlation of expression and differentiation with undermethylation[END_REF] .

Finally, when we compared profiles with two peaks (bimodality), where we had noticed differences between CDS (derived from genomes) and dbEST/CleanEST (mRNA) for three invertebrate and one plant species (Crassostrea gigas, Nasonia vitripennis, Nematostella vectensis and Oryza sativa) (Supplementary figures 5-8): mRNA derived profiles showed a higher peak in genes predicted to be methylated. Gene body methylation is suspected to increase transcription [START_REF] Bewick | Gene body DNA methylation in plants[END_REF][START_REF] He | Regulation and function of DNA methylation in plants and animals[END_REF] . The principal differences between CDS and EST data is that for the former only one FASTA sequence per gene is considered while for the later potentially several FASTA sequences for a gene could be present. We therefore hypothesized that RNA abundance induced the bias in EST data. To test this hypothesis, we performed a RNA-seq analysis in these four species. We found that genes under the low CpGo/e peak (presumably hypermethylated) show higher median RNA amounts than genes under the high

CpGo/e peak (this presumably hypomethylated). mRNA FPKM medians are 1.95 to 5.45 higher in presumably hypermethylated gene bodies (Supplementary figure 5678). We conclude that this expression difference is probably the origin of the bias in EST datasets.

In summary, dbEST and cleanEST have the advantage of being large repositories with data for many species, but for the purpose of our study we considered them too noisy. A complete list of species for each dataset is in supplementary file 4.

CpGo/e clustering identifies four types of gene body DNA methylation

After having firmly established that cDNA provides an unbiased data basis, CpG o/e clustering was carried out on the 142 species for which CDS or cDNA were available. Parameters for mode number (n), mode positions (Mo), skewness (sk), and standard deviation (SD) of CpGo/e values were iteratively changed using species with known gene body DNA methylation. For further analysis, we used the following features that produced four clusters of CpGo/e: (cluster 1) species with one mode Mo ≥ 0.69 and SD < 0.12, (cluster 2) species with one mode Mo ≥ 0.69 and SD ≥ 0.12, (cluster 3) species with one mode Mo < 0.69, and (cluster 4) species with (a) two modes or (b) one mode and a skewness smaller than -0.04. Results are in supplementary file 5 and supplementary figure 9. We then associated the four clusters with known methylation types.

Fourteen species (9.72 %), from different phylogenetic groups (e.g. Ascomycota, Apicomplexa, Basidiomycota, Plathyhelminthes and Arthropoda) constitute the cluster 1. All the species have a

CpGo/e mode position mode above 0.69 (the mean CpGo/e peak position is 1.00), a weak negative skewness (mean absolute Q50 skewness = -0.0019) and a narrow standard deviation (mean SD = 0.11). For 4 species (29%), experimental data on DNA methylation was available in the literature. All these species showed either absence of DNA methylation in the gene bodies or extremely low levels (Supplementary file 5). The latter was found in only one species (Chlamydomonas reinhardtii) where WGBS revealed methylation in exons but still it was 20-30 times weaker compared to other plant species in the same study. We qualify cluster 1 as "ultra-low gene body methylation" (type 1).

It could be argued that absence of gene body methylation is simply a consequence of absence of enzymatic methylation activity. We therefore performed a metanalysis of existing literature data concerning DNMT presence. In cluster ). We qualify this cluster as with "gene body methylation" (type 3).

Finally, cluster 4 contains species that show bimodality or are strongly negatively skewed CpGo/e distributions (mean absolute Q50 skewness = -0.0424, mean CpGo/e position of mode 1 is 0.54, of mode 2 0.85). Twenty-seven species (18.75 %) from different phylogenetic groups compose this cluster (apicomplexa, cnidaria, nematoda, arthropoda, mollusca, tunicata, embryophyta). Five species are strongly negatively skewed and 15 species are bimodal. We found information on DNA methylation for 10 species (50%). All species show a mosaic type of methylation with DNA regions of ultra-low methylation interspersed with regions of strong methylation. Eleven species out of 15 that were studied have the three DNMT (1, 2 and 3), two had DNMT1 and 2, and two DNMT1 and 3 with uncertainty about DNMT2 (Supplementary file 5). Species in this cluster were considered as "mosaic type DNA methylation" (type 4).

Decision criteria are summarized in Figure 1.

DNA methylation types do not follow the tree of life but are consistent within major clades Among unicellular organisms, kinetoplastids are unmethylated, while alveolate protists are generally methylated (Figure 2) with secondary loss of methylation that can lead to weaker methylation level or mosaic methylation. In flowering plant species, we differentiate either high

(

Discussion

Evolution is based on the selection of phenotypic variants that must (i) confer a reproductive advantage to the individual, and (ii) are heritable, i.e. information how to generate the phenotypic variants in response to an environment are passed from parents to offspring. Heritability has traditionally been thought to be exclusively genetic, i.e. based on variations in the DNA sequence.

In this view, genetic information is then expressed under influence of environmental cues to bring about the phenotype, a process known as G x E→ P [START_REF] Bowman | Genotype × environment interactions[END_REF] . During the last 30 years it became however clear, that a substantial amount of heritable phenotypic variance can be coded by non-genetic means [START_REF] Danchin | Beyond DNA: integrating inclusive inheritance into an extended theory of evolution[END_REF] . We had earlier conceptualized this view as a systems approach to inheritance, that includes genetic, epigenetic, cytoplasmic and microbial elements that are interrelated by forward and reverse interaction 2 . These elements interact mutually, and with the environment, to give raise to the phenotype. In this concept, genetic information (the genotype) is only one of many elements that as part of an inheritance system providing heritable information, that the environment will shape into a phenotype. We define here 'inheritance system' as a system that is able to write and store, transmit, and receive hereditary information [START_REF] Lamm | The Standford Encyclopedia of Philosophy[END_REF] . The concept implies also that genotype and epigenotype cannot exist independent of each other, and are interrelated by forward action and feedback. This is different from the idea that sees the genome as hard-wired information that is controlled by the epigenome [START_REF] Laland | Does evolutionary theory need a rethink?[END_REF] . In the latter, the epigenome is conceptually closer to the (molecular) phenotype (i.e. product of the genotype) than to an element of the inheritance system itself.

The introduction of the epigenotype notion did not really solve the question, since theoretically each phenotype could just be the visible expression of its underlying epigenotype. Given the multiple facets phenotypes can acquire in living organisms, it is remarkable that, with very few exceptions, the genetic material and the genetic code remains extremely constant and thus universal [START_REF] Koonin | Origin and evolution of the genetic code: The universal enigma[END_REF] . In other words, there exists a single 'type' of genome. The origin of this universality of the genetic code remains enigmatic and controversial but whatever the origin is, it allows to transmit coded information from one generation to the next. These generations can understand the code since it uses a universal and constant key.

Given the presumably close relation between genotype and epigenotype we and others reasoned that the epigenotype and the epigenetic code should equally possess universality. The high conservation of histones and histone marks, and the conservation of methylation of cytosines suggests indeed this. Nevertheless, one could argue that the epigenetic code is simply entirely genetically determined. If this were true, we would expect that different DNA methylation types would correspond to the clades in taxonomical tree that are based on DNA sequence similarity. Our results do not support this view. Alternatively, DNA methylation types could entirely be determined by environmental conditions. In this case, similar environments should impose similar DNA methylation types. Neither our results, nor recent analyses of DNA methylation in invertebrates provide evidence of this. E.g. a very comprehensive study of DNA methylation in insects [START_REF] Bewick | Evolution of DNA Methylation across Insects[END_REF] did not find relations of methylation types to social behavior and the authors concluded that DNA methylation must have "more ubiquitous function". However, compared to the tremendous amount of genomic data that is available, epigenomic data is relatively sparse and biased, which is an obstacle to answer the question conclusively. In the present study, we coped with this caveat by using a hybrid approach in which we combined available experimental data on DNA methylation with results coming from a newly developed software that predicts gene body DNA methylation types with CpG o/e ratios. Our algorithm (based on the number of species positive predicted and true positives based on the literature) allowed for including species for which no experimental DNA methylation data existed. The PPV of the algorithm is excellent for mosaic methylation (PPV=1), and methylated gene bodies (PPV=0.875), but decreases then to 0.75 (low methylated) and 0.5 for ultra-low gene body methylation. This is due to the fact that our algorithm does not differentiate well between low and ultra-low methylation. If we consider the dataset as a whole, out of the 54 species with known DNA methylation types, 41 were predicted correctly (total PPV = 0.76).

There are some particularly interesting cases of "wrong" prediction. Cryptosporium parvum is a monoxenic unicellular parasite of vertebrates. It belongs to the apicomplexan its exact phylogenetic position is controversial. Exysted oocysts are the only stage that can be used to produce host DNA free genomic DNA preparations. Notos predicts clearly high gene body methylation but LC-ESI-MS did not detect 5mC in exysted oocysts purified from infected cattle [START_REF] Gissot | Toxoplasma gondii and Cryptosporidium parvum lack detectable DNA cytosine methylation[END_REF] . Genome analysis of C.

hominis to which C. parvum has only 3-5% sequence divergence [START_REF] Xu | The genome of Cryptosporidium hominis[END_REF] , showed that the number of genes is reduced (3,952 genes) compared to other apicomplexan, relying heavily on host gene activity. The genome shows also traces of integration of genes by lateral transfer. We hypothesize that either the progenitor DNA was methylated, or that cryptosporidium methylates DNA in the intracellular stages using the vertebrate host DNMTs.

Another peculiar case is the ciliate Tetrahymena thermophile. Also for this species Notos predicted high methylation while radioisotope labeling showed that Tetrahymena contains only N6-methyladenine but not 5mC [START_REF] Hattman | Comparative study of DNA methylation in three unicellular eucaryotes[END_REF] . T. thermophila and other ciliates use DNA elimination to remove approximately one-third of the genome, when the somatic macronucleus differentiates from the germline micronucleus. Histone 3 lysine 9 trimethylation (H3K9me3) is deposited on DNA destined for this elimination (reviewed in Bracht [START_REF] Bracht | Beyond transcriptional silencing: Is methylcytosine a widely conserved eukaryotic DNA elimination mechanism?[END_REF] ). Interestingly, in other ciliates, DNA methylation is used for the tagging of DNA to be eliminated. It might therefore be that Tetrahymena had used DNA methylation in the past and has lost this capacity relatively recently, so that we still see traces in the CpG o/e ratio.

In summary, Notos predicts very reliable mosaic and high gene body methylation without being entirely error free. We had earlier [START_REF] Bulla | Notos -a Galaxy tool to analyze CpN observed expected ratios for inferring DNA methylation types[END_REF] [START_REF] Nanty | Comparative methylomics reveals gene-body H3K36me3 in Drosophila predicts DNA methylation and CpG landscapes in other invertebrates[END_REF] and 0.80 [START_REF] Suzuki | CpG methylation is targeted to transcription units in an invertebrate genome[END_REF] , and Nematostella vectensis, 0.70 [START_REF] Nanty | Comparative methylomics reveals gene-body H3K36me3 in Drosophila predicts DNA methylation and CpG landscapes in other invertebrates[END_REF] , or Apis mellifera, 1.0 [START_REF] Nanty | Comparative methylomics reveals gene-body H3K36me3 in Drosophila predicts DNA methylation and CpG landscapes in other invertebrates[END_REF] ). It is conceivable that Mo could be slightly different for each major phylogenetical clade, and using more sophisticated clustering algorithms such as support vector machine clustering that can use multiple thresholds could still improve the PPV of our method. In addition, more experimental data on a wide range of organisms is urgently needed.

We find that there are four types of gene body methylation. Despite a much wider data basis in terms of phylogenetic clades, our results confirm earlier findings that concluded on three to four DNA methylation types [START_REF] Suzuki | DNA methylation landscapes: provocative insights from epigenomics[END_REF][START_REF] Yi | Computational approaches for understanding the evolution of DNA methylation in animals[END_REF] . This could be the result of a "frozen accident" situation in which methylation (e.g. type 1 and type 3) occurred randomly in early ancestors (since 5mC is coding neutral that would not have had an impact on translation), but with the establishment of a chromatin structure 5mC was recruited as epigenetic information carrier, and any change in DNA methylation type would have had a strong impact on genome function and thus fitness and was therefore maintained. Nevertheless, switching of methylation type has occurred in evolutionary time scales.

Our findings indicate that there were at least three large events of secondary loss of DNA methylation: in archaeplastida (the "true" plants) where we find one branch with high methylation and another with mosaic methylation (in monocotyledons), ultra-low or mosaic methylation in the apicomplexa branch of "protists", and one transition to ultra-low gene body methylation in Diptera (Figure 2). For D. melanogaster in the dipteran branch it was shown experimentally that only the 'writing' capacity of the epigenetic inheritance element was lost, not the receiving ('reading') capacity [START_REF] Lyko | Mammalian (cytosine-5) methyltransferases cause genomic DNA methylation and lethality in Drosophila[END_REF] . The reason for evolutionary switching between methylation types is not clear and arguments are controversial.

It has been proposed that secondary loss of DNA methylation occurs because its mutational costs outweighed its adaptive value [START_REF] Zemach | Genome-Wide Evolutionary Analysis of Eukaryotic DNA Methylation[END_REF] . Indeed, in mosaic type methylation it is the evolutionary stable "old" genes that are in the methylated compartments meaning that there must be stabilizing mechanisms that prevent mutations there. Therefore, it might not be the mutational costs but the costs of maintaining such mechanisms that becomes an evolutionary burden. It was an early observation that CpG containing codons are used much less in coding sequences of vertebrates, and mutations due to CpG methylation was considered a major cause for such codon bias [START_REF] Dixon | Evolutionary Consequences of DNA Methylation in a Basal Metazoan[END_REF] and therein.

Codon bias was observed also recently in the reef-building coral Acropora millepora [START_REF] Dixon | Evolutionary Consequences of DNA Methylation in a Basal Metazoan[END_REF] , and linked to mosaic methylation in this species. Again, phylogenetically old genes which are constitutively expressed are methylated and CpG depleted. The authors conclude that CpG methylation leads to mutations that establish a set of preferred codons in constitutively expressed genes. Once such codon bias is fixed, then alleles that control the abundance of appropriate tRNAs could have stronger effects more amenable to natural selection. The authors hypothesize that an advantage of mutation-driven codon bias that it would be beneficial for organisms with small population size or otherwise inefficient selection. Still another explanation for mosaic methylation was advanced by Gavery and Roberts 32 who speculated that hypo-methylated regions (here in the pacific oyster Crassostrea gigas) might have greater epigenetic flexibility and higher regulatory control than hyper-methylated ones. Mosaic methylation could also be the result of whole genome duplication (WGD) events as suggested for Oryza sativa [START_REF] Zemach | Genome-Wide Evolutionary Analysis of Eukaryotic DNA Methylation[END_REF] . In addition, we have shown that environmental conditions can influence on germ-line methylation in C. gigas that possess mosaic methylation, and that blocks of CpG methylation are added or removed preferentially in or around genes [START_REF] Rondon | Effects of a parental exposure to diuron on Pacific oyster spat methylome[END_REF] . One should keep in mind that DNA methylation is only one of many bearers of epigenetic information.

Another one, and probably the most difficult to capture is the topology of the interphase nucleus.

Using Hi-C data, Lieberman-Aiden et al. [START_REF] Lieberman-Aiden | Comprehensive Mapping of Long-Range Interactions Reveals Folding Principles of the Human Genome[END_REF] established that the human genome is divided into two compartments (A-B) with pairs of loci in compartment B showing higher interaction frequency at a given genomic distance than pairs of loci in compartment A. They concluded that compartment B is more densely packed (heterochromatic) than compartment A. Higher average DNA methylation was later found to be a good predictor for the open compartment A in human cell lines [START_REF] Fortin | Reconstructing A/B compartments as revealed by Hi-C using longrange correlations in epigenetic data[END_REF] but that link could be broken in cancer cells. This cannot be interpreted as DNA methylation being decisive for topologically associated domains (TAD) establishment since DNA methylation free organisms such as D. melanogaster also presents canonical A-B domains [START_REF] Sexton | Three-dimensional folding and functional organization principles of the Drosophila genome[END_REF] . But in drosophila, such TAD organization is not driven by long-lived interactions but rather relies on the formation of transient, low-frequency contacts [START_REF] Cattoni | Single-cell absolute contact probability detection reveals that chromosomes are organized by multiple, low-frequency yet specific interactions[END_REF] . We hypothesize therefore that DNA methylation actually impacts on the relative dynamics of formation of contacts in A and B compartments, possibly stabilizing them. It is tempting to speculate that one consequence of compartmentation of genomes dynamics by methylation is that this might create additional units of selection. Results from tunicates support this idea: Ciona CpGo/e ratios have different profiles (bimodal for C. intestinalis and unimodal for C. savignyi). The C. intestinalis methylome is predicted to be mosaic that corresponds to experimental observations [START_REF] Suzuki | Identical sets of methylated and nonmethylated genes in Ciona intestinalis sperm and muscle cells[END_REF] . Our prediction for C. savignyi is low methylation (cluster 2). Both species diverged from each other 184 (±15) Mya [START_REF] D'onofrio | How fast is the sessile Ciona?[END_REF] and their genomes are very different [START_REF] Small | Extreme genomic variation in a natural population[END_REF][START_REF] Kourakis | An organismal perspective on C. intestinalis development, origins and diversification[END_REF] . For instance, analysis of 18S rRNA sequences shows that the pairwise divergence of the two ciona species is slightly greater than that between human and e.g. birds [START_REF] Johnson | Noncoding regulatory sequences of Ciona exhibit strong correspondence between evolutionary constraint and functional importance[END_REF] . This is puzzling since developmental features, body plan, effective population size and environment are very similar, and even hybrids can be generated to the tadpole stage [START_REF] Byrd | Mechanism of the block to hybridization and selfing between the sympatric ascidiansCiona intestinalis andCiona savignyi[END_REF] . However, C. savignyi shows a genome wide average Single Nucleotide Polymorphism (SNP) heterozygosity of 4.5% while C. intestinalis, that has mosaic methylation, is genetically less polymorphic (1.5%) (reviewed in Veeman et al. [START_REF] Veeman | Vertebrate Embryogenesis: Embryological, Cellular, and Genetic Methods[END_REF] ). It is conceivable that the methylated C. intestinalis genome can generate sufficiently stable TADs so that genome x epigenome interactions can serve as heritable unit of selection, while in C. savignyi TADs are more dynamic because the relative weight of DNA methylation in the generation of stable heritable phenotypic variants is less important. Our prediction concords with very recent results showing that stress-induced DNA methylation changes in C. savignyi can occur but are highly ephemeral (<48-120 h), and thus not maintained through germline [START_REF] Huang | Rapid response to changing environments during biological invasions: DNA methylation perspectives[END_REF] .

In conclusion, our findings indicate that initially there were three types of gene body DNA methylation: 'primary no methylation', 'primary whole genome methylation', and 'primary mosaic methylation' that produced by secondary loss 'weak methylation', or 'secondary no methylation'.

These findings are in concordance with the idea that DNA methylation in gene bodies (i) uses three types of universal codes (low, high and mosaic)), and (ii) that it is an element of the inheritance system and not a molecular phenotype that results from genotype x environment interaction. This has immediate practical consequences: e.g. since there are three types of methylation codes, panspecies conclusions about the potential function of DNA methylation can only be drawn within the type (e.g. functional tests in vertebrates with high gene body methylation cannot be used to conclude on methylation function in mosaic type mollusks). In addition, if DNA methylation is part of the inheritance system then heritable phenotypic diversity can be produced by DNA methylation changes without changes in the DNA sequence. The notion that everything that is heritable is necessarily genetic should be abandoned.

Methods

Origin of sequences, data cleaning and Notos parameters

In this study, coding sequences (CDS) or cDNA sequences of 147 species were downloaded from Ensembl and VEGA databases. Expressed sequence tags (ESTs) were downloaded from two different databases: dbEST 70 (605 species) and CleanEST 71 (110 species) (Supplementary file 3).

We used Notos [START_REF] Bulla | Notos -a Galaxy tool to analyze CpN observed expected ratios for inferring DNA methylation types[END_REF] for the calculation and modelling of CpGo/e distribution with the three datasets, with a minimal length L=200 bp and formula 1 72 :

[Formula 1] ⁄ = × × -1
All the values outside the interval, and all the values with a score of 0 were removed. For each species, the number of mode, the position of mode(s), the Q50 skewness coefficient and the standard deviation (SD) were calculated.

Blast searches and gene ontology analysis

Database searches were done by Blastx searches against a local instance (ncbi-blast-2.2.30+) of non-redundant 'nr' with 20 maximum hits, an E-value of 0.001, and other parameters as default values. Gene ontology searches were performed with blast2go [START_REF] Conesa | Blast2GO: A universal tool for annotation, visualization and analysis in functional genomics research[END_REF] .

RNA seq analysis

RNAseq datasets for Nematostella vectensis, Nasonia vitripennis, Crassostrea gigas and Oryza sativa japonica were downloaded as fastq files from the European Nucleotide Archive and NCBI (details in supplementary file 3). For each dataset, the reads were filtered with a Fred quality score ≥26. Filtered reads were mapped on their reference genomes (downloaded from Ensembl, details in supplementary file 3) with RNA STAR [START_REF] Dobin | STAR: Ultrafast universal RNA-seq aligner[END_REF] on a local Galaxy instance (v2.4.0d-2). Resulting BAM files and the gff files (downloaded from Ensembl, details in supplementary file 3) with the coding sequences were used for FPKM estimations with Cufflinks [START_REF] Trapnell | Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation[END_REF] . Annotation gff-files were used to extract CDS in fasta format from their genomes and we calculated the CpGo/e ratios with Notos [START_REF] Bulla | Notos -a Galaxy tool to analyze CpN observed expected ratios for inferring DNA methylation types[END_REF] and detected modes (peaks). To compare FPKM for the genes under the peaks, a bandwidth of 0.2 (±0.1 around mode maximum) was arbitrarily chosen for the CpG o/e ratio. FPKMs were extracted and used for statistical analysis of expression level in gene bodies with low and high predicted methylation. Mood's median test was used with R [START_REF]R: A Language and Environment for Statistical Computing[END_REF] .

Meta-analysis of DNA methylation using literature data

For each species for which data was available in the above-mentioned databases, we searched the literature on Google scholar (as of April-June 2016) with the following keywords: DNA methylation, 5-methyl-cytosine, gene body, mosaic methylation, global methylation, DNA methylation pattern. Articles were obtained from Bib CNRS (https://bib.cnrs.fr/) and manually curated to obtain gene body methylation, and presence of DNMTs.

Clustering

To identify distinct subgroups within the 147 analyzed species, we generated descriptive analyses, considering both the KDE of the CpGo/e ratios and aggregating statistics based on it. The statistics we used were (1) the number of modes of the KDE, (2) the position of the modes, (3) the standard deviation SD of the CpGo/e ratios, (4) absolute Q50 mode skewness of the CpGo/e ratios, i.e. ,

[Formula 2]

+ 2 -

with Q 1 and Q 3 the 25 % and 75 % quantile of the CpGo/e ratios, respectively, and Mo the global mode of the KDE. We investigated several formulas for the skewness, deeming the absolute Q50 mode skewness the most informative for our analysis [START_REF] Bulla | Notos -a Galaxy tool to analyze CpN observed expected ratios for inferring DNA methylation types[END_REF] . For the sake of readability, we refer to the absolute Q50 mode skewness as "skewness" in what follows.

The four clusters into which we classify the species are specified in the result section. The values of the three thresholds used in the definition of the clusters were determined by evaluating the prediction performance of our approach depending on these three values, using 54 species for which the true methylation type had been determined experimentally (given in Supplementary file 5).

Hereby, the clusters correspond to the patterns used in that file like follows: Cluster 1 -not methylated; Cluster 2 -low methylated; Cluster 3 -(high) methylated / global methylation; Cluster 4 -mosaic methylation.

To determine the optimal threshold values, we employed a two-step approach. First, we searched the whole parameter space (i.e. the space of all possible values the thresholds can assume) using a Metropolis-Hastings algorithm to ensure that strongly deviating from the threshold values we chose manually always leads to a poor prediction performance. Second, we systematically searched the parameter space around the manually chosen values. That is, we evaluated the prediction performance on a grid of size 21³ = 9261, covering the following threshold values: skewness from -0.08 to -0.04 in steps of 0.002; peak position from 0.69 to 0.79 in steps of 0.005; SD from 0.11 to 0.21 in steps of 0.005. For the present work we used: skewness -0.04, peak position 0.69, and SD 0.12.

Due to the scarcity of the data, the optimal prediction (76 %, 41 out of 54 true) is achieved for a rather large set of threshold values. To judge the performance of our algorithm, it should be noted that for 7 out the 13 misclassified species the true and the predicted classifications are "not methylated" and "low methylated", or vice versa respectively. That is, the mistake made by the algorithm is rather small in these cases. The remaining four wrong predictions are actually peculiar cases that were discussed above.

The clustering was implemented using R version 3.4.0 (supplementary files 6 and 7). For visualizing the clustering, the R package dendextend has been used. Parameters were Rscript cluster.r -0.04 0.69 0.12 input_file_notos input_file_notos_bootstrap with first parameter

[skewness], second [Mo], and third [SD], and Notos outputfiles as input. Further details on our method can be found in [START_REF] Bulla | Notos -a galaxy tool to analyze CpN observed expected ratios for inferring DNA methylation types[END_REF] .

Tree of life

We recovered the taxonomic IDs of all investigated species from the NCBI taxonomy database (https://www.ncbi.nlm.nih.gov/taxonomy) and created a tabular file (.txt). We used this file to generate a Phylip tree file based on the classification in the NCBI taxonomy database with the NCBI common taxonomy tree online tool (https://www.ncbi.nlm.nih.gov/Taxonomy/CommonTree/wwwcmt.cgi) and designed a tree with the interactive Tree of life (version 4.0.2) (https://itol.embl.de) [START_REF] Letunic | Interactive tree of life (iTOL) v3: an online tool for the display and annotation of phylogenetic and other trees[END_REF] .

Figure legends

Notos parameters

Species with a single mode, weak asymmetry and a low standard deviation.

Species with a single mode and left asymmetry, or bimodality Species with a single mode, weak asymmetry and a large standard deviation.

Species with a single mode and CpGo/e peak position below 0.69. Ch lam yd om on as rei nh ard tii [1] Glycine max [3] Bo m by x m or i [4 ] Cavia porcellu s [3] Canis lupus familiaris [3] Ornithorhynchus anatinus [3] G a s te ro s te u s a c u le a tu s [2 ] [2 ] A s p e r g il lu s fl a v u s G a ll u s g a ll u s [3 ] A p is m e lli fe ra [4 ] Prunus persica [3] [4] Crassostrea gigas

Cluster profile description

Figure 1 :Figure 2 :

 12 Figure 1: Summary of decision grid for clustering of CpG o/e ratio distributions on species level. Figure 2: Schematic representation of the "Tree of Life" for 147 species, associated with the four different types of DNA methylation that were identified in this work. Numbers in brackets indicate DNA methylation types ("clusters") for each species. Line colors correspond to methylation types.

[ 4 ][ 2 ]

 42 with PM ≥ 0.1 Mode with 0.05 ≤ PM < 0.1 Mode with 0.01 ≤ PM < 0.05 Strigamia_maritima.Smar1.29.cds.all.fa_final Gaussian kernel with band width 1with PM ≥ 0.1 Mode with 0.05 ≤ PM < 0.1 Mode with 0.01 ≤ PM < 0.05 with PM ≥ 0.1 Mode with 0.05 ≤ PM < 0.1 Mode with 0.01 ≤ PM < 0.05 Nematostella vectensis [2 ] S c h is to s o m a ja p o n ic u m Blat tella germ anic a [3] Rhodnius prolixus [3] Nem atost ella vecte nsis G ad us m or hu a [3 ] [2 ] A n c y lo s to m a c e y la n ic u m Ec hi no co cc us gr an ul os us [3 ] S tro ng yl oc en tro tu s pu rp ur at us T r y p a n o s o m a c r u z i [2 ] [3] Solan um tuber osum T a e n io p y g ia g u tt a ta [3 ] A e d e s a e g y p ti [2 ] [ 2 ] W u c h e r e r ia b a n c r o f t i Zo ote rm op sis ne va de ns is [4] M e lo id o g y n e h a p la [2 ] N ip p o s tr o n g y lu s b r a s il ie n s is G o ri lla go ri lla [3 ] Ba be sia eq ui [3 ] [2 ] C o rd y c e p s m il it a ri s A n o p h e le s d a r li n g i [ 1 ]

[ 2 ]

 2 S c h iz o s a c c h a r o T ri tic u m a e st iv u m T he ile ria pa rv a[3 ]

2 ]

 2 Ec hin oc oc cu s mu ltil oc ula ris [ 2 ] A s p e r g il lu s t e r r e u s P la s m o d iu m fa lc ip a ru m [3 s o p h il a a n a n a s s a e [2 ]

[ 2 ][ 3 ][[ 2 ]X[ 3 ][ 3 ][ 2 ][ 3 ][[ 2 ]

 23233232 P u c c in ia g ra m in is O ry zi a s la tip e s[3 ] T o x o p la s m a g o n d ii[1 ] Ara bid op sis tha lian a C lo no rc hi s si ne ns is[2 ] F u s a ri u m o x y s p o ru m C u le x q u in q u e fa s c ia tu s[2 ] [ 2 ] L o a lo a [2 ] S u il lu s lu te u s [ 1 ] A s p e r g il lu s o r y z a e Ov is arie s [3] [2 ] T ri c h u ri s m u ri s [1 ] C o p r in o p s is c in e r e a N a e g le r ia g r u b e r ix [3 ] [1 ] N e u r o s p o r a c ra s s a Tae nia sol ium T h e ile ri a o ri e n ta lis [2 ] T r y p a n o s o m a b r u c e i [2 ] E im e ri a m a x im a [4 ] M e ta rh iz iu m a n is o p lia e Y a r r o w ia li p o ly t ic a E im e ri a a c e rv u li n a [4 ] Br as si ca ra pa A n o li s c a r o li n e n au s pl ex ip pu s [4 ] C h o n d r u s c r is p u s [ 2 ri s ti o n c h u s p a c if ic u s M o r g a n e ll a Sol anu m lyco per sicu m Bo s tau rus [3] [2 ] T u b e r m e la n o s p o r u m H om o sa pi en s [3 ] As ty an ax m ex ic an us [3 ] S o le n o p s is in v ic ta [2 ] C ry p to s p o ri d iu m p a rv u m [3 ] [4 ] S te in e r n e m a c a r p o c a p s a e M ac ac a m ul at ta [3 ] Sel agi nel la mo elle ndo rffii [2] Sch midt ea med iterr ane a A n o p h e A m ph im ed on qu ee ns la nd ic a Th ei le ria an nu la ta [3 on iu s m el po m en e [4 ] [ 2 ] M e lo id o g y n e in c o g n it a D ro s o p h il a m e la n o g a s te r [2 Br as sic a ole ra ce a [2] Capitella teleta [2 ] P o s ti a p la c e n ta A n a s p la ty r h y n c h o s [3 ] [4 ] S e ta ri a it a lic a A tt a c e p h a lo te s [2 ] [4 ] Z ea m ay s Or yct ola gu s cun icu lus [3] [4 ] P a r a s tr o n g y lo id e s tr ic h o s u r i N e c a to r a m e r ic a n u s g o c a n d id a [2 ] G ia r d ia in t e s t in a li s [ 2 ] M e le a g ri s g a ll o p a v o [3 ] [1 ] T ri c h o d e rm a a tr o v ir id e N a s o n ia v it ri p e n n is [4 ] ri c h in e ll a s p ir a li s O re o c h ro m is n ilo ti c u s [3 ] [2] Tribolium castaneum E n t a m o e b a h is t o ly ti c a [ 3 ] [2 ] A n c y lo s to m S c h is to s o m a m a n s o n O ry z a s a ti v a [1 ] C la vi ce p s p u rp u re a [ 2 ] A s p e r g il lu s n ig e Mn em iop sis lei dy i [2 ] H a e m o n c h u s c o n to r tu s [1 ] M a g n a p o r th e o r y z a e [ 2 ] H e t e r o r h a b d it is b a c t e r io p h o r a C io n a s a v ig n y i [ 2 ] type 1: ultra-low gene body methylation type 2: low gene body methylation type 3: gene body methylation type 4: mosaic type DNA methylation h i g h t o l o w o r m o s a i c h ig h to m o s a ic l o w t o u l t r a -l o w

  1, 6 out of 14 species have DNMT2 or TRDMT1, and one specie has DNMT1 and DNMT2. Only the de-novo methylase DNMT3 is absent in this cluster. /e position is 0.92) and a mean absolute Q50 skewness of 0.0012. However, in contrast to cluster 1, a wide standard deviation (mean SD = 0.18) has been observed. Literature data were available for 18 species (30%). For 6 species (Schizosaccharomyces pombe, Aspergillus flavus,

	Brugia malayi, Meloidogyne incognita, Tribolium castaneum, Drosophila melanogaster) no
	methylation was reported. Methylation in 6 species (Schistosoma mansoni, Schistosoma japonicum,
	Fasciola hepatica, Petromyzon marinus, Caenorhabditis elegans and Saccharomyces cerevisiae) is
	controversial since different authors come to different conclusions. Nevertheless, methylation
	seems to be very low. Only three species (Trichinella spiralis, Solenopsis invicta, Physcomitrella
	patens) showed DNA methylation in gene bodies. DNMTs were studied in 37 species. In this
	cluster, 16 species have just DNMT2. One species has DNMT1 or TRDMT1 only, and 6 species
	have all the DNMT (DNMT1, 2 and 3). Interestingly, two species (Selaginella moellendorffii,
	Absence of methylation does therefore not indicate necessarily absence of DNMT genes
	(Supplementary file 5).

Cluster 2 is constituted by 60 species (41.67 %), also from different phylogenetic groups (apicomplexa, oobionta, rhodonbionta, ascomycota, basidiomycota, nematoda, platyhelminthes, annelida, arthropoda, ctenophora, chordata, embryophyta). As in the cluster 1, species present in the second cluster have a mode position > 0.69 with a mean mode position very close to the first cluster (mean CpGoPhyscomitrella patens) has just DNMT1 and 3, and two just DNMT1 and 2 (Tribolium castaneum, Gasterosteus aculeatus). Also, for Strigamia martima DNMT1 and DNMT3 were found, but DNMT2 was not searched for (Supplementary file 5). We consider cluster 2 as "low gene body methylation" (type 2).

Species with a mode position ≤ 0.69 form the cluster 3 (mean CpGo/e position is 0.45). The skewness is positive and larger than in the two first clusters (mean absolute Q50 skewness = 0.0483). The standard deviation (mean SD = 0.19) is wider than the cluster 1 and 2. Forty-three species

(29.86 %)

are present in this cluster, belonging to various phylogenetic groups (apicomplexa, sponges and nematodes, arthropoda, a large panel of chordata, and embryophytes). Many of them are important model organisms. For 26 (60%) literature data on DNA methylation was found. Gene bodies are methylated. In 19 species DNMTs were analyzed. Twelve species have all three DNMTs. One species (Naegleria gruberi) has DNMT1 and 2, but methylation of DNA was not yet studied. Four species have only DNMT2. Two species has DNMT1 and 3, and one only DNMT1 (Supplementary file 5

  probably global) methylation in dicotyledons, and mosaic methylation in poaceae and potentially all monocotyledons. Fungi process in general ultra-low to weak methylation in the gene bodies. Also, platyhelminthes are characterized by low methylation in the coding regions. Gene bodies of

deuterostomes are in general strongly methylated. There are, however, peculiar cases, e.g. two tunicate species (Ciona intestinalis and C. savignyi) that diverged from each other 184 (±15) Mya are in two different clusters (types 4 and 2, respectively)

[START_REF] D'onofrio | How fast is the sessile Ciona?[END_REF]

with mosaic and weak methylation, probably due to secondary loss of methylation. Within lophotrochozoa, annelids show low gene body methylation, and all studied mollusks are of the mosaic type. Nematodes have in general weak gene body methylation. A particular interesting and heterogenic clade concerning methylation types are arthropods. All tested diptera (and potentially all antliophora) belong to the low methylation clusters 1 and 2, certainly due to secondary loss of methylation after splitting from its insect sister clades. All other orders show weak to high methylation with occasionally mosaic type, probably through secondary gain of local methylation.

Acknowledgements

The work on this tool was initiated during a meeting that had received funding of the French-Norwegian travel program AURORA. This work has been supported by Campus France and the Norges forskningsrad (program AURORA, nr. 34040YK) to C. Grunau and J. Bulla, the grant Felleslegat til fordel for biologisk forskning ved Universitetet i Bergen to J. Bulla, the ANR grant ANR-10-BLAN-1720 (EpiGEvol) to C. Grunau, a PhD grant for disabled students by the French Ministry of Education and Research to B. Aliaga, and a DFG return grant to I. Bulla (BU 2685/4-1). The authors are grateful to Céline Cosseau and Cristian Chaparro for helpful discussion.

Author Contributions Statement

CG and DD designed the study. BA performed the experiments. IB developed the mathematical models, wrote the associated code, and did the clustering. GM and BA worked on the phylogeny. All authors wrote and reviewed the manuscript.

Additional information Competing interests

The author(s) declare no competing interests.


Ferenc Denes

Xavier Joulia

Peter Lang
email: lang@mail.bme.hu

Validation of a New Double-Column System for Heteroazeotropic Batch Distillation by Experiments and Dynamic Simulation

Keywords: batch heteroazeotropic distillation, closed system, experiment, dynamic simulation

The separation of a binary heteroazeotropic mixture (1-butanol water) by batch distillation is studied in a new Double-Column System (DCS) by pilot plant experiments and dynamic simulation performed with a professional flowsheet simulator. The DCS is operated in closed mode. At the end of the process products of high purity were obtained with high recoveries in the reboilers. However the limited facility of setting accurately the ratio of the reboiler heat duties caused difficulties concerning the optimal operation of the new system.

Introduction

If the components of a mixture form a heteroazeotrope, or by the addition of an entrainer a heteroazeotrope can be formed, the separation can be performed by heteroazeotropic distillation (HAD). So far the batch heteroazeotropic distillation was applied in the industry only in batch rectifiers (equipped with a decanter, BR) in open operation mode. Different operation modes of the batch HAD and their strategies are presented in the book of [START_REF] Luyben | Design and Control of Distillation Systems for Separating Azeotropes[END_REF]. The batch HAD is studied experimentally among others by Ooms et al. (2013, ethyl acetate isooctane by using methanol and acetonitrile as entrainer).

We developed a new batch double-column HAD configuration (double-column system, DCS, Figure 1) operated in closed mode (without continuous product withdrawal), [START_REF] Denes | New double column system for heteroazeotropic batch distillation[END_REF]. The results of feasibility studies and dynamic simulation calculations showed that comparing with the BR significant reduction of the energy demand can be achieved. Then the DCS was experimentally validated for the separation of a binary heteroazeotropic mixture in a simple laboratory glass equipment [START_REF] Denes | Experimental validation of a new double-column system for heteroazeotropic batch distillation[END_REF]. After that the DCS was extended to a generalised double-column system (GDCS, [START_REF] Denes | Generalised closed double-column system for batch heteroazeotropic distillation[END_REF] for the batch HAD separation of binary homoazeotropic mixtures by using entrainer. The GDCS has three further degrees of freedom compared with the DCS (locations of two continuous feedings, reflux ratio in Column ) and is more advantageous for the separation of binary homoazeotropic mixtures with a heterogeneous entrainer. In the next a pilot plant was installed, by which comparing with the laboratory column much more information can be acquired about the process (temperature profiles, pressure drop etc.), and whose operation is more flexible (variable reflux ratio) and appropriate (higher decanter holdup). The goal of this paper is -to validate the DCS by pilot plant experiments, -to evaluate the results with dynamic simulation.

Pilot plant experiment

The pilot plant equipment and its operation

The columns (D = 80 mm, H = 4x1 m) of the pilot plant equipment (Figure 2) are filled with structured packing: Sulzer CY (Column ), Kühni Rombopak (Column ). Both columns have an own condenser, reflux divider and aftercooler. Hence both condensates can be refluxed bypassing the decanter, too. The heating medium (oil) is circulated and heated by electricity (maximal performances: Q = 6 kW, Q = 4 kW). The temperature of the oil entering into the reboiler heating spiral is set.

The charge (of 18.6 °C) contains 1-butanol (A, 7.5 dm 3 ) and water (B, 9.5 dm 3 ). Into the reboilers 8.0 dm 3 saturated organic and 9.0 dm 3 aqueous phases are filled. The decanter is filled with heterogeneous liquid (Figure 2). The two columns are operated simultaneously. There are three operational steps:

1. Boiling-up: until the liquid in the reboiler begins to boil. 2. Heating-up: until the vapour arrives at the top of the column.

3. Purification: until the end of the process.

In the boiling-up and heating-up steps both columns are operated under total reflux Since Heater generates too high vapour flow rate, therefore 60 % of the condensate was directly refluxed in order to avoid the too short residence time in the decanter.

The purification steps of both columns were started at the same time. However, the durations of the boiling-up and heating-up steps were different because of the limited controllability of the reboiler heat duties.

The charge in Reboiler would have boiled up much earlier than in the other reboiler therefore its heating was strongly reduced until the heating-up of the other column approached to the end (for 76 min). Table 1 shows the lengths of the different operation steps. (The optimal lengths if they are different are in parentheses.) The experiment was finished when either reboiler temperature did not significantly increase more.

Experimental results

The reboiler temperatures (Figure 3) verify that products of very high purity were obtained at the end of the purification step with very high recoveries: 1-butanol: 99.9 mass%, 7.4 dm 3 ; water: 99.9 mass%, 9.5 dm 3 . At the end of heating-up both reboiler temperatures are near to the azeotropic one.

Posterior simulation

Simulation method

We performed rigorous simulation applying the following simplifying assumptions: theoretical trays, constant volumetric liquid holdup on the trays and in the decanter, negligible vapour holdup. The model equations to be solved are non-linear ODEs (material and heat balances) and non-linear algebraic equations (VLLE, hold-up and physical property models etc.). The phase equilibrium is described by the NRTL model (BIPs, K -1 : 215.427 and 1468.34; = 0.3634). For the solution of the above equations the dynamic simulator of ChemCAD 6.4 (CCDCOLUMN) is applied. The model of the pilot plant is shown in Figure 4.

Input data

The number of theoretical stages of the columns (of holdup 25 cm 3 /stage) without condenser, decanter and reboiler is estimated to 20. The levels of the liquid phases (of 60 °C) are kept constant in the decanter. The net heat duties of the reboilers ( are so chosen to give the same durations of the operation steps as in the experiment (Figure 5).The calculations are started with dry plates. The boiling-up step of Column was started by 76 minutes later in order to avoid the waiting period of the experiment (Figure 5, Table 2).

Results

The posterior simulation verified the production of high purity products: 1-butanol: 99.9 mass%, 7.72 dm 3 ; water: 99.9 mass%, 9.48 dm 3 . The product quantities and qualities are in good agreement with those of the experiment. Figure 6 shows the evolution of the reboiler temperatures.

Conclusions

The separation of the mixture 1-butanol water by batch heteroazeotropic distillation is studied in a new closed Double-Column System (DCS) by pilot plant experiments and posterior dynamic simulation with ChemCAD 6.4. At the end of the process products of high purity were obtained in the reboilers with high recoveries. However the experiment highlighted the importance of the accurate control of the reboiler heat duties in order to operate the DCS at optimal conditions and to avoid the loss of time and energy. This is even more important for the GDCS used for the separation of a binary homoazeotropic mixture with a heterogeneous entrainer where the composition of the condensate may get outside the heterogeneous region if the ratio of the heat duties is not appropriate.

Figure 1 .

 1 Figure 1. Theoretical scheme of the DCS for the separation of a binary heteroazeotropic mixture

Figure 2 .

 2 Figure 2. DCS pilot plant equipment for the separation of 1-butanol water

Figure 3 .Figure

 3 Figure 3. Evolution of the reboiler temperatures (experiment)

Figure 5 .

 5 Figure 5. Operational parameters in the different operation steps

Table 1 .

 1 Duration of the different operation steps[min]

	Step	Column Time at the end	Duration	Column Time at the end	Duration
	Boiling up	148	148 (53)	93	93
	Heating-up	175	27	175	82 (63)
	Purification	357	182	357	182