Title: a Natural Communication System on Genome Evolution

Title: A Natural Communication System on genome evolution

suggests a punctuated equilibria pattern

Authors:Qi Wua, YadiWangb, Yun Dinga, c, ShuaiMaa, d, ZongminWub, *, FuwenWeia, *

Affiliations:

aKey Laboratory of Animal Ecology and Conservation Biology, Institute of Zoology, Chinese Academy of Sciences, Beijing, China.

bShanghai Key Laboratory for Contemporary Applied Mathematics, School of Mathematical Sciences, Fudan University, Shanghai, China.

cSchool of Life Sciences, Nanchang University, Nanchang, China.

dCollege of Life Sciences, University of the Chinese Academy of Sciences, Beijing, China.

*Correspondence to: Professor Fuwen Wei () or Professor Zongmin Wu ().

Supplementary Information

1. Materials and Methods
2. Construction of geodesic equations for the information curve surface.
3. Numerical computation approaches for arc length of the geodesic andevolutionary rate
4. Supplementary Results
5. Supplementary Discussions

Figure S1. Patterns of the information evolutionary rate

Figure S2. Number of punctuated pattern in virtual evolutionary process in tested genome data.

Table S1 Classification of evolution pattern base on ratio of evolutionary rate

Table S2Web sources of genome data.

Table S3Taxonomic relations of the selected species.

1. Materials and Methods

1.1Construction of geodesic equations for the information curve surface

Since DNA sequence consists of four kinds of nucleotides, there are 16 kinds ofnucleotide pairs, which could be shown as a set in a form of 4×4 table like

,

and

Consider the reverse compliment rule, the non-palindrome pairs could be merged. For

example, AA in plus strand of DNA is the same as TT in minus strand of DNA. So

there are only ten pairs left. one got

For convenience sake, we usex1, x2, x3 and x4to denote frequencies ofA, C, Gand Trespectively. So the frequencies of four nucleotides xi and the join frequencies of the

ten nucleotide pairs xij could be defined as matrices

,

and

Then the mutual information functional could be defined as

(1)

Subject to

1)

2) (2)

3), i, j, k=1, 2, 3, 4; 和 xij, xki=(M2)ij

Although xij=1 andxij=0 is mathematically possible, it never occurs in biology and the domain of definition excludes the cases. Since only nine of the ten nucleotide pairs are independent variables, the domain of definition is a nine dimensional space. And the mutual information function is then a nine dimensional curved surface in a ten dimensional space, named the mutual information curved surface (MICS).

The process of information transmission from the source to the destination (or the

evolutionary process of the source to the destination) can be modeled as a process of a

point moving along the geodesic from point S representing the source to point D

representing the destination on the information curved surface. Because of the

assumption of high efficiency, the movement must find the shortest distance on the

information curved surface. Then the length between the source point and the

destination point on the geodesic can be used to quantify the distance from the source

to the destination. Therefore the problem is simplified to calculating the length

between two points on the geodesic on the nine dimensional curved surface in a ten

dimensional space.

1.2 Numerical computation approaches for arc length of the geodesic andevolutionary rate

We want to get the temporal channel which is the shortest curve between two different

points according to the assumption of high efficiency. In other words, the manifold

M(x, I) is constructed of the proportion of four kinds of nucleotides and information.

A(xA, IA) and B(xB, IB) represent two different species and the channel is a geodesic

between A and B.

Now we solve the geodesic with the theory that the intersection angle between the

normal vector to the curve and the normal vector to the manifold is zero at the point

on the geodesic. This gives an iterative algorithm as follows[1]

1. Select a proper initial curve C between A and B on the manifold M. Then we get a series of initial points from curve C, such that:

(S18)

E. g.

(S19)

1. Define the error ej between the normal vector to the curve Cnc,jand the normal vector of M nm,jat the point Pj. If the curve C is a geodesic, ej = 0 for all j = 0, 1, … , n. If not, there exists one point Pj such that ej is not equal to zero at least.
2. Compute the normal to the manifold M with PDEs and the normal to the curve Cwith interpolation property at the pointPjfor all j = 0, 1, … , n. Then compute the error ej.
3. Get a new series of points{ Pj*(xj* , Ij*)| xj*=xj+T(ej),Tis a projection}. Use this new series as the initial points in the next iteration which starts from step 3 to step 5.
4. Get a new curve C* with the series of points{ Pj*|j = 0, 1, … , n } by interpolation.
5. With an increase in the number of iterations, while the error ej approximates zerofor all j = 0, 1, … , n, the curve C* approximates the geodesic. If all ejare lessthan a given error bound, we can end the iteration and obtain the requiredgeodesic.

