Estimation of the True Evolutionary Distance
under the Fragile Breakage Model^{†}^{†}thanks: The work is supported by the National Science Foundation under the grant No. IIS1462107.
Washington, DC, USA
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 . These results demonstrate how drastically the two distances can differ and justify the use of true evolutionary distance in phylogenomic studies.
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 blocks as a genome graph composed of directed block edges encoding blocks and their strands, and undirected adjacency edges encoding adjacencies between blocks.
Let and be genomes on the same set of blocks. We assume that in their genome graphs, the adjacency edges of are colored black (Fig. 1a) and the adjacency edges of are colored red (Fig. 1b). The breakpoint graph is the superposition of the genome graphs of and with the block edges removed (Fig. 1c). The black and red adjacency edges in form a collection of alternating blackred cycles.
We say that a blackred cycle is an cycle if it contains black edges (and red edges) and let be the number of cycles in . We refer to cycles as trivial^{1}^{1}1In 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. and to the other cycles as nontrivial. The vertices of nontrivial cycles are called breakpoints.
A DCJ in genome replaces any pair of red adjacency edges , with either a pair of edges , or a pair of edges , . We say that such a DCJ operates on the edges , and their endpoints . A DCJ in genome transforming it into a genome corresponds to a transformation of the breakpoint graph into the breakpoint graph (Fig. 2). Each DCJ in the breakpoint graph can merge two blackred cycles into one (if edges , belong to distinct cycles), split one cycle into two or keep the number of cycles intact (if edges , belong to the same cycle). The DCJ distance between genomes and is the minimum number of DCJs required to transform into . It can be evaluated as , where is half the number of breakpoints and is the number of nontrivial cycles in the breakpoint graph [1].
Evolutionary Model
To estimate the true evolutionary distance between genomes and on the same set of blocks, we view the evolution between them as a discrete Markov process that transforms genome into genome with a sequence of DCJs occurring independently. The process starts at genome and ends at , and corresponds to a transformation of the breakpoint graphs starting at (formed by a collection of trivial cycles) and ending at . The number of DCJs in this transformation represents the true evolutionary distance between genomes and .
We remark that under the FBM, the number of trivial cycles (if any) in 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 and ) by chance.^{2}^{2}2In contrast, the method of Lin and Moret [2] considers only the latter option, which corresponds to the RBM. Furthermore, there may exist such conserved fragile regions within the blocks forming and and thus such regions cannot be observed in . To account for all these possibilities, we assume that and are composed of a large unknown number 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 and are formed by (invisible) solid regions and conserved fragile regions. Let and denote these representations of and as sequences of the solid regions. We view genome as obtained from with a sequence of 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 of solid regions, the breakpoint graphs and have the same cycle structure, except for trivial cycles. That is, we have for every , implying, in particular, that , , and . Indeed, if genomes and are obtained from and by replacing a single block with two consecutive smaller blocks , then can be obtained from by adding one trivial cycle (corresponding to the shared adjacency ). Since the genomes and can be obtained from and with a number of such operations, the breakpoint graphs and may differ only in the number of trivial cycles.^{3}^{3}3In contrast to , the value of is rather arbitrary and thus is ignored in our model.
In our evolutionary model, the following parameters are observable:

for any , i.e., the number of cycles in ;

, the number of broken fragile regions between and , which is also the number of synteny blocks between and , or half of the total length of all nontrivial cycles in ;

, the DCJ distance between and ;
while the following parameters are hidden:

, the number of trivial cycles in ;

, the number of fragile regions in each of genomes and , half the total length of all cycles in ;

