Estimation of the true evolutionary distance under the fragile breakage model
 Nikita Alexeev^{1}Email author and
 Max A. Alekseyev^{1}
Published: 24 May 2017
Abstract
Background
The ability to estimate the evolutionary distance between extant genomes plays a crucial role in many phylogenomic studies. Often such estimation is based on the parsimony assumption, implying that the distance between two genomes can be estimated as the rearrangement distance equal the minimal number of genome rearrangements required to transform one genome into the other. However, in reality the parsimony assumption may not always hold, emphasizing the need for estimation that does not rely on the rearrangement distance. The distance that accounts for the actual (rather than minimal) number of rearrangements between two genomes is often referred to as the true evolutionary distance. While there exists a method for the true evolutionary distance estimation, it however assumes that genomes can be broken by rearrangements equally likely at any position in the course of evolution. This assumption, known as the random breakage model, has recently been refuted in favor of the more rigorous fragile breakage model postulating that only certain “fragile” genomic regions are prone to rearrangements.
Results
We propose a new method for estimating the true evolutionary distance between two genomes under the fragile breakage model. We evaluate the proposed method on simulated genomes, which show its high accuracy. We further apply the proposed method for estimation of evolutionary distances within a set of five yeast genomes and a set of two fish genomes.
Conclusions
The true evolutionary distances between the five yeast genomes estimated with the proposed method reveals that some pairs of yeast genomes violate the parsimony assumption. The proposed method further demonstrates that the rearrangement distance between the two fish genomes underestimates their evolutionary distance by about 20%. These results demonstrate how drastically the two distances can differ and justify the use of true evolutionary distance in phylogenomic studies.
Keywords
Background
Genome rearrangements are evolutionary events that shuffle genomic architectures. Most frequent genome rearrangements are reversals (that flip segments of a chromosome), translocations (that exchange segments of two chromosomes), fusions (that merge two chromosomes into one), and fissions (that split a single chromosome into two). These four types of rearrangements can be modeled by DoubleCutandJoin (DCJ) operations [1], which break a genome at two positions and glue the resulting fragments in a new order.
The ability to estimate the evolutionary distance between extant genomes plays a crucial role in many phylogenomic studies. Often such estimation is based on the parsimony assumption, implying that the distance between two genomes can be estimated as the rearrangement distance equal the minimal number of genome rearrangements required to transform one genome into the other. However, in reality the parsimony assumption may not always hold, emphasizing the need for estimation that does not rely on the (minimal) rearrangement distance. The evolutionary distance that accounts for the actual (rather than minimal) number of rearrangements between two genomes is often referred to as the true evolutionary distance.
While there exists a method for estimation of the true evolutionary distance [2], it however implicitly assumes that genomes can be broken by rearrangements equally likely at any position in the course of evolution. This assumption, known as the random breakage model (RBM) of chromosome evolution [3, 4], has been refuted in favor of the more rigorous fragile breakage model (FBM) [5] postulating that only certain “fragile” genomic regions are prone to rearrangements. The FBM is supported by many recent studies of various genomes (e.g., see references in [6]). The RBM can be viewed as an extremal case of the FBM, where every genomic region is fragile.
In the current study, we propose a new method for estimating the evolutionary distance between two genomes with high accuracy under the FBM. We assume that the given genomes are represented as sequences of the same blocks (synteny blocks or orthologous genes) and evolved from a common ancestral genome with a number of DCJs. Our method estimates the total number of DCJs on the evolutionary path between such genomes. The results of using our method on a simulated dataset show a high level of precision. We also analyze yeast genome data and show that some, but not all pairs of yeast genomes fall under the parsimony assumption.
The subtle difference between the RBM and FBM from the perspective of true evolutionary distance estimation is the (in)ability to count the number of fragile genomic regions. While breakpoints (block adjacencies present in one genome and absent in the other) definitely represent fragile regions, the shared block adjacencies are treated differently under the two models. Namely, under the RBM a shared block adjacency still represents a fragile region, which just happened to remain conserved across the two genomes by chance. Under the FBM, a shared block adjacency may or may not be fragile.
Methods
Breakpoint graphs and DCJs
We start our analysis with circular genomes (i.e., genomes with circular chromosomes) and address linear genomes later. We represent a genome with n blocks as a genome graph composed of n directed block edges encoding blocks and their strands, and n undirected adjacency edges encoding adjacencies between blocks.
We say that a blackred cycle is an ℓcycle if it contains ℓ black edges (and ℓ red edges) and let c _{ ℓ }(P,Q) be the number of ℓcycles in G(P,Q). We refer to 1cycles as trivial ^{1} and to the other cycles as nontrivial. The vertices of nontrivial cycles are called breakpoints.
Evolutionary model
To estimate the true evolutionary distance between genomes P and Q on the same set of blocks, we view the evolution between them as a discrete Markov process that transforms genome P into genome Q with a sequence of DCJs occurring independently. The process starts at genome X=P and ends at X=Q, and corresponds to a transformation of the breakpoint graphs starting at G(P,P) (formed by a collection of trivial cycles) and ending at G(P,Q). The number of DCJs k in this transformation represents the true evolutionary distance between genomes P and Q.
We remark that under the FBM, the number of trivial cycles (if any) in G(P,Q) is an obscure parameter, since they may correspond to solid (nonfragile) regions as well as to fragile regions that just happen to remain conserved (in both P and Q) by chance.^{2} Furthermore, there may exist such conserved fragile regions within the blocks forming P and Q and thus such regions cannot be observed in G(P,Q). To account for all these possibilities, we assume that P and Q are composed of a large unknown number n of solid regions interspersed with the same number of fragile regions, some of which remain conserved by chance. In other words, we assume that the given blocks of genomes P and Q are formed by (invisible) solid regions and conserved fragile regions. Let P _{ n } and Q _{ n } denote these representations of P and Q as sequences of the solid regions. We view genome Q _{ n } as obtained from P _{ n } with a sequence of k DCJs each operating on two randomly selected fragile regions (i.e., adjacencies between solid regions).
The crucial observation is that while we do not know the number n of solid regions, the breakpoint graphs G(P,Q) and G(P _{ n },Q _{ n }) have the same cycle structure, except for trivial cycles. That is, we have c _{ ℓ }(P _{ n },Q _{ n })=c _{ ℓ }(P,Q) for every ℓ≥2, implying, in particular, that b(P _{ n },Q _{ n })=b(P,Q), c(P _{ n },Q _{ n })=c(P,Q), and d(P _{ n },Q _{ n })=d(P,Q). Indeed, if genomes P ^{′} and Q ^{′} are obtained from P and Q by replacing a single block a with two consecutive smaller blocks a _{1},a _{2}, then G(P ^{′},Q ^{′}) can be obtained from G(P,Q) by adding one trivial cycle (corresponding to the shared adjacency a _{1},a _{2}). Since the genomes P _{ n } and Q _{ n } can be obtained from P and Q with a number of such operations, the breakpoint graphs G(P,Q) and G(P _{ n },Q _{ n }) may differ only in the number of trivial cycles.^{3}