1. Result of test the model with genome data

We collected 94 genomes from public databases (Supplementary Tables S2 and S3) covering multicellular animals, green plants, some multicellular fungi and some unicellular eukaryotes. We calculated the evolutionary rate using these genome data and our model. A rigorous test requires genome data for all ancestral species, which is not available. Instead, we calculated evolutionary rate using every species as a virtual ancestral species and the other 93 species as virtual descendent species and obtained 94×(94-1)=8742 virtual evolutionary processes. An implicit assumption was that the distribution of genome data points of selected species in the MICS was the same as the distribution of all real ancestral species. Of all the processes, we found three kinds of patterns of concave curves. First pattern: the left part of the curve was flat, suggesting the evolutionary rate remained stable the majority of evolutionary time since the ancestor. As the curve approached the descendent point of the right part it increased sharply, meaning that the evolutionary rate increased significantly and reached the maximum when evolution approached the end. Second pattern: the primal trends were similar but with the difference that it reached the maximum just before the end of the evolution process. Third pattern: the primal trends were completely different and the evolutionary rates have two kinds of values in both positive and negative directions (Figure S1).

Firstly, we consider nucleotide pair with k=0, or the directly neighbored nucleotide pairs. We calculated the ratio of stable rate the the highest rate in the process of one direction evolution. Of the whole matrix, 7839 of 8742 grids have one direction (pattern one and two), which suggested that the majority of the geodesic in the MICS doesn’t meet any “ridges” or “peaks”. Of the 7839 one-direction processes, there are 1486 processes whose end rate are one order of magnitudes higher that the stable rate (ratio value equals 10) (Table S1). The result showed that a fairly large amount of one direction evolution processes (nearly one fifth) share a one-order-of-magnitudes-higher rate difference, which could be explained by evolution pattern of punctuated equilibrium. There are two virtual evolution processes of every two species. One has the direction of mutual information increasing while the other has an opposite one of mutual information decreasing. The result showed that all the 1486 processes are of the direction of mutual information decreasing.

Previous researches have manifested that the statistical correlation is short range dominant within one genome[12]. And evolution indicated an mutual-information increasing pattern through a coarse-grained average[12]. In this work we tested the influence of k value to the ratio of evolution pattern of punctuated equilibrium, combined with the two cases of mutual information increasing and decreasing respectively. The summary results are in Figure S2. It could be seen that as the k distance extends, the total number of the virtual evolution processes with one-order-of-magnitudes-higher rate tends to reductive. As the same time, more processes with one-order-of-magnitudes-higher rate showed the mutual information decreasing direction. Theses results are accordant with the previous researches and suggested such an evolutionary insight that positive selection tends to increasing the mutual information, negative selection tends to keep the mutual information stationary, while neutrality tends to dissipate the mutual information to either positive or negative direction.

1. Discussion

The re-measured time by information satisfies additivity. Consider three affairshappening sequentially A, B and C. Let TAB represent the time of the process herebyA happened and then B happened, TBC represents the time of the process whereby Bhappened and then C happened, and TABC represented the time of the process herebyA happened and then B happened and then C happened. Therefore:

(S20)

However, in the re-measured time definition, the time for the process of A happened

and then B happened and then C happened is different from the time for the process

that A happened and then C happened without B happening, or:

(S21)

This is because the traditional time is defined in one dimension while re-measured

time is defined in a multidimensional information curved surface. So time has

different routes to accommodate different intermediate affairs.