, the number of DCJs in the Markov process, the true evolutionary distance between and .
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 and , while the number of fragile regions is at least ), 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 for various in terms of and , 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 (which is our goal) and the number of fragile regions (as a byproduct). In our analysis, we consider only relatively small and assume that and are sufficiently large (see Theorem 1).
Theoretical Analysis
Theorem 1.
Let genome be a genome with fragile regions and genome be obtained from with random DCJs for some . Then, for any fixed , the proportion of edges that belong to cycles in is
where denotes a term stochastically bounded^{4}^{4}4We remind that a sequence of random variables is stochastically bounded by a deterministic sequence , denoted , if for all , there exists such that for all , . Similarly, if for all , ..
To prove Theorem 1, we will need some lemmas.
Let us consider a sequence of random DCJs transforming the graph into the graph . Each DCJ in 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 if in the breakpoint graph these edges form an cycle resulted entirely from merging DCJs. We remark that each cycle in is a result of merging DCJs (we denote the number of such cycles by ) or is formed with at least one splitting or preserving DCJ. The following lemmas provide an asymptotic for and an upper bound for .
Lemma 1.
The mean value of has the following asymptotic:
Proof.
Let us fix a set of black edges and find the probability that is proper for . We remark that this probability does not depend on the content of but only on its size , and denote it by . There are ways to arrange black edges into an cycle, since there exist circular permutations of length and ways to assign directions to their elements. For any cycle, there are ways to obtain it as a result of DCJs [8].
We represent as the union , where the subsequence contains DCJs connecting edges from into an cycle, and the subsequence contains the rest of the DCJs. Since and , there are ways to choose positions for elements of in . The DCJs from operate on the red edges that are not incident to any black edge from , and for each pair of red edges there are two possible ways to recombine them with a DCJ. This gives
possible subsequences . So, there are
transformations such that is proper for . The total number of step transformations is equal to Then the probability can be found as follows
(1) 
Since there are ways to choose the set of black edges, the average number of cycles obtained by only merging DCJs equals
(2) 
∎
Lemma 2.
The variance of is bounded: .
Proof.
For a transformation and a set of black edges, we define a random variable :
By definition, where the sum is taken over all possible sets of black edges. The variance of is equal to
Since
we will find for each pair of sets of black edges:
where is defined by (1). This implies
Finally, using (2), we obtain
∎
For a positive integer , we call a splitting DCJ splitting, if it splits some cycle into an cycle and an cycle for some . Similarly, we call a preserving DCJ preserving if it operates on an cycle for some .
Lemma 3.
Let be any positive integer. Then (i) the number of splitting DCJs in is stochastically dominated by a Poisson random variable with parameter ; (ii) the number of preserving DCJs in is stochastically dominated by a Poisson random variable with parameter .
Proof.
To prove (i), we notice that the probability that a DCJ from splits a fixed cycle into an cycle and an cycle is (if ) or (if ). The probability that a DCJ splits a cycle into an cycle with and another cycle can be bounded by
This bound implies that a number of splitting DCJs in is stochastically dominated by the random variable that equals with the probability:
Since represents the probability that a Poisson random variable with parameter is equal to , the proof of (i) is completed.
To prove (ii), we notice that the probability that a DCJ from operates on two red edges from the fixed cycle and does not split this cycle equals . Summing over all , we bound the probability that a fixed DCJ is preserving as
and the proof is completed similarly to the proof of (i) above. ∎
Now we can prove Theorem 1.
Proof of Theorem 1..
Lemma 3 implies that , the number of cycles obtained with at least one splitting or preserving DCJ, is of order . 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 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 .
We remark that for , the sequence
defines a probability mass function, which characterizes a Borel distribution [9]. So, if , for a randomly chosen edge, the length of the cycle that contains it follows a Borel distribution with parameter . Moreover, the results of Erdös and Rényi [10] imply that
which can also be seen empirically in Fig. 3a. It follows that for , and the DCJ distance closely approximates the true distance. While corresponds to the parsimony assumption, we say that the process is in the parsimony phase as soon as .
From Theorem 1, we have
Corollary 1.
Proof.
Indeed, and . ∎
Corollary 2.
Proof.
By definition,
For any fixed , the number of cycles, where , is bounded by . By Theorem 1, for any fixed , we have , where is some random variable with and . Hence, for any fixed we have:
Let , then
From the definition of it follows that the random variable has the mean value and the variance . Let
so that . Then
The lower bound is obvious
Since the series converges, for each one can choose and then choose in such a way that
∎
The estimations given in Corollaries 1 and 2 are very precise, as shown in Fig. 3. All the simulations are performed for (see Discussion).
(a)  (b)  



Estimation Algorithm
From Corollaries 1 and 2, we obtain the following approximation for the ratio :
(3) 
and then estimate with the bisection method. The plot of as a function of is shown in Fig. 4a. As one can see, this function is increasing, and each value of in the interval uniquely determines the value of . In particular, corresponds to . That is, if , then the process is in the parsimony phase and the true distance is accurately approximated by the DCJ distance .
Our simulations demonstrate that this estimation method is robust for . For larger values of , the estimator is too sensitive to small random deviations.
(a)  (b)  