c _{ ℓ }=c _{ ℓ }(P _{ n },Q _{ n })=c _{ ℓ }(P,Q) for any ℓ≥2, i.e., the number of ℓcycles in G(P,Q);

\(b=b(P_{n},Q_{n}) = b(P,Q) = \sum _{\ell \geq 2} \ell \cdot c_{\ell }\), the number of broken fragile regions between P and Q, which is also the number of synteny blocks between P and Q, or half of the total length of all nontrivial cycles in G(P,Q);

\(d = d(P_{n},Q_{n}) = d(P,Q) = b  \sum _{\ell \geq 2} c_{\ell }\), the DCJ distance between P and Q;

c _{1}=c _{1}(P _{ n },Q _{ n }), the number of trivial cycles in G(P _{ n },Q _{ n });

n=n(P)=n(Q), the number of fragile regions in each of genomes P and Q, half the total length of all cycles in G(P _{ n },Q _{ n });

k=k(P,Q), the number of DCJs in the Markov process, the true evolutionary distance between P and Q.
Extension to linear genomes
To analyze linear genomes, we add one artificial adjacency between telomeres for each chromosome and consider the resulting circular genomes. Since the number of chromosomes is often negligible as compared to the number of fragile regions (for example, in the yeast genomes that we analyze in the Discussion section, the number of chromosomes ranges between 6 and 8, while the number of fragile regions is at least 710), these artificial edges do not significantly affect the estimation. We refer to [7] for a discussion of subtle differences in the analysis of circular and linear genomes.
Results
We propose a new method for estimating the evolutionary distance between two genomes with high accuracy under the FBM. A key component of our method is the analytical estimation of c _{ ℓ } for various ℓ≥2 in terms of k and n, which we describe in the Theoretical analysis subsection below. From this estimation, in the Estimation algorithm subsection we solve the inverse problem to find the true evolutionary distance k (which is our goal) and the number of fragile regions n (as a byproduct). In our analysis, we consider only relatively small ℓ and assume that n and k are sufficiently large (see Theorem 1).
Theoretical analysis
Theorem 1
To prove Theorem 1, we will need some lemmas.
Let us consider a sequence \(\mathcal {D}\) of k random DCJs transforming the graph G(P _{ n },P _{ n }) into the graph G(P _{ n },Q _{ n }). Each DCJ in \(\mathcal {D}\) either merges two cycles (merging DCJ), splits one cycle into two (splitting DCJ), or preserves the cycle structure (preserving DCJ). We call set of ℓ black edges proper for the transformation \(\mathcal {D}\) if in the breakpoint graph G(P _{ n },Q _{ n }) these edges form an ℓcycle resulted entirely from ℓ−1 merging DCJs. We remark that each ℓcycle in G(P _{ n },Q _{ n }) is a result of ℓ−1 merging DCJs (we denote the number of such cycles by \(\tilde {c}_{\ell }\)) or is formed with at least one splitting or preserving DCJ. The following lemmas provide an asymptotic for \(\tilde {c}_{\ell }\) and an upper bound for \(c_{\ell }  \tilde {c}_{\ell }\).
Lemma 1
Proof
Let us fix a set A of ℓ black edges and find the probability that A is proper for \(\mathcal {D}\). We remark that this probability does not depend on the content of A but only on its size A=ℓ, and denote it by \(\tilde {p}_{n,k,\ell }\). There are 2^{ ℓ−1}(ℓ−1)! ways to arrange ℓ black edges into an ℓcycle, since there exist (ℓ−1)! circular permutations of length ℓ and 2^{ ℓ−1} ways to assign directions to their elements. For any ℓcycle, there are ℓ ^{ ℓ−2} ways to obtain it as a result of ℓ−1 DCJs [8].
□
Lemma 2
The variance of \(\tilde {c}_{\ell }\) is bounded: \(\text {Var} \left (\tilde {c}_{\ell }\right) = O(n)\).
Proof
For a positive integer M, we call a splitting DCJ Msplitting, if it splits some mcycle into an icycle and an (m−i)cycle for some i≤M. Similarly, we call a preserving DCJ Mpreserving if it operates on an mcycle for some m≤M.
Lemma 3
Let M be any positive integer. Then (i) the number of Msplitting DCJs in \(\mathcal {D}\) is stochastically dominated by a Poisson random variable with parameter γ M/2; (ii) the number of Mpreserving DCJs in \(\mathcal {D}\) is stochastically dominated by a Poisson random variable with parameter γ M ^{2}/4.
Proof
Since \(e^{\gamma M/2} \frac {(\gamma M/2)^{j}}{j!}\) represents the probability that a Poisson random variable with parameter γ M/2 is equal to j, the proof of (i) is completed.
Now we can prove Theorem 1.
Proof
Lemma 3 implies that \(c_{\ell }  \tilde {c}_{\ell }\), the number of ℓcycles obtained with at least one splitting or preserving DCJ, is of order O(1). Indeed, each such ℓcycle uniquely determines the last splitting or preserving DCJ that participates in the formation of this cycle (i.e., operates on some of the cycle vertices). Clearly, this last splitting (resp. preserving) DCJ is ℓsplitting (resp. ℓpreserving). Since such last DCJ correspond to at most two cycles, it follows that the number \(c_{\ell }  \tilde {c}_{\ell }\) does not exceed twice the number of ℓsplitting and ℓpreserving DCJs. By Lemma 3, the number of ℓsplitting and ℓpreserving DCJs is bounded (stochastically dominated by a Poisson random variable), we have \(c_{\ell }  \tilde {c}_{\ell } = O(1)\).
From Theorem 1, we have
Corollary 1
Proof
Indeed, b=n−c _{1} and \(\frac {c_{1}}{n} = e^{\gamma }+o_{p}\,(1)\,\). □
Corollary 2
Proof
Estimation algorithm
Our simulations demonstrate that this estimation method is robust for γ<2. For larger values of γ, the estimator is too sensitive to small random deviations.
We report a value of k _{ e } as our estimation for the true evolutionary distance.
Discussion
We evaluated the proposed method on simulated genomes, and further applied it for estimation of the evolutionary distances within a set of five yeast genomes and a set of two fish genomes.
Simulated genomes
Since the function in the r.h.s. is continuous and increasing, a unique solution exists for any value of the l.h.s. From the estimated value γ _{ e }, we compute n _{ e } and k _{ e } as in (5).
Instead of having 6 terms in the denominator of (6), one can take a smaller number m≥2. While the behavior of the sum with a larger m is more stable generally, our simulations showed very close results for all m between 2 and 10.
Instead of (6), we also used the approximate Eq. (3), which involves only two observable parameters: the DCJ distance d and half the number of breakpoints b. While the resulting estimation is quite similar to the previous one, it is more accurate when γ<1.6 and less accurate γ>2.
Yeast genomes
For each pair of yeast genomes, x : y gives the estimated true evolutionary distance x and the rearrangement distance y
K. lactis  K. thermotolerans  S. kluyveri  Z. rouxii  