In the model of NCS, evolutionary novelties related with punctuated equilibrium were assumed to be determined by the statistical correlation of the whole DNA sequence of a genome instead of certain genes. While in present genetic theory it is gene, with the chemical basis of a trench of DNA, that determines the heritable feature of living creatures[2,3]. We admit that it is a classical approach for biologists using alignment-based method to identify gene by finding sequence similarity and then attributing genetic features to certain genes. However, it has also been a long-debated issue that the features of evolutionary novelty are too complicated to be explained by gradually accumulated variation of a few genes[4,5]. Researchers in development biology or system biology tended to introduce network consist of a number of genes to explain complicated features, which may also decrease the importance of one single gene[6]. There have also been some literature arguing a holistic perspective on genome in evolutionary theory[3,7]. Actually, nucleotide statistical correlations of genome DNA sequence reflect some over-all effects of variation distribution along chromosomes of a genome, the distribution of which may caused by genome replication, genome rearrangement and so force[8]. Evidences have been accumulated that genome rearrangement as well as copy number variation are responsible for the punctuated phynotypic variation in both evolution and tumorgenesis[9, 10]. And the approaches of nucleotide statistical correlations have been widely accepted as a genome signature of different species and have been applied widely in evolutionary biology[8,11-13].

As for the model itself, one of the most inconceivable idea could be that different affair sequences have different times, which suggests that different evolutionary processes have different times or that evolution time is lineage specific, even if the descendent species share a common ancestor. One can argue that the astronomic cycle for year and day is the same for all lives; however, animals with different body sizes have different time[14]. Although species cannot escape physical time measured by the day-and-night and annual seasons, these cycles have different meaning; for example, because of small and large body sizes, long and short generation times, etc. It is reasonable to deduce that time is different in the evolution of species passing a series of speciations and species known as living fossils.

Time measured by physical cycles is a dimensional quantity, but the channel time defined by probability is not because of the definition of information used to quantify time. We used mutual information and the perennity of time taken for duplicating the ratio as a unit. When using another definition, such as Laxton’s diversity[15], time may be a dimensional quantity. For geneticists, neither mutual information[8,13] nor Laxton’s diversity[15] are familiar methods. Sequence alignment[16]and related methods[17]are either difficult to use directly in information quantification or when missing genome-scale information (e.g. distribution of variation among chromosomes). It remains an essential challenge for biology to develop new statistical methods to explain and quantify variation in information at the genome scale.

References:

[1] RaviKumara G. V. V., Srinivasan, P., DevarajaHollaa V., Shastrya K. G. & PrakashcB. G. Geodesic curve computations on surfaces. Computer Aided Geometric Design. 20, 119-133 (2003).

[2] Krebs J. E., Goldstein E. S. & Kilpatrick S. T. in Gene X (Higher Education Press, 2010)

[3] Hartl D. L. and Clark A. G. in Principles of Population Genetics (Sinauer Associates, 2007)

[4] Heng H. H. Q. The genome-centric concept: resynthesis of evolutionary theory. Bioessays. 31, 512-525 (2009).

[5] Gould S. J. and Eldredge N. Punctuated equilibrium comes of age. Nature. 366, 223-227 (1993).

[6] Kitano H. Computational systems biology. Nature 420, 206-210 (2002).

[7] Tenreiro Machado J. A. T. Accessing complexity from genome information. Commun Nonlinear Sci Numer Simulat. 17, 2237-2243 (2012).

[8] Karlin S. and Ladunga. Comparisons of eukaryotic genomic sequences. Proc. Natl. Acad. Sci. USA. 91, 12832-12836 (1994).

[9] Lonnig W-E and Saedler H. Chromosome rearrangements and transposable elements. Annu. Rev. Genet. 36, 389–410 (2002).

[10] Hastings, P. J., Lupski, J. R., Rosenberg, S. M. & Ira, G. Mechanisms of change in gene copy number. Nat Rev Genet. 10, 551-564 (2009).

[11] Song K., Ren J., Reinert G., Deng M., Waterman M. & Sun F. New developments of alignment-free sequence comparison: measures, statistics and next-generation sequencing. Brief Bioinform. 15:343-353 (2014).

[12] Xu Z. & Hao B. CVTree update: a newly designed phylogenetic study platform using composition vectors and whole genomes. Nucleic Acids Res. 37(Web Server issue):W174-W178 (2009).