Alternatively, the value of can be approximated by for some . In our simulations, the best results were observed for . Again, the function
(4) 
is increasing (Fig. 4b). The applicability limits for this estimator depend on the value of . Namely, if is close to zero, then the estimator becomes sensitive to small deviations. In our simulations with (see Discussion), the estimator is very precise for values with the relative error below 10% in 95% of observations.
Once we obtain an estimated value for , it is easy to estimate the values of and as follows:
(5) 
We report a value of 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
We performed a simulation with a fixed number of blocks and various values of . In each simulation, we started with a genome on blocks and applied a number of DCJs, until we reached the upper value of (in our case this upper value is 2.5). We denote the resulting genome by and estimate , , and from the genomes and . First we estimate by solving the approximate equation
(6) 
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 , we compute and as in (5).
The result of this procedure is shown in Fig. 5 with a boxplot diagram of relative error of our estimation for each with a step . As one can see, the relative error increases with the increase of , but even for large values of , the median of the relative error is small (e.g., for the median of the relative error is only ); and for all , the interquartile range is less than .
Instead of having terms in the denominator of (6), one can take a smaller number . While the behavior of the sum with a larger is more stable generally, our simulations showed very close results for all between and .
Yeast genomes
We analyzed a set of five yeast genomes: A. gossypii, K. lactis, K. thermotolerans, S. kluyveri, and Z. rouxii, represented as sequences of the same synteny blocks [11]. For each pair of genomes, we circularized their chromosomes, constructed the breakpoint graph, and independently estimated the true evolutionary distance between them. The results in Table 1 demonstrate that some but not all pairs of yeast genomes fall under the parsimony phase.
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 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 , while the true evolutionary distance is about .
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 . 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 . At the same time, in comparative genomics studies the ratio 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 . 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 (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]).
References
 [1] Yancopoulos, S., Attie, O., Friedberg, R.: Efficient sorting of genomic permutations by translocation, inversion and block interchange. Bioinformatics 21(16), 3340–3346 (2005). doi:10.1093/bioinformatics/bti535
 [2] Lin, Y., Moret, B.M.E.: Estimating true evolutionary distances under the DCJ model. Bioinformatics 24(13), 114–122 (2008). doi:10.1093/bioinformatics/btn148
 [3] Ohno, S.: Evolution by Gene Duplication. Springer, New York, NY (1970)
 [4] Nadeau, J.H., Taylor, B.A.: Lengths of Chromosomal Segments Conserved since Divergence of Man and Mouse. Proceedings of the National Academy of Sciences 81(3), 814–818 (1984). doi:10.1073/pnas.81.3.814
 [5] Pevzner, P.A., Tesler, G.: Human and mouse genomic sequences reveal extensive breakpoint reuse in mammalian evolution. Proceedings of the National Academy of Sciences 100, 7672–7677 (2003). doi:10.1073/pnas.1330369100
 [6] Alekseyev, M.A., Pevzner, P.A.: Comparative Genomics Reveals Birth and Death of Fragile Regions in Mammalian Evolution. Genome Biology 11(11), 117 (2010). doi:10.1186/gb20101111r117
 [7] Alekseyev, M.A.: Multibreak rearrangements and breakpoint reuses: from circular to linear genomes. Journal of Computational Biology 15(8), 1117–1131 (2008). doi:10.1089/cmb.2008.0080
 [8] Ouangraoua, A., Bergeron, A.: Combinatorial structure of genome rearrangements scenarios. Journal of Computational Biology 17(9), 1129–1144 (2010). doi:10.1089/cmb.2010.0126
 [9] Tanner, J.: A derivation of the Borel distribution. Biometrika, 222–224 (1961)
 [10] Erdös, P., Rényi, A.: On the evolution of random graphs. Publ. Math. Inst. Hung. Acad. Sci 5, 17–61 (1960)
 [11] Tannier, E.: Yeast Ancestral Genome Reconstructions: The Possibilities of Computational Methods. In: Ciccarelli, F.D., Miklós, I. (eds.) Proceedings of the 7th Annual Research in Computational Molecular Biology Satellite Workshop on Comparative Genomics (RECOMBCG). Lecture Notes in Computer Science, vol. 5817, pp. 1–12. Springer, Berlin Heidelberg (2009). doi:10.1007/9783642047442_1
 [12] Jaillon, O., Aury, J.M., Brunet, F., Petit, J.L., 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 431(7011), 946–957 (2004). doi:10.1038/nature03025
 [13] Jones, F.C., Grabherr, M.G., Chan, Y.F., Russell, P., Mauceli, E., Johnson, J., Swofford, R., Pirun, M., Zody, M.C., White, S., et al.: The genomic basis of adaptive evolution in threespine sticklebacks. Nature 484(7392), 55–61 (2012). doi:10.1038/nature10944
 [14] Alexeev, N., Aidagulov, R., Alekseyev, M.A.: A computational method for the rate estimation of evolutionary transpositions. In: Ortuño, F., Rojas, I. (eds.) Proceedings of the 3rd International WorkConference on Bioinformatics and Biomedical Engineering (IWBBIO). Lecture Notes in Computer Science, vol. 9043, pp. 471–480. Springer, Switzerland (2015). doi:10.1007/9783319164830_46
 [15] Berestycki, N., Durrett, R.: A phase transition in the random transposition random walk. Probability theory and related fields 136(2), 203–233 (2006). doi:10.1007/s0044000504797
 [16] SpringPearson, S.M., Stone, J.K., Doyle, A., Allender, C.J., Okinaka, R.T., Mayo, M., Broomall, S.M., Hill, J.M., Karavis, M.A., Hubbard, K.S., et al.: Pangenome analysis of burkholderia pseudomallei: Genome evolution preserves gene order despite high recombination rates. PloS one 10(10), 0140274 (2015). doi:10.1371/journal.pone.0140274
 [17] Attie, O., Darling, A.E., Yancopoulos, S.: The rise and fall of breakpoint reuse depending on genome resolution. BMC bioinformatics 12(Suppl 9), 1 (2011). doi:10.1186/1471210512S9S1
 [18] Alexeev, N., Pologova, A., Alekseyev, M.A.: Generalized Hultman Numbers and Cycle Structures of Breakpoint Graphs. Journal of Computational Biology 24(2), 93–105 (2017). doi:10.1089/cmb.2016.0190