A. gossypii  375 : 359  260 : 249  228 : 217  336 : 319 
K. lactis  280 : 272  253 : 240  352 : 344  
K. thermotolerans  75 : 75  198 : 198  
S. kluyveri  162 : 162 
Fish genomes
We also analyzed two fish genomes: Tetraodon nigroviridis [12] and Gasterosteus aculeatus [13] represented as sequences of the same 6132 genes. Our estimation for the true evolutionary distance between these genomes shows that they do not fall under the parsimony phase: their rearrangement distance equals 3705, while the true evolutionary distance is about 4500.
Conclusions
In the current study, we address the problem of estimating the true evolutionary distance between two genomes under the fragile breakage model. Similarly to our previous study on estimation of the proportion of evolutionary transpositions [14], we model evolution as a sequence of random DCJs and track how these random DCJs change the cycle structure of the breakpoint graph. We show that, while the number of DCJs is less than half the number of fragile regions in the genome, the parsimony assumption holds, and in this case we prove that lengths of alternating cycles are distributed according to a Borel distribution. In this sense our process in the parsimony phase is closely related to the evolution of a random Erdös–Rényi graph [10] in the subcritical regime. A similar process was also analyzed by Berestycki and Durrett [15]. They studied the cycle structure of a permutation obtained from the identity permutation with a sequence of random algebraic transpositions. Our results are consistent with their findings.
We provide estimators for the true evolutionary distance, which show high accuracy on simulated genomes. Our analysis of five yeast genomes shows that some but not all genome pairs fall under the parsimony phase. We also analyzed two fish genomes and revealed that the rearrangement distance between them underestimate the true evolutionary distance by about 20%. This data show how drastically the two distances can differ and emphasize the importance of using the evolutionary (rather than rearrangement) distance in comparative genomics. In contrast to the method of Lin and Moret [2], our method does not rely on the number of common gene adjacencies across two given genomes. Since some genes can form conserved clusters [16], treating all gene adjacencies as fragile regions can lead to a huge bias in the estimation. Our estimator is based only on breakpoints (but not on conserved adjacencies) across two genomes, and so it avoids this issue.
We further remark that our method is based on accurate estimation of the ratio \(\frac {d}{b}\). At the same time, in comparative genomics studies the ratio \(\frac {2d}{b}\) is known as the breakpoint reuse rate [5, 17], representing the average number of breakpoint “uses” by rearrangements in the course of evolution between two genomes. Our parsimony phase condition can therefore be restated as the breakpoint reuse rate being below 1.58. We argue however that the conventional definition of the breakpoint reuse rate does not accurately capture its meaning, in the same way how the rearrangement distance does not quite capture the meaning of evolutionary distance. This inspires a notion of the true breakpoint reuse rate defined as \(\frac {2k}{b}\) (i.e., with the true evolutionary instead of rearrangement distance in the numerator), which can be easily estimated with our method.
In further development of our method, we plan to address the even more accurate turnover fragile breakage model [6], where the set of fragile regions changes with time. We believe this model can better explain the cycle structure of the breakpoint graphs of real genomes (for a recent discussion of associated combinatorial aspects, see [18]).
Endnotes
^{1} In the breakpoint graph constructed on synteny blocks of two genomes, there are no trivial cycles since no adjacency is shared by both genomes. However, the breakpoint graph constructed on orthologous genes or multigenome synteny blocks may contain trivial cycles.
^{2} In contrast, the method of Lin and Moret [2] considers only the latter option, which corresponds to the RBM.
^{3} In contrast to c _{1}(P _{ n },Q _{ n }), the value of c _{1}(P,Q) is rather arbitrary and thus is ignored in our model.
^{4} We remind that a sequence of random variables {X _{ n }} is stochastically bounded by a deterministic sequence {a _{ n }}, denoted X _{ n }=O _{ p }(a _{ n }), if for all ε>0, there exists C such that for all n, Pr{X _{ n }/a _{ n }>C}<ε. Similarly, X _{ n }=o _{ p }(a _{ n }) if for all C>0, \({\lim }_{n\to \infty } \text {Pr}\{\left X_{n}/a_{n}\right  > C\} =0\).
Declarations
Acknowledgements
A short abstract of the present work appeared in the proceedings of the 5th IEEE International Conference on Computational Advances in Bio and Medical Sciences (ICCABS 2015).
Funding
The work is supported by the National Science Foundation under the grant No. IIS1462107, which has also funded the publication costs.
Availability of data and materials
The data is available upon request.
Authors’ contributions
Both authors wrote the manuscript, designed the experiments, and analyzed the data. NA also run the experiments and prepared the illustrations. All authors read and approved the final manuscript.
Competing interests
The authors declare that they have no competing interests.
Consent for publication
Not applicable.
Ethics approval and consent to participate
Not applicable.
About this supplement
This article has been published as part of BMC Genomics Volume 18 Supplement 4, 2017: Selected articles from the Fifth IEEE International Conference on Computational Advances in Bio and Medical Sciences (ICCABS 2015): Genomics. The full contents of the supplement are available online at https://bmcgenomics.biomedcentral.com/articles/supplements/volume18supplement4.
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
Authors’ Affiliations
References
 Yancopoulos S, Attie O, Friedberg R. Efficient sorting of genomic permutations by translocation, inversion and block interchange. Bioinformatics. 2005; 21(16):3340–346. doi:10.1093/bioinformatics/bti535.View ArticlePubMedGoogle Scholar
 Lin Y, Moret BME. Estimating true evolutionary distances under the DCJ model. Bioinformatics. 2008; 24(13):114–22. doi:10.1093/bioinformatics/btn148.View ArticleGoogle Scholar
 Ohno S. Evolution by Gene Duplication. New York, NY: Springer; 1970.View ArticleGoogle Scholar
 Nadeau JH, Taylor BA. Lengths of Chromosomal Segments Conserved since Divergence of Man and Mouse. Proc Nat Acad Sci. 1984; 81(3):814–18. doi:10.1073/pnas.81.3.814.View ArticlePubMedPubMed CentralGoogle Scholar
 Pevzner PA, Tesler G. Human and mouse genomic sequences reveal extensive breakpoint reuse in mammalian evolution. Proc Nat Acad Sci. 2003; 100:7672–677. doi:10.1073/pnas.1330369100.View ArticlePubMedPubMed CentralGoogle Scholar
 Alekseyev MA, Pevzner PA. Comparative Genomics Reveals Birth and Death of Fragile Regions in Mammalian Evolution. Genome Biol. 2010; 11(11):117. doi:10.1186/gb20101111r117.View ArticleGoogle Scholar
 Alekseyev MA. Multibreak rearrangements and breakpoint reuses: from circular to linear genomes. J Comput Biol. 2008; 15(8):1117–131. doi:10.1089/cmb.2008.0080.View ArticlePubMedGoogle Scholar
 Ouangraoua A, Bergeron A. Combinatorial structure of genome rearrangements scenarios. J Comput Biol. 2010; 17(9):1129–144. doi:10.1089/cmb.2010.0126.View ArticlePubMedGoogle Scholar
 Tanner J. A derivation of the Borel distribution. Biometrika. 1961; 48(1/2):222–24.View ArticleGoogle Scholar
 Erdös P, Rényi A. On the evolution of random graphs. Publ Math Inst Hung Acad Sci. 1960; 5:17–61.Google Scholar
 Tannier E. Yeast Ancestral Genome Reconstructions: The Possibilities of Computational Methods. In: Ciccarelli FD, Miklós I, editors. Proceedings of the 7th Annual Research in Computational Molecular Biology Satellite Workshop on Comparative Genomics (RECOMBCG). Lecture Notes in Computer Science. Berlin Heidelberg: Springer: 2009. p. 1–12. doi:10.1007/9783642047442_1.Google Scholar
 Jaillon O, Aury JM, Brunet F, Petit JL, StangeThomann N, Mauceli E, Bouneau L, Fischer C, OzoufCostaz C, Bernot A, et al.Genome duplication in the teleost fish tetraodon nigroviridis reveals the early vertebrate protokaryotype. Nature. 2004; 431(7011):946–57. doi:10.1038/nature03025.View ArticlePubMedGoogle Scholar
 Jones FC, Grabherr MG, Chan YF, Russell P, Mauceli E, Johnson J, Swofford R, Pirun M, Zody MC, White S, et al.The genomic basis of adaptive evolution in threespine sticklebacks. Nature. 2012; 484(7392):55–61. doi:10.1038/nature10944.View ArticlePubMedPubMed CentralGoogle Scholar
 Alexeev N, Aidagulov R, Alekseyev MA. A computational method for the rate estimation of evolutionary transpositions. In: Ortuño F, Rojas I, editors. Proceedings of the 3rd International WorkConference on Bioinformatics and Biomedical Engineering (IWBBIO). Lecture Notes in Computer Science. Switzerland: Springer: 2015. p. 471–80. doi:10.1007/9783319164830_46.Google Scholar
 Berestycki N, Durrett R. A phase transition in the random transposition random walk. Probab Theory Relat Fields. 2006; 136(2):203–33. doi:10.1007/s0044000504797.View ArticleGoogle Scholar
 SpringPearson SM, Stone JK, Doyle A, Allender CJ, Okinaka RT, Mayo M, Broomall SM, Hill JM, Karavis MA, Hubbard KS, et al.Pangenome analysis of burkholderia pseudomallei: Genome evolution preserves gene order despite high recombination rates. PloS One. 2015; 10(10):0140274. doi:10.1371/journal.pone.0140274.View ArticleGoogle Scholar
 Attie O, Darling AE, Yancopoulos S. The rise and fall of breakpoint reuse depending on genome resolution. BMC Bioinformatics. 2011; 12(Suppl 9):1. doi:10.1186/1471210512S9S1.View ArticleGoogle Scholar
 Alexeev N, Pologova A, Alekseyev MA. Generalized Hultman Numbers and Cycle Structures of Breakpoint Graphs. J Comput Biol. 2017; 24(2):93–105. doi:10.1089/cmb.2016.0190.View ArticlePubMedGoogle Scholar