[13] Luo L. F., Lee W. J., Jia L. J., Ji F. M. & Tsai L. Statistical correlation of nucleotides in a DNA sequence. Phy Rev E.58, 861-871 (1998).

[14] Schmidt-Nielsen K. in Scaling: Why Is Animal Size So Important? (Cambridge University Press, 1984)

[15] Laxton R. R. The measure of diversity. J. Theor. Biol. 70, 51-67 (1978).

[16] Karlin S. and Altschul S. F. Methods for assessing the statistical significance of molecular sequence features by using general scoring schemes. Proc. Natl. Acad. Sci. USA. 87, 2264-2268 (2005).

[17] Wakeley J. in Coalescent Theory: An Introduction (Roberts & Company Publishers, 2009)

Figure S1. Patterns of the information evolutionary rate. The x axis is non-dimensional channel time; the y axis is the information evolutionary rate. The source point is the virtual ancestral genome in the ICS. The destination point is the virtual descendent genome in the ICS. A. An example of pattern with genome of rifleman (Acanthisitta chloris) as virtual ancestral species and genome of cow (Bos taurus) as descendent species, showing the highest value of rate at the point of the destination. The curve rises, indicating a faster rate of increase. At the point that the curve reaches the highest point it drops to zero, meaning that evolution has finished at the destination point. B. An example shows a pattern with the highest value of rate just before the point of the destination. Similarly to pattern A, the curve rises, indicating a faster rate of increase. After the point that the curve reaches the highest point it slows down a little and then drops to zero, suggesting that evolution has met some upslope on the ICS before it has finished at the destination point.The ancestral species for this example was rifleman and the descendent species was button mushroom (Agaricus bisporus (H97)). C. An example of a pattern showing both positive and negativevaluesfor information evolutionary rate, which means that evolution has crossed some peak or ridge on the ICS. The ancestral species for this example was rifleman and the descendent species was Aspergillus nidulans.

Figure S2. Number of punctuated pattern in virtual evolutionary process in tested genome data. The positive and negative direction of Y axis represents virtual evolutionary process with mutual information increasing and decreasing respectively. As mentioned in the main text, all the virtual evolution processes for given mutual information direction for given k value are divided into three groups, the group of bidirectional processes, the group that the ratio of end rate to stable rate is no higher than ten, the group that the ratio is higher than ten. The gray color denotes the group of bidirectional processes of both mutual information direction. The brown color denotes the group with the ratio no higher than ten in mutual information increasing direction. The yellow color denotes the group with the ratio higher than ten in mutual information increasing direction. The dark blue color denotes the group with the ratio no higher than ten in mutual information decreasing direction. The light blue color denotes the group with the ratio higher than ten in mutual information decreasing direction.

Table S1Web sources of genome data.

Species Name / source / website
Homo sapiens / ENSEMBL /
Pan troglodytes / ENSEMBL /
Gorilla gorilla / ENSEMBL /
Macaca mulatta / ENSEMBL /
Mus musculus / ENSEMBL /
Rattus norvegicus / ENSEMBL /
Cricetulus griseus / NCBI / ftp://ftp.ncbi.nlm.nih.gov/genomes/Cricetulus_griseus/
Bos grunniens mutus / NCBI / ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF_000298355.1_BosGru_v2.0/
Bos taurus / ENSEMBL /
Ovis aries / ENSEMBL /
Sus scrofa / ENSEMBL /
Lipotes vexillifer / NCBI / ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA_000442215.1_Lipotes_vexillifer_v1/
Balaenoptera acutorostrata / UCSC / ftp://hgdownload.soe.ucsc.edu/goldenPath/balAcu1/
Ailuropoda melanoleuca / ENSEMBL /
Ursus maritimus / BGI /
Aliurus fulgen / - / Unpublished sequencing data from Wei’s group
Mustela putorius furo / ENSEMBL /
Canis lupus familiaris / ENSEMBL /
Odobenus rosmarus / NCBI /
Panthera tigris / BGI /
Equus caballus / ENSEMBL /
Trichechus manatus latirostris / UCSC / ftp://hgdownload.soe.ucsc.edu/goldenPath/triMan1/
Dasypus novemcinctus / ENSEMBL /
Monodelphis domestica / ENSEMBL /
Sarcophilus harrisii / ENSEMBL /
Ornithorhynchus anatinus / ENSEMBL /
Anas platyrhynchos / BGI /
Aptenodytes forsteri / BGI /
Calypte anna / BGI /
Chaetura pelagica / BGI /
Charadrius vociferus / BGI /
Columba livia / BGI /
Corvus brachyrhynchos / BGI /
Cuculus canorus / BGI /
Egretta garzetta / BGI /
Falco peregrinus / BGI /
Gallus gallus / ENSEMBL /
Geospiza fortis / BGI /
Haliaeetus leucocephalus / BGI /
Manacus vitellinus / BGI /
Meleagris gallopavo / BGI /
Melopsittacus undulatus / UCSC / ftp://hgdownload.soe.ucsc.edu/goldenPath/melUnd1/
Nipponia nippon / BGI /
Ophisthocomus hoazin / BGI /
Picoides pubescens / BGI /
Pygoscelis adeliae / BGI /
Struthio camelus / BGI /
Taeniopygia guttata / ENSEMBL /
Tinamus guttatus / BGI /
Anolis carolinensis / ENSEMBL /
Chrysemys picta / UCSC / ftp://hgdownload.soe.ucsc.edu/goldenPath/chrPic1/
Xenopus tropicalis / ENSEMBL /
Latimeria chalumnae / ENSEMBL /
Danio rerio / ENSEMBL /
Oryzias latipes / ENSEMBL /
Takifugu rubripes / ENSEMBL /
Tetraodon nigroviridis / ENSEMBL /
Petromyzon marinus / ENSEMBL /
Ciona intestinalis / ENSEMBL /
Ciona savignyi / ENSEMBL /
Branchiostoma floridae / JGI /
Saccoglossus kowalevskii / NCBI / ftp://ftp.ncbi.nlm.nih.gov/genomes/Saccoglossus_kowalevskii/
Strongylocentrotus purpuratus / NCBI / ftp://ftp.ncbi.nlm.nih.gov/genomes/Strongylocentrotus_purpuratus/
Drosophila melanogaster / ENSEMBL /
Daphnia pulex / JGI /
Metaseiulus occidentalis / NCBI / ftp://ftp.ncbi.nlm.nih.gov/genomes/Metaseiulus_occidentalis/
Caenorhabditis elegans / ENSEMBL /
Ascaris suum / wormbase /
Capitella teleta / JGI /
Helobdella robusta / JGI /
Aplysia californica / UCSC / ftp://hgdownload.soe.ucsc.edu/goldenPath/aplCal1/
Lottia gigantea / JGI /
Schistosoma japonicum / SGST / ftp://lifecenter.sgst.cn:2121/subjectData/schistosoma/
Schistosoma mansoni / sanger / ftp://ftp.sanger.ac.uk/pub/project/pathogens/Schistosoma/mansoni/
Hydra magnipapillata / NCBI / ftp://ftp.ncbi.nlm.nih.gov/genomes/Hydra_magnipapillata/
Nematostella vectensis / JGI /
Amphimedon queenslandica / NCBI / ftp://ftp.ncbi.nlm.nih.gov/genomes/Amphimedon_queenslandica/
Trichoplax adhaerens Grell-BS-1999 / JGI /
Monosiga brevicollis / JGI /
Agaricus bisporus (H97) / JGI /
Coprinopsis cinerea okayama7#130 / JGI /
Pleurotus ostreatus / JGI /
Schizophyllum commune H4-8 / JGI /
Batrachochytrium dendrobatidis JAM81 / JGI /
Aspergillus niger / JGI /
Saccharomyces cerevisiae / ENSEMBL /
Phanerochaete carnosa / JGI /
Encephalitozoon cuniculi / NCBI / ftp://ftp.ncbi.nlm.nih.gov/genomes/Fungi/Encephalitozoon_cuniculi_uid155/
Arabidopsis thaliana / ENSEMBL /
Oryza sativa / ENSEMBL /
Phyllostachys heterocycla / NCGR /
Physcomitrella patens / ENSEMBL /
Ostreococcus lucimarinus / ENSEMBL /
Chlamydomonas reinhardtii / ENSEMBL /

Table S2Taxonomic relations of the selected species.