# Patterns of intron gain and conservation in eukaryotic genes

- Liran Carmel
^{1}, - Igor B Rogozin
^{1}, - Yuri I Wolf
^{1}and - Eugene V Koonin
^{1}Email author

**7**:192

https://doi.org/10.1186/1471-2148-7-192

© Carmel et al; licensee BioMed Central Ltd. 2007

**Received:**18 June 2007**Accepted:**12 October 2007**Published:**12 October 2007

## Abstract

### Background:

The presence of introns in protein-coding genes is a universal feature of eukaryotic genome organization, and the genes of multicellular eukaryotes, typically, contain multiple introns, a substantial fraction of which share position in distant taxa, such as plants and animals. Depending on the methods and data sets used, researchers have reached opposite conclusions on the causes of the high fraction of shared introns in orthologous genes from distant eukaryotes. Some studies conclude that shared intron positions reflect, almost entirely, a remarkable evolutionary conservation, whereas others attribute it to parallel gain of introns. To resolve these contradictions, it is crucial to analyze the evolution of introns by using a model that minimally relies on arbitrary assumptions.

### Results:

We developed a probabilistic model of evolution that allows for variability of intron gain and loss rates over branches of the phylogenetic tree, individual genes, and individual sites. Applying this model to an extended set of conserved eukaryotic genes, we find that parallel gain, on average, accounts for only ~8% of the shared intron positions. However, the distribution of parallel gains over the phylogenetic tree of eukaryotes is highly non-uniform. There are, practically, no parallel gains in closely related lineages, whereas for distant lineages, such as animals and plants, parallel gains appear to contribute up to 20% of the shared intron positions. In accord with these findings, we estimated that ancestral introns have a high probability to be retained in extant genomes, and conversely, that a substantial fraction of extant introns have retained their positions since the early stages of eukaryotic evolution. In addition, the density of sites that are available for intron insertion is estimated to be, approximately, one in seven basepairs.

### Conclusion:

We obtained robust estimates of the contribution of parallel gain to the observed sharing of intron positions between eukaryotic species separated by different evolutionary distances. The results indicate that, although the contribution of parallel gains varies across the phylogenetic tree, the high level of intron position sharing is due, primarily, to evolutionary conservation. Accordingly, numerous introns appear to persist in the same position over hundreds of millions of years of evolution. This is compatible with recent observations of a negative correlation between the rate of intron gain and coding sequence evolution rate of a gene, suggesting that at least some of the introns are functionally relevant.

## Keywords

- Intron Position
- Gain Rate
- Intron Loss
- Intron Gain
- Eukaryotic Evolution

## Background

The presence of spliceosomal introns and the concurrent splicing machinery are one of the principal distinctive features of eukaryotic genomes [1–3]. Indeed, even early branching eukaryotes that were once suspected of being intronless have been shown in recent years to posses at least a few introns [4–7]. This suggests that the evolution of introns is tightly linked to the central aspects of the evolution of eukaryotes, and it even has been proposed that introns were the driving force behind the emergence of the eukaryotic nucleus and other features of the eukaryotic cell [8, 9].

Thanks to their ubiquity and potential major role in the evolution of eukaryotes, intron evolution has drawn considerable attention [1–3, 10]. It is generally accepted that introns are units of evolution such that their presence/absence pattern is a result of stochastic processes of loss and gain. The details of these processes, however, remain elusive. Only in recent years, with the accumulation of genomic data, evolution of introns has become amenable to a systematic, genome-wide analysis. However, different attempts to accomplish such a study led to incongruent conclusions regarding the prevalence, rates, and timing of intron loss and gain during the evolution of eukaryotes [11–18]. Apparently, these discrepancies are due, primarily, to incomplete underlying evolutionary models and biased estimation techniques [10].

Recently, we have obtained more conclusive results by developing a comprehensive probabilistic model of intron evolution, and by compiling a data set that is considerably larger than any one previously used [19, 20]. The probabilistic models used so far to study intron evolution can be classified into two groups: branch-specific and gene-specific ones. The branch-specific models assume that the processes of intron gain and loss along a branch are determined only by the properties of the branch, regardless of the particular gene in question [17, 18, 21]. Conversely, gene-specific models assume that these processes are determined only by the particular gene, independent on the branch in question [12]. Obviously, in reality, the characteristics of intron gain and loss processes vary considerably both across genes and across branches. Thus, each of the models employed thus far seem to describe only one facet of a more complex reality. By contrast, our model allows for the variability of intron gain and loss characteristics with respect to both genes and branches such that any of the previously suggested models can be shown to be a special case of this comprehensive model [20]. Moreover, this model is even more realistic in that it also includes rate variability between sites, with respect to both intron gain and intron loss. In order to estimate the model parameters, we devised an expectation-maximization (EM) algorithm that can also be used to reconstruct ancestral states [22]. Combining this algorithm with the profile likelihood technique allows one, in addition to all of the above, to compute confidence intervals of the model parameters.

Here, we describe, in detail, the developed model of intron evolution and the derivation of an improved version of the EM algorithm used to estimate its parameters. Previously, we have applied this comprehensive model to a set of 391 conserved genes from 19 eukaryotic species, to investigate evolutionary trends in gene structure both at the lineage level [20] and at the gene level [19]. The results of this analysis suggest that introns invaded eukaryotic genomes at early stages of eukaryogenesis in a nearly neutral process. At early times, during periods of major transitions in the eukaryotic evolution that led to population bottlenecks, these introns seem to have vastly proliferated in the ancient genomes. Gradually, a considerable fraction of the introns appear to become involved in various cellular functions, mostly, regulation of gene expression [19].

In the present work, we focus on the process of intron gain, and address the causes of the high level of intron position conservation between eukaryotic taxa. It had been already noticed that many intron positions are shared between distant eukaryotic taxa [11, 13]. For example, plants and animals share up to 25% of the intron positions [13]. However, these findings can be explained by either remarkable conservation of ancient introns or by parallel, independent, intron gains in the same positions, or (perhaps, most likely) by a combination of both these factors. The previous analyses that have attempted to quantify the relative contributions of evolutionary conservation and parallel gain in intron position sharing have differed widely, with estimates of the extent of parallel gain ranging from nearly 0% to nearly 100% [11–13, 18, 23]. Using our probabilistic model, we developed a rigorous measure for assessing the amount of parallel gain of introns. We found that, overall, parallel gain is responsible for ~8% of the shared intron positions, with the rest due to shared ancestry. However, we also demonstrate substantial heterogeneity in the extent of parallel gain, with almost none in closely related lineages, but up to ~20% in distant ones, such as plants and unikonts. On the whole, these results support the notion that intron positions are highly conserved during evolution.

## Results and discussion

### Notation

The primary input component in the study of intron evolution consists of *G* sets of aligned protein-coding sequences of orthologous genes from *S* species. Each nucleotide in these alignments is substituted by 0 or 1 depending on whether or not an intron is present following the respective position. We allow for missing data by using a third symbol (*) to indicate lack of knowledge about the presence or absence of an intron. Consequently, every site in the alignments, called *pattern*, is a vector of length *S* over the alphabet (0,1,*) and is denoted by *ω*. Let Ω be the total number of unique patterns in the entire set of *G* alignments, and let *n*_{
gp
}be the count of the number of times pattern *ω*_{
p
}(*p* = 1, ..., Ω) is found in the multiple alignment of gene *g*. Assuming that the sites evolve independently, the set *M*_{
g
}= (*n*_{g 1}, ..., *n*_{g Ω}) fully characterizes the multiple alignment of the *g* th gene.

Let *T* be a rooted, bifurcating phylogenetic tree with *S* leaves (terminal nodes), describing the evolutionary relationships between the *S* species above. The total number of nodes in *T* is *N* = 2*S* - 1, and we index them by *t* = 0,1, ..., *N* - 1, with the convention that zero is the root node. The state of node *t* is described by the variable *q*_{
t
}, which can take the values 0 and 1 (and * in leaves). We use *V*_{
t
}for the set of all leaves such that node *t* is among their ancestors. The entire collection of leaves is, obviously, *V*_{0}. The parent node of *t* is denoted P(*t*). We use the notations ${q}_{t}^{P}$ and ${V}_{t}^{P}$ for *q*_{P(t)}and *V*_{P(t)}, respectively. Analogously, the two direct descendents of the node *t* are denoted L(*t*) and R(*t*), and we use the notations ${q}_{t}^{L}$, ${q}_{t}^{R}$, ${V}_{t}^{L}$, and ${V}_{t}^{R}$ for *q*_{L(t)}, *q*_{R(t)}, *V*_{L(t)}, and *V*_{R(t)}, respectively. The branches are indexed by the node into which they are leading, and Δ_{
t
}denotes the length (in time units) of the *t* th branch. Hereinafter we assume that the tree topology, as well as all the branch lengths Δ_{1}, ..., Δ_{N-1}, are known.

### The probabilistic model

*π*

_{ i }= Pr(

*q*

_{0}=

*i*) to describe the prior probability of the root, and

*T*

_{ ij }(

*g, t*) = Pr(

*q*

_{ t }=

*j*|${q}_{t}^{P}$ =

*i, g*) to describe the transition probability for gene

*g*along branch

*t*. In our model, we assume that this transition probability depends on both the gene and the branch, and that it takes the explicit form

Here, *η*_{
g
}and *θ*_{
g
}are non-negative parameters which determine, respectively, the intron gain and loss rates per site for gene *g*. That is, along branch *t* the gene's contribution to intron gain and retention probabilities per site is $1-{e}^{-{\eta}_{g}{\Delta}_{t}}$ and ${e}^{-{\theta}_{g}{\Delta}_{t}}$, respectively. We assume that each branch is characterized by an intrinsic *branch-specific intron gain coefficient*, *ξ*_{
t
}, as well as an intrinsic *branch-specific intron loss coefficient*, *ϕ*_{
t
}, both of which are bounded, 0 ≤ *ξ*_{
t
}, *ϕ*_{
t
}≤ 1.

*rate variable*,

*r*, which scales, for each site, the time units of the phylogenetic tree, Δ

_{ t }←

*r*·Δ

_{ t }. This rate variable is a random variable that is distributed according to a distribution function with non-negative domain and unit mean, typically, the unit-mean gamma distribution. The rate variability captures rate variations among sites of the same gene. Specifically, there are fast-evolving sites (

*r*>> 1), as well as slow-evolving ones (

*r*<< 1). In our model of intron evolution, we extend this idea by assuming that the gain and loss processes are subject to rate variability, independently of each other. Hence, a site can have any combination of gain and loss rates, for example, it can be fast to gain introns but slow to lose them. To implement this approach, we use two independent rate variables,

*r*

^{ η }and

*r*

^{ θ }, that are used to scale, for each site, the gene-specific gain rate,

*η*

_{ g }←

*r*

^{ η }·

*η*

_{ g }, and the gene-specific loss rate,

*θ*

_{ g }←

*r*

^{ θ }·

*θ*

_{ g }, respectively. We further assume that the distributions of these rate variables are independent of the genes, and are explicitly given by

Here, Γ(*x*; *λ*) is the unit-mean gamma distribution of the variable *x* with the shape parameter *λ*, *δ*(*x*) is the Dirac delta-function, and *ν* is the fraction of sites that are assumed to have zero gain rate. The existence of these *zero sites* reflects the notion that introns cannot be gained at any location within genes, but rather are preferentially inserted at specific locations, contingent on particular sequence motifs known as proto-splice sites [27–29], the density of other introns in the neighborhood, the chromatin exposure, and more. According to this interpretation, 1 - *ν* measures the density of potential intron insertion sites. Importantly, using the same value of *ν* for the entire tree does not mean that the proto-splice sites are constant throughout the evolution, or are identical for different lineages. It only means that, on the average, the fraction of potential insertion sites, whatever is their concrete nature, is similar across the lineages throughout the course of eukaryotic evolution. The incorporation of such invariant sites in a rate variability model appears natural for intron evolution and has proved beneficial also in other fields of molecular evolution [30–32]. By contrast, intron loss does not have an invariant counterpart because the assumption is that, once an intron is gained, it can always be lost. Therefore, the loss rate variable is assumed to be distributed according to a gamma distribution, which is by far the most popular distribution for describing rate variability [24, 33].

In practice, the rate distributions in (3) are rendered discrete [34]. We assume that the gain rate variable can take *K*_{
η
}discrete values ${r}_{1}^{\eta}=0,{r}_{2}^{\eta},\dots ,{r}_{{K}_{\eta}}^{\eta}$ with probabilities ${f}_{1}^{\eta}=\nu ,{f}_{2}^{\eta}\dots ,{f}_{{K}_{\eta}}^{\eta}$ such that ${\sum}_{k=1}^{{K}_{\eta}}{f}_{k}^{\eta}}=1$. Analogously, we assume that the loss rate variable can take *K*_{
θ
}discrete values ${r}_{1}^{\theta},\dots ,{r}_{{K}_{\theta}}^{\theta}$ with probabilities ${f}_{1}^{\theta},\dots ,{f}_{{K}_{\theta}}^{\theta}$ such that ${\sum}_{k=1}^{{K}_{\theta}}{f}_{k}^{\theta}}=1$. For a particular gain rate value ${r}_{k}^{\eta}$, we denote the actual gain rate ${r}_{k}^{\eta}$·*η*_{
g
}by *η*_{
kg
}. Similarly, for a particular loss rate value ${r}_{k}^{\theta}$, we denote the actual loss rate ${r}_{k}^{\theta}$·*θ*_{
g
}by *θ*_{
kg
}.

For notational clarity, we aggregate the model's parameters into a small number of sets. To this end, let Ξ_{
t
}= {*ξ*_{
t
}, *ϕ*_{
t
}} be the set of parameters that are specific to branch *t*, and let Ξ = (Ξ_{1}, ..., Ξ_{N-1}) be the set of all branch-specific parameters. Similarly, let Ψ_{
g
}= (*η*_{
g
}, *θ*_{
g
}) be the set of parameters that are specific for gene *g*, and let Ψ = (Ψ_{1}, ..., Ψ_{
G
}) be the set of all gene-specific parameters. Additionally, we denote by Λ = (*π*_{0}, *ν*, *λ*_{
η
}, *λ*_{
θ
}) the "global" parameters that determine the rate variability and the prior probability of an intron absent in the root. When the distinction between the different sets of parameters is irrelevant, we shall use Θ = (Ξ, Ψ, Λ) as the set of all the model's parameters. We achieve further succinctness in notations by denoting the actual gene-specific rate values for particular values ${r}_{k}^{\eta}$ and ${r}_{k\text{'}}^{\theta}$ of the rate variables as Ψ_{kk'g}= (*η*_{
kg
}, *θ*_{k'g}).

### The Expectation-Maximization algorithm

We estimate the parameters of the model using maximum likelihood. As the probability model includes observed random variables (state of the tree leaves) as well as hidden random variables (state of internal nodes in the tree and value of actual rate variables), the expectation-maximization algorithm is a natural tool to use [35]. This is a hill climbing iterative algorithm that requires two steps in each iteration – the expectation (E) step followed by the maximization (M) step. The details of the algorithm are given under Methods.

The total number of parameters in the model is 2*G* + 4*S*, where *G* is the number of genes and *S* is the number of species. For data sets in the hundreds of genes, this sums up to >1000 parameters. For infinite data, maximum likelihood estimators are known to be nonbiased and efficient. However, as each gene typically goes through a small number of intron-related events during its lifetime, the information content of our data is limited, and cannot straightforwardly support the estimation of such a large number of parameter. To overcome this, we adopt a two-phase approach in the data analysis. In the first phase, that we denoted *homogeneous evolution*, it is assumed that all the genes have identical intron gain and loss rates, formally, that *θ*_{
g
}= *θ*_{0} and *η*_{
g
}= *η*_{0} for any gene *g*. In the second phase, denoted *heterogeneous evolution*, all the global and branch-specific parameters are fixed, which allows for estimation of the gene-specific parameters *θ*_{
g
}and *η*_{
g
}that can now take different values for different genes. Except for the rate variability within genes, the model of evolution under the homogeneous phase resembles the branch-specific models, and consequently the EM algorithm used in this part has similar structure to the EM algorithm developed by Nguyen *et al*. [18, 36].

Using simulations, we showed that this approach yields highly accurate evolutionary reconstructions: a relative error of 1%, 3%, and 11% in estimating the number of introns in internal nodes, the number of loss events along each branch, and the number of gain events along each branch, respectively [20]. The current analysis is pattern-centered rather than gene-centered and thus we have used the results of the homogeneous phase only. While slightly less accurate than the heterogeneous counterparts, the reconstructions after the homogeneous phase are still highly accurate: a relative error of 2%, 4%, and 12% in estimating the number of introns in internal nodes, the number of losses along each branch, and the number of gains along each branch, respectively [20]. Also notable is that, apart from the fraction of zero sites *v*, the rate variability within genes can be ignored without having any significant effect on the results (Ref. [20] and Additional file 1), indicating that those sites that are capable of gaining introns do so at comparable rates.

The EM algorithm was applied to the set of 391 orthologous genes from 19 eukaryotic species ([20] and see Methods), under the homogeneous evolution assumption, and the profile likelihood technique was used to estimate the 95% confidence interval for each parameter (Additional file 1). In computing the confidence intervals for a particular parameter, we allow all other parameters to vary. Sometimes, different combinations of parameters yield similar likelihood values, an observation that has already been made formal for simpler models [18]. This occasionally leads to wide single-parameter confidence intervals, especially for the intron gain and loss rates of deep nodes. Importantly, this does not reflect similarly large errors in the ensuing inference, as was demonstrated by an exhaustive simulation study [19, 20].

### One in 7 nucleotides is a potential intron insertion site

The maximum-likelihood estimate for the fraction of zero sites is *ν* = 0.862, suggesting a potential intron insertion site every 7 nucleotides. Taking into account the 95% confidence interval for *ν*, [0.631, 0.928], the density of potential insertion sites is estimated to be 1 site per 3-14 nucleotides. Based on a smaller data set, Nguyen *et al*. [18] computed a 95% confidence interval of 1 potential insertion site in 9–14 nucleotides. Although this estimate falls within our confidence interval, our results suggest, generally, a denser population of potential insertion sites. The present estimate is also compatible with the high density of intron insertion sites (one site per 9.7 nucleotides) observed for some large protein families, e.g., small G proteins [37]. The problem of estimating the number of zero sites is of fundamental importance in maximum likelihood techniques, so both Csuros [17] and Nguyen *et al*. [18] developed different heuristic procedure to handle it. We took the alternative and, arguably, more natural approach of integrating this estimate into the model as part of the gain rate variability distribution (see above) such that additional, *ad hoc* computations are not required.

One should be cautious about identifying potential intron insertion sites with proto-splice sites. More realistically, the density of potential insertion sites reflects the overall average (across species) of the impact of multiple factors such as differential tendencies to be inserted into different proto-splice sites, local densities of pre-existing introns, and degree of chromatin exposure. Furthermore, this density estimation strongly depends on the specific data set as the parameter *ν* "absorbs" the information on zero sites. Thus, if intronless genes are added to the data, only this parameter will be heavily affected.

With the total number of sites *N*_{0} = 289,902 our calculations yield 39,962 sites that can gain introns (95% confidence interval 20,891 to 106,901). Accordingly, the 5,755 sites that are actually occupied by introns in the genes analyzed here comprise 14.4% of all sites that can potentially gain introns (95% confidence interval 5.4% to 27.5%). Thus, even when data from 19 species are combined, the density of the sites actually occupied by introns is still far below the density of sites capable of hosting introns, which emphasizes the inadequacy of any analysis that considers only those sites in which introns are, actually, observed.

### The number of shared intron positions is much greater than expected by chance

*N*

_{ i }and

*N*

_{ j }be the number of intron positions in species

*i*and

*j*, respectively; let

*S*

_{ ij }be the number of intron positions shared by the two species; and let

*N*be the total number of sites capable of gaining introns. Then, approximately, we would expect that

*N*

_{ i }

*N*

_{ j }/

*N*positions will be shared between the two species by chance alone. We define the

*conservation level*

*S*

_{ ij }exceeds random expectation, whereas negative values indicate that

*S*

_{ ij }is below expectation. Even if we take for

*N*its lower 95% confidence interval bound of 20,891 (see the preceding section), almost all pairs show a positive value, namely, a greater than random number of shared positions (Fig. 1). The only exception is the pair

*S. cerevisiae*–

*P. falciparum*(the two most intron-poor species in our set) that do not have even a single shared position.

The number of positions shared by chance between two species is a random variable distributed, approximately, according to a binomial distribution with the probability of success *p* = *N*_{
i
}*N*_{
j
}/*N*^{2}, and *N* experiments. Therefore, it is easy to associate a p-value with any observed number of shared positions, *S*_{
ij
}measuring how improbable it is to obtain by chance this value or a greater one. An overall significance level of 0.05 is equivalent to a Bonfferoni-corrected significance level of 0.0003 (overly conservative as we assume that all pairs are independent). The calculations indicate that only 20 species pairs out of 171, all involving the intron-poor *S. cerevisiae* or one of the apicomplexans, had a number of shared positions that was indistinguishable from the random expectation; all other pairs had a significant excess of shared intron positions (Additional file 2).

### Twelve out of thirteen shared intron positions reflect common ancestry

The significant excess of shared intron positions over the random expectation can be explained by one of two factors or a combination thereof. The first explanation is that this observation reflects genuine evolutionary conservation; the implication is that, once an intron is gained, it is hard to lose, perhaps, due to functional importance of introns or because the deletion event itself is destructive. The second explanation is that different lineages gain introns at the same position, independently. The chance of such parallel gain is not necessarily as low as it seems at first sight because introns are, apparently, preferentially inserted into proto-splice sites (see above). Accordingly, regions of high sequence conservation might gain introns at exactly the same position.

The two explanations have opposing impacts on our understanding of the evolution of eukaryotic genes. If introns are persistent, many must have been gained early in the eukaryotic evolution, and the later evolution involved multiple losses [16, 38]. By contrast, if parallel gain is the dominant mechanism, many of the extant introns are evolutionarily young, and the eukaryotic evolution involves, primarily, multiple introns gains [12].

In order to estimate the extent of parallel gain, we examined all the sites that host introns in at least two species. Assuming that these sites are not invariant, the probability was estimated that the corresponding pattern had arisen from parallel gain, by computing the probability of the last common ancestor of the intron-bearing species to lack an intron. This approach assumes that, if the last common ancestor had an intron in a particular position, no parallel gain occurred, i.e., highly unlikely scenarios involving at least two gains and one loss at the same site are neglected. The probabilities of parallel gain for all the patterns observed in our data set are summarized in Additional file 3. Overall, the data set harbors 3176 sites with shared intron positions (741 unique patterns) out of which 317 positions (10%) are expected to result from parallel intron gain. As inferences involving the root of the tree are prone to significantly elevated standard errors [20], the calculation was repeated using only patterns that do not have the root as the last common ancestor of all intron-bearing species. This calculation yielded 2913 sites with shared intron positions (568 unique patterns) of which 229 (~7.9%) are expected to be due to parallel intron gain. Each of these calculations provides a rough estimate of the level of errors introduced into some of the recent studies that explicitly excluded the possibility of parallel gain [13, 16, 21].

Importantly, however, the contributions of parallel gains to the emergence of different patterns of intron sharing are widely different; in particular, some rare patterns are explained (almost) entirely by parallel gain and do not reflect evolutionary conservation (Additional file 3). For example, for the single site that harbors introns only in humans and *N. crassa*, the probability is >0.99 that it results from parallel gain. Considering somewhat more frequent patterns, 11 sites harbor introns in *C. intestinalis*, *A. thaliana* and *O. saliva*, with the probability of parallel gain ~0.875, and another 12 sites harbor introns only in *C. elegans* and *C. intestinalis*, with the probability of parallel gain ~0.8.

*t*(excluding the root of the tree), all patterns that have 1s in the two sub-clades stemming from

*t*, ${V}_{t}^{L}$ and ${V}_{t}^{R}$ were tallied, and the probability that

*t*is in state zero was computed. By this analysis, in which 728 unique patterns were included, we found that nearly 20% of the shared intron positions between plants and unikonts, thought to have diverged more than a billion and a half years ago [39], are due to parallel gain. In fungi and metazoa, diverged ~1.4 billion years ago, >10% of the shared positions are estimated to derive from parallel gains. In contrast, many recently diverged clades show almost no parallel gain (e.g., humans versus rodents, birds versus mammals,

*Aspergillus*versus

*Neurospora*, and flies versus mosquitoes). Table 2 lists all patterns for which the estimated contribution to parallel gain was greater than two sites. The estimated total number of sites with parallel gain is 248. Out of the 728 unique patterns, the 22 in this list (3%) account for 116 parallel gains, i.e., ~47% of the total estimated number.

The estimated number of parallel gains on the branches stemming out of each of the internal nodes in the phylogenetic tree (excluding the root; see Additional file 4)

internal node | Subclade_1 | subclade_2 | total number of shared sites | total number of parallel gains [95% confidence inrterval] | % parallel gains [95% confidence inrterval] |
---|---|---|---|---|---|

AME | Unikonts | Magnoliophyta | 630 | 122.8 [38.5 – 229.3] | 19.5 [6.1 – 36.4] |

Unikonts | Dicdi | Opisthokonts | 212 | 4.9 [1.6 – 15.3] | 2.3 [0.7 – 7.2] |

Opisthokonts | Metazoa | Fungi | 606 | 70.7 [23.0 – 123.1] | 11.7 [3.8 – 20.3] |

Metazoa | Caeel | Coelomata | 374 | 24.0 [7.8 – 38.7] | 6.4 [2.1 – 10.4] |

Coelomata | Deuterostomia | Diptera | 350 | 7.0 [2.4 – 11.5] | 2.0 [0.7 – 3.3] |

Deuterostomia | Strpu | Chordata | 1395 | 4.7 [1.6 – 8.2] | 0.3 [0.1 – 0.6] |

Diptera | Drome | Anoga | 192 | 0.0 [0.0 – 0.1] | 0.0 [0.0-0.0] |

Fungi | Cryne | Ascomycota | 223 | 7.0 [2.4 – 13.7] | 3.1 [1.1 – 6.1] |

Ascomycota | Schpo | ScAfNc | 82 | 0.3 [0.0 – 0.5] | 0.3 [0.0 – 0.7] |

ScAfNc | Sacce | Pezizomycotina | 5 | 0.0 [0.0 – 0.1] | 0.5 [0.1 – 1.2] |

Magnoliophyta | Arath | Orysa | 1337 | 0.2 [0.1 – 0.5] | 0.0 [0.0-0.0] |

Chordata | Cioin | Vertebrata | 822 | 2.5 [0.9 – 4.2] | 0.3 [0.1 – 0.5] |

Vertebrata | Danre | Amniota | 1701 | 0.3 [0.1 – 0.5] | 0.0 [0.0-0.0] |

Apicomplexa | Thepa | Plafa | 113 | 3.9 [0.4 – 8.0] | 3.4 [0.3 – 7.1] |

Pezizomycotina | Aspfu | Neucr | 221 | 0.1 [0.0 – 0.3] | 0.1 [0.0 – 0.1] |

Amniota | Galga | Mammals | 1659 | 0.1 [0.0 – 0.1] | 0.0 [0.0-0.0] |

Mammals | Homsa | Roden | 1448 | 0.0 [0.0-0.0] | 0.0 [0.0-0.0] |

Patterns that are estimated to contribute more than two sites to the total count of parallel gains

Pattern | total number of patterns | Estimated number of parallel gains |
---|---|---|

Caeel, Arath, Orysa | 20 | 14.9 |

Cryne, Arath, Orysa | 28 | 12.6 |

Cioin, Arath, Orysa | 11 | 9.6 |

Caeel, Cioin | 12 | 9.6 |

Strpu, Danre, Galga, Homsa, Arath, Orysa, Roden | 39 | 7.5 |

Strpu, Cioin, Danre, Galga, Homsa, Arath, Orysa, Roden | 32 | 6.1 |

Strpu, Cryne | 16 | 5.8 |

Strpu, Arath, Orysa | 13 | 5.7 |

Danre, Galga, Homsa, Arath, Orysa, Roden | 13 | 5.3 |

Cioin, Cryne | 5 | 4.5 |

Caeel, Cryne | 6 | 4.3 |

Strpu, Danre, Galga, Homsa, Cryne, Roden | 32 | 4.1 |

Thepa, Plafa | 65 | 3.3 |

Caeel, Thepa | 23 | 3.0 |

Caeel, Strpu, Cioin, Danre, Galga, Homsa, Arath, Orysa, Roden | 4 | 2.9 |

Aspfu, Arath | 3 | 2.6 |

Dicdi, Strpu, Danre, Galga, Homsa, Aspfu, Neucr, Roden | 4 | 2.6 |

Caeel, Strpu, Danre, Galga, Homsa, Drome, Anoga, Arath, Orysa, Roden | 7 | 2.3 |

Caeel, Strpu, Danre, Galga, Homsa, Schpo | 14 | 2.3 |

Arath, Orysa, Thepa, Plafa | 3 | 2.3 |

Cioin, Danre, Galga, Homsa, Arath, Orysa, Roden | 7 | 2.2 |

Danre, Galga, Homsa, Drome, Anoga | 3 | 2.0 |

Given that most of the other studies report overall estimates on a smaller (8 species) data set [13], a direct comparison of their results with the present ones is not feasible due to the apparent non-uniformity in the extent of parallel gain in different parts of the tree. Nevertheless, an overall parallel gain of ~10% has been computed also for the smaller data set. Rogozin *et al*. [13] used simulations to assess the extent of parallel gain and showed that, under various assumptions regarding the density of proto-splice sites, parallel gains could be responsible to ~2–40% of the shared intron positions. Accounting more accurately for the likely density of proto-splice sites, Sverdlov et al. estimated the contribution of parallel gains to the observed sharing of intron positions to be in the range of 5–10% [23]. An even higher estimate, 18.5% parallel gains, was obtained by Nguyen et al. using a branch-specific maximum-likelihood model [18].

### Intron retention

*B*denote the set of branches that comprise a path in the tree. Then, intron retention probability along this path is

where $(1-{\phi}_{t}){e}^{-{\theta}_{k\text{'}g}{\Delta}_{t}}$ is the retention probability along branch *t*, and ${f}_{k\text{'}}^{\theta}$ is the probability of being at loss-rate category *k'*.

Probabilities of an intron to be retained along selected paths in the phylogenetic tree

If an intron was present in | chances are | that it would be present in | 95% confidence interval |
---|---|---|---|

AME | 0.63 | Homsa | [0.57 – 0.68] |

Metazoa | 0.83 | Homsa | [0.79 – 0.84] |

Deuterostomia | 0.86 | Homsa | [0.85 – 0.88] |

Vertebrata | 0.95 | Homsa | [0.94 – 0.96] |

Mammals | 0.96 | Homsa | [0.95 – 0.97] |

Fungi | 0.01 | Sacce | [0.00 – 0.01] |

Fungi | 0.26 | Aspfu | [0.24 – 0.27] |

Apicomplexa | 0.26 | Plafa | [0.20 – 0.33] |

Apicomplexa | 0.69 | Thepa | [0.57 – 0.80] |

AME | 0.57 | Arath | [0.50 – 0.64] |

AME | 0.57 | Orysa | [0.50 – 0.65] |

*C. elegans*,

*H. sapiens*and

*D. melanogaster*has a probability of 0.44, 0.69, and 0.68, respectively, to have been present in the last common ancestor of metazoa.

Probability of an intron in extant species to be inherited from an ancestral node

intron is present in | Probability | that it is also present in | 95% confidence interval |
---|---|---|---|

Dicdi | 0.71 | AME | [0.46 – 0.92] |

Caeel | 0.19 | AME | [0.10 – 0.74] |

Strpu | 0.31 | AME | [0.17 – 0.79] |

Cioin | 0.22 | AME | [0.12 – 0.75] |

Danre | 0.29 | AME | [0.16 – 0.78] |

Galga | 0.29 | AME | [0.16 – 0.78] |

Homsa | 0.30 | AME | [0.17 – 0.78] |

Roden | 0.30 | AME | [0.17 – 0.78] |

Drome | 0.30 | AME | [0.17 – 0.78] |

Anoga | 0.28 | AME | [0.16 – 0.78] |

Cryne | 0.29 | AME | [0.17 – 0.78] |

Schpo | 0.39 | AME | [0.24 – 0.82] |

Sacce | 0.27 | AME | [0.16 – 0.77] |

Aspfu | 0.27 | AME | [0.15 – 0.77] |

Neucr | 0.21 | AME | [0.11 – 0.75] |

Arath | 0.31 | AME | [0.18 – 0.79] |

Orysa | 0.30 | AME | [0.17 – 0.78] |

Thepa | 0.30 | Apicomplexa | [0.10 – 0.77] |

Plafa | 0.46 | Apicomplexa | [0.21 – 0.81] |

## Conclusion

Despite the extensive attention given to the evolution of eukaryotic gene structure over the last three decades, the fundamental characteristics of this process remain controversial. In particular, depending on the methods and data sets used, different researchers have reached opposite conclusions on the causes of the remarkably high fraction of shared introns in orthologous genes from distant eukaryotic species. Some attribute it (almost) entirely to a remarkable evolutionary conservation of intron positions and others, largely, to parallel gain of introns. To resolve these contradictions, it is important to analyze the evolution of introns by using a probabilistic model that minimally relies on arbitrary assumptions. To this end, we developed a model that allows for variability of intron gain and loss rates over branches of the phylogenetic tree, individual genes, and individual sites. Applying this model to an extended set of conserved eukaryotic genes, we found that parallel gain, on average, accounts for only ~8% of the shared intron positions. However, the distribution of parallel gains over the phylogenetic tree of eukaryotes is highly non-uniform such that there are, practically, no parallel gains in closely related lineages, whereas for distant lineages, such as animals and plants, parallel gains might have contributed up to 20% of the shared intron positions. Given the distinctly non-uniform distribution of the inferred gain events over the phylogenetic tree of eukaryotes [20], most of the recently diverged lineages have amassed very few gains during their separate evolution; by contrast, deeply diverged lineages, such as animals and plants appear to have gone through independent stages of extensive intron gain, which would explain the greater share of parallel gains.

We estimated that ancestral introns have a high probability to be retained in extant genomes, and conversely, many of the extant introns are ancient. The reasons for this remarkable endurance of a substantial fraction of the introns are not clear. One possibility is mechanistic, i.e., removing existing introns might be imprecise and hence potentially deleterious. Another possibility is functional, i.e., introns might have acquired many functional roles since entering eukaryotic genomes. The latter possibility is compatible with the recent observation of a negative correlation between the rate of intron gain and coding sequence evolution rate of a gene [19] which suggests that at least some of the introns are functionally relevant.

## Methods

### The data set

The methods and criteria used in compiling the data set have been described previously [19, 20]. Briefly, the analyzed data set consisted of the reliable alignments of 391 genes from 19 eukaryotic species (a total of 289,902 sites). These include 9 animals (*Caenorhabditis elegans*, *Strongylocentrotus purpuratus*, *Ciona intestinalis*, *Danio rerio*, *Gallus gallus*, *Homo sapiens*, rodents (*Mus musculus* and *Rattus norvegicus* combined), *Drosophila melanogaster*, *Anopheles gambiae*); 5 fungi (*Cryptococcus neoformans*, *Schizosaccharomyces pombe*, *Saccharomyces cerevisiae*, *Aspergillus fumigatus*, *Neurospora crassa*); two plants (*Arabidopsis thaliana*, *Oryza sativa*); two apicomplexans (*Theileria parva*, *Plasmodium falciparum*); and the amoebozoan *Dictyostelium discoideum*.

### Phylogenetic tree topology

Throughout this paper we assumed the traditional "crown-group" tree topology (Additional file 4). Specifically, the root position is between the Apicomplexa and the common ancestor of multicellular eukaryotes (plants and animals) [40], as opposed to the alternative Unikont-Bikont division [41]. We furthermore assume the Coelomata topology (Deuterostomia and insects are grouped together to the exclusion of nematodes) [42, 43] as opposed to the Ecdysozoa topology (insects and nematodes are grouped together to the exclusion of Deuterostomia) [44, 45]. The results, however, are not sensitive to the exact tree topology, as explicitly shown elsewhere [19, 20].

### The EM Algorithm

For each site, the *S* leaves form a set of observed random variables, their states being described by the corresponding pattern *ω*_{
p
}. The states of all the internal nodes, denoted *σ*, form a set of hidden random variables, that is, random variables whose states are not observed. In order to account for rate variability across sites, we associate with each pattern two hidden random variables, ${\rho}_{p}^{\eta}$ and ${\rho}_{p}^{\theta}$, that determine the value of the rate variables in that site. To sum up, the observed random variables are *ω*_{
p
}, and the hidden random variables are (*σ*, ${\rho}_{p}^{\eta}$, ${\rho}_{p}^{\theta}$).

*L*(

*M*

_{1}, ...,

*M*

_{ G }|Θ) is guaranteed to increase as long as we maximize the auxiliary function

*w*

_{gpkk'}and

*Q*

_{gpkk'}the first and second square brackets, respectively, this expression becomes

#### The E-step

In this step, the function *Q*(Θ, Θ^{0}) or, equivalently, the set of coefficients *w*_{gpkk'}and *Q*_{gpkk'}is computed by using inward-outward recursion on the tree.

#### The inward (*γ*) recursion

*t*(except for the root) a vector ${\gamma}_{i}^{gpkk\text{'}}(t)=\mathrm{Pr}\phantom{\rule{0.25em}{0ex}}({V}_{t}|{q}_{t}^{P}=i,{\Xi}^{0},{\Psi}_{gkk\text{'}}^{0})$. This is the probability of observing the nodes

*V*

_{ t }(which are a subset of the pattern

*ω*

_{ p }) for a gene

*g*, when the gain and loss rate variables are ${r}_{k}^{\eta}$ and ${r}_{k\text{'}}^{\theta}$, respectively, and the parent node of

*t*is known to be in state

*i*. By definition, this function is initialized at all leaves (

*t*ϵ

*V*

_{0}) by

*γ*. For all internal nodes (except for the root),

*γ*is computed using the recursion

*γ*

_{ j }(L(

*t*))

*γ*

_{ j }(R(

*t*)). This is easy to see, as

The first term is, simply, the definition of *T*_{
ij
}(*g, t*). Given *q*_{
t
}, ${V}_{t}^{L}$ is independent of ${q}_{t}^{P}$, thus the second term is just Pr(${V}_{t}^{L}$ | *q*_{
t
}= *j*) = *γ*_{
j
}(*t*^{
L
}) . By similar argument, the third term is simply Pr(${V}_{t}^{R}$ | *q*_{
t
}= *j*) = *γ*_{
j
}(*t*^{
R
}). Substituting these results in (14), we recover the recursion formula (13).

*γ*-recursion allows for computing the likelihood of any observed pattern

*ω*

_{ p }, given the values of the rate variables:

*q*

_{0}, ${V}_{0}^{R}$ is independent of ${V}_{0}^{L}$, and so

#### The outward (*α*-*β*) recursion

*γ*-recursion is computed, we can use it to compute a second, complementary, recursion. To this end, let us associate with each node

*t*(except for the root node) a matrix ${\alpha}_{ij}^{gpkk\text{'}}(t)=\mathrm{Pr}\phantom{\rule{0.25em}{0ex}}({q}_{t}=j,{q}_{t}^{P}=i|{\omega}_{p},{\Xi}^{0},{\Psi}_{gkk\text{'}}^{0})$. It is convenient to define for each node

*t*(except for the root node) a vector ${\beta}_{j}^{gpkk\text{'}}(t)={\displaystyle {\sum}_{i=0}^{1}{\alpha}_{ij}^{gpkk\text{'}}(t)=}\mathrm{Pr}\phantom{\rule{0.25em}{0ex}}({q}_{t}=j|{\omega}_{p},{\Xi}^{0},{\Psi}_{gkk\text{'}}^{0})$. Upon the computation of

*α*,

*β*is readily computed, too. Again, omitting the superscripts,

*α*can be initialized from its definition on the two direct descendants of the root,

*α*is computed using the outward-recursion

*α*,

*V*

_{0}=

*V*

_{ t }+ ${\overline{V}}_{t}$, with ${\overline{V}}_{t}$ being the set of all leaves such that node

*t*is not among their ancestors. But, given ${q}_{t}^{P}$, the state of node

*t*is independent on ${\overline{V}}_{t}$, and therefore (19) becomes

*q*

_{ t },

*V*

_{ t }is independent of P(

*t*) and therefore

These inward-outward recursions are the phylogenetic equivalent of the backward-forward recursions known from hidden Markov models, and other versions of this method have been developed previously [47, 48]. The version developed here can be shown to be the realization of the junction tree algorithm [49] on rooted bifurcating trees. The junction tree algorithm is a scheme to compute marginal probabilities of maximal cliques on graphs by means of belief propagation on a modified junction tree. Indeed, the matrix *α* computes marginal probabilities of pairs (*t*, P(*t*)), but such pairs are nothing but maximal cliques on rooted bifurcating trees.

#### Computing the coefficients *w*_{
gpkk'
}

*γ*-recursion is sufficient to compute the coefficients

*w*

_{gpkk'}. From the definition, ${w}_{gpkk\text{'}}=\mathrm{Pr}\phantom{\rule{0.25em}{0ex}}({\rho}_{p}^{\eta}=k,{\rho}_{p}^{\theta}=k\text{'}|{\omega}_{p},{\Xi}^{0},{\Psi}_{g}^{0},{\Lambda}^{0})$. Using the Bayes formula Pr(

*x, y*|

*z*) = Pr(

*x, y, z*)/Σ

_{x, y}Pr(

*x, y, z*), we can rewrite it as

*w*

_{gpkk'}reduces to

The function $\mathrm{Pr}\phantom{\rule{0.25em}{0ex}}({\omega}_{p}|{\Xi}^{0},{\Psi}_{gkk\text{'}}^{0})$ is the likelihood of observing pattern *ω*_{
p
}for gain and loss rate variables ${r}_{k}^{\eta}$ and ${r}_{k\text{'}}^{\theta}$, respectively. This is readily computed upon completion of the *γ*-recursion, using (15).

#### Computing the coefficients *Q*_{gpkk'}

*α*,

*β*-recursion. By definition,

*ω*

_{ p },

*σ*| Ξ, Ψ

_{gkk'}) is the likelihood of a particular realization of the tree, thus from (1)

*δ*(

*a, b*) is the Kronecker delta function, which is 1 for

*a*=

*b*and 0 otherwise. Denote the expectation over $\mathrm{Pr}\phantom{\rule{0.25em}{0ex}}(\sigma |{\omega}_{p},{\Xi}^{0},{\Psi}_{gkk\text{'}}^{0})$ by

*E*

_{ σ }. Applying it to (25), we get

But ${E}_{\sigma}[\delta ({q}_{0},i)]=\mathrm{Pr}\phantom{\rule{0.25em}{0ex}}({q}_{0}=i|{\omega}_{p},{\Xi}^{0},{\Psi}_{gkk\text{'}}^{0})={\beta}_{i}(0)$, and similarly ${E}_{\sigma}[\delta ({q}_{t},j)\delta ({q}_{t}^{P},i)]={\alpha}_{ij}(t)$.

*Q*

_{gpkk'}is given by

### The M-step

Actually, any increase in *Q* is sufficient to guarantee an increase in the likelihood, suggesting that the precise maximization of *Q* is not particularly important. Therefore, we speed up the computations by performing low-tolerance maximization with respect to each of the parameters individually. Except for the parameters *λ*_{
η
}and *λ*_{
θ
}, it is easy to differentiate *Q* twice with respect to any parameter. This lends itself to using simple zero-finding algorithms of which we chose the Newton-Raphson algorithm [50].

Maximizing *Q* with respect to the shape parameters *λ*_{
η
}and *λ*_{
θ
}is more involved, as *Q* depends on these parameters only through the discrete approximation of the rate variability distributions (3). In our implementation, we used Yang's quantile method [34] to compute the discrete levels of the gamma distributions such that each level has equal probability. Formally, ${f}_{1}^{\eta}=\nu ,{f}_{k}^{\eta}=(1-\nu )/({K}_{\eta}-1)$ for *k* = 2, ..., *K*_{
η
}and ${f}_{k}^{\theta}=1/{K}_{\theta}$ for *k* = 1, ..., *K*_{
θ
}. To perform the maximization in this case, we used Brent's maximization algorithm that does not require derivatives [50].

#### Reconstruction of Ancestral States and Events

Given the *α*, *β*, *γ*-recursions on the tree with the final model parameters, it is straightforward to reconstruct the history of intron evolution, and to assign gains and losses to specific branches. The number of introns in an internal node *t* for a gene *g*, assuming gain rate variable ${r}_{k}^{\eta}$ and loss rate variable ${r}_{k\text{'}}^{\theta}$, given that the observed pattern is *ω*_{
p
}, is ${n}_{gp}{\omega}_{gpkk\text{'}}{\beta}_{1}^{gpkk\text{'}}(t)$. Similarly, the number of loss events along the branch *t* is ${n}_{gp}{\omega}_{gpkk\text{'}}{\alpha}_{10}^{gpkk\text{'}}(t)$, and the number of gain events along this branch is ${n}_{gp}{\omega}_{gpkk\text{'}}{\alpha}_{01}^{gpkk\text{'}}(t)$.

#### Confidence Intervals

*L*(Θ) = max

_{Θ}Pr(Θ), then the likelihood ratio

under the hypothesis that ϑ = ϑ_{0} is distributed according to the *χ*^{2} distribution with one degree of freedom. In order to find the 95% confidence interval of the parameter ϑ around its optimal value ϑ_{0}, we find the value of ϑ for which the likelihood ratio (28) exceeds the value 3.84 (95%-percentile of the *χ*^{2}(1) distribution). This value is found numerically using Ridder's method [50].

## Declarations

### Acknowledgements

We thank Jacob Goldberger for useful discussions regarding aspects of the EM algorithm. This work was supported by the Intramural Research Program of the National Institutes of Health, National Library of Medicine.

## Authors’ Affiliations

## References

- Lynch M, Richardson AO: The evolution of spliceosomal introns. Curr Opin Genet Dev. 2002, 12 (6): 701-710. 10.1016/S0959-437X(02)00360-X.View ArticlePubMedGoogle Scholar
- Roy SW, Gilbert W: The evolution of spliceosomal introns: patterns, puzzles and progress. Nat Rev Genet. 2006, 7 (3): 211-221.PubMedGoogle Scholar
- Rodríguez-Trelles F, Tarrío R, Ayala FJ: Origins and evolution of spliceosomal introns. Annu Rev Genet. 2006, 40: 47-76. 10.1146/annurev.genet.40.110405.090625.View ArticlePubMedGoogle Scholar
- Nixon JE, Wang A, Morrison HG, McArthur AG, Sogin ML, Loftus BJ, Samuelson J: A spliceosomal intron in Giardia lamblia. Proc Natl Acad Sci U S A. 2002, 99 (6): 3701-3705. 10.1073/pnas.042700299.PubMed CentralView ArticlePubMedGoogle Scholar
- Simpson AG, MacQuarrie EK, Roger AJ: Eukaryotic evolution: early origin of canonical introns. Nature. 2002, 419 (6904): 270-10.1038/419270a.View ArticlePubMedGoogle Scholar
- Vanacova S, Yan W, Carlton JM, Johnson PJ: Spliceosomal introns in the deep-branching eukaryote Trichomonas vaginalis. Proc Natl Acad Sci U S A. 2005, 102 (12): 4430-4435. 10.1073/pnas.0407500102.PubMed CentralView ArticlePubMedGoogle Scholar
- Russell AG, Shutt TE, Watkins RF, Gray MW: An ancient spliceosomal intron in the ribosomal protein L7a gene (Rpl7a) of Giardia lamblia. BMC Evol Biol. 2005, 5: 45-10.1186/1471-2148-5-45.PubMed CentralView ArticlePubMedGoogle Scholar
- Koonin EV: The origin of introns and their role in eukaryogenesis: A compromise solution to the introns-early versus introns-late debate?. Biol Direct. 2006, 1 (1): 22-10.1186/1745-6150-1-22.PubMed CentralView ArticlePubMedGoogle Scholar
- Martin W, Koonin EV: Introns and the origin of nucleus-cytosol compartmentalization. Nature. 2006, 440 (7080): 41-45. 10.1038/nature04531.View ArticlePubMedGoogle Scholar
- Rogozin IB, Sverdlov AV, Babenko VN, Koonin EV: Analysis of evolution of exon-intron structure of eukaryotic genes. Brief Bioinform. 2005, 6 (2): 118-134. 10.1093/bib/6.2.118.View ArticlePubMedGoogle Scholar
- Fedorov A, Merican AF, Gilbert W: Large-scale comparison of intron positions among animal, plant, and fungal genes. Proc Natl Acad Sci U S A. 2002, 99 (25): 16128-16133. 10.1073/pnas.242624899.PubMed CentralView ArticlePubMedGoogle Scholar
- Qiu WG, Schisler N, Stoltzfus A: The evolutionary gain of spliceosomal introns: sequence and phase preferences. Mol Biol Evol. 2004, 21 (7): 1252-1263. 10.1093/molbev/msh120.View ArticlePubMedGoogle Scholar
- Rogozin IB, Wolf YI, Sorokin AV, Mirkin BG, Koonin EV: Remarkable interkingdom conservation of intron positions and massive, lineage-specific intron loss and gain in eukaryotic evolution. Curr Biol. 2003, 13 (17): 1512-1517. 10.1016/S0960-9822(03)00558-X.View ArticlePubMedGoogle Scholar
- Roy SW, Fedorov A, Gilbert W: Large-scale comparison of intron positions in mammalian genes shows intron loss but no gain. Proc Natl Acad Sci U S A. 2003, 100 (12): 7158-7162. 10.1073/pnas.1232297100.PubMed CentralView ArticlePubMedGoogle Scholar
- Roy SW, Gilbert W: The pattern of intron loss. Proc Natl Acad Sci U S A. 2005, 102 (3): 713-718. 10.1073/pnas.0408274102.PubMed CentralView ArticlePubMedGoogle Scholar
- Roy SW, Gilbert W: Complex early genes. Proc Natl Acad Sci U S A. 2005, 102 (6): 1986-1991. 10.1073/pnas.0408355101.PubMed CentralView ArticlePubMedGoogle Scholar
- Csuros M: Likely scenarios of intron evolution. Comparative Genomics Lecture Notes in Computer Science. 2005, 3678: 47-60.View ArticleGoogle Scholar
- Nguyen HD, Yoshihama M, Kenmochi N: New maximum likelihood estimators for eukaryotic intron evolution. PLoS Comput Biol. 2005, 1 (7): e79-10.1371/journal.pcbi.0010079.PubMed CentralView ArticlePubMedGoogle Scholar
- Carmel L, Rogozin IB, Wolf YI, Koonin EV: Evolutionarily conserved genes preferentially accumulate introns. Genome Res. 2007, 17: 1045-1050. 10.1101/gr.5978207.PubMed CentralView ArticlePubMedGoogle Scholar
- Carmel L, Wolf YI, Rogozin IB, Koonin EV: Three distinct modes of intron dynamics in the evolution of eukaryotes. Genome Res. 2007, 17: 1034-1044. 10.1101/gr.6438607.PubMed CentralView ArticlePubMedGoogle Scholar
- Roy SW, Gilbert W: Rates of intron loss and gain: implications for early eukaryotic evolution. Proc Natl Acad Sci U S A. 2005, 102 (16): 5773-5778. 10.1073/pnas.0500383102.PubMed CentralView ArticlePubMedGoogle Scholar
- Carmel L, Rogozin IB, Wolf YI, Koonin EV: An expectation-maximization algorithm for analysis of evolution of exon-intron structure of eukaryotic genes. Comparative Genomics Lecture Notes in Computer Science. 2005, 3678: 35-46.View ArticleGoogle Scholar
- Sverdlov AV, Rogozin IB, Babenko VN, Koonin EV: Conservation versus parallel gains in intron evolution. Nucleic Acids Res. 2005, 33 (6): 1741-1748. 10.1093/nar/gki316.PubMed CentralView ArticlePubMedGoogle Scholar
- Nei M, Chakraborty R, Fuerst PA: Infinite allele model with varying mutation rate. Proc Natl Acad Sci U S A. 1976, 73 (11): 4164-4168. 10.1073/pnas.73.11.4164.PubMed CentralView ArticlePubMedGoogle Scholar
- Uzzell T, Corbin KW: Fitting discrete probability distributions to evolutionary events. Science. 1971, 172 (988): 1089-1096. 10.1126/science.172.3988.1089.View ArticlePubMedGoogle Scholar
- Felsenstein J: Inferring Phylogenies. 2004, Sunderland, MA , Sinauer AssociatesGoogle Scholar
- Dibb NJ: Proto-splice site model of intron origin. J Theor Biol. 1991, 151 (3): 405-416. 10.1016/S0022-5193(05)80388-1.View ArticlePubMedGoogle Scholar
- Dibb NJ, Newman AJ: Evidence that introns arose at proto-splice sites. Embo J. 1989, 8 (7): 2015-2021.PubMed CentralPubMedGoogle Scholar
- Sverdlov AV, Rogozin IB, Babenko VN, Koonin EV: Reconstruction of ancestral protosplice sites. Curr Biol. 2004, 14 (16): 1505-1508. 10.1016/j.cub.2004.08.027.View ArticlePubMedGoogle Scholar
- Gu X, Fu YX, Li WH: Maximum likelihood estimation of the heterogeneity of substitution rate among nucleotide sites. Mol Biol Evol. 1995, 12 (4): 546-557.PubMedGoogle Scholar
- Hasegawa M, Kishino H, Yano T: Man's place in Hominoidea as inferred from molecular clocks of DNA. J Mol Evol. 1987, 26 (1-2): 132-147. 10.1007/BF02111287.View ArticlePubMedGoogle Scholar
- Mayrose I, Friedman N, Pupko T: A Gamma mixture model better accounts for among site rate heterogeneity. Bioinformatics. 2005, 21 Suppl 2: ii151-ii158. 10.1093/bioinformatics/bti1125.PubMedGoogle Scholar
- Jin L, Nei M: Limitations of the evolutionary parsimony method of phylogenetic analysis. Mol Biol Evol. 1990, 7 (1): 82-102.PubMedGoogle Scholar
- Yang Z: Maximum likelihood phylogenetic estimation from DNA sequences with variable rates over sites: approximate methods. J Mol Evol. 1994, 39 (3): 306-314. 10.1007/BF00160154.View ArticlePubMedGoogle Scholar
- Dempster AP, Laird NM, Rubin DB: Maximum likelihood from incomplete data via the EM algorithm. Journal of the Royal Statistical Society, Series B. 1977, 39: 138-Google Scholar
- Nguyen HD, Yoshihama M, Kenmochi N: Authors' reply. PLoS Comput Biol. 2006, 2: e83-10.1371/journal.pcbi.0020083.PubMed CentralView ArticleGoogle Scholar
- Stoltzfus A, Logsdon JM, Palmer JD, Doolittle WF: Intron "sliding" and the diversity of intron positions. Proc Natl Acad Sci U S A. 1997, 94 (20): 10739-10744. 10.1073/pnas.94.20.10739.PubMed CentralView ArticlePubMedGoogle Scholar
- Roy SW: Intron-rich ancestors. Trends Genet. 2006, 22 (9): 468-471. 10.1016/j.tig.2006.07.002.View ArticlePubMedGoogle Scholar
- Hedges SB, Kumar S: Genomic clocks and evolutionary timescales. Trends Genet. 2003, 19: 200-206. 10.1016/S0168-9525(03)00053-2.View ArticleGoogle Scholar
- Hedges SB: The origin and evolution of model organisms. Nat Rev Genet. 2002, 3 (11): 838-849. 10.1038/nrg929.View ArticlePubMedGoogle Scholar
- Stechmann A, Cavalier-Smith T: The root of the eukaryote tree pinpointed. Curr Biol. 2003, 13 (17): R665-6. 10.1016/S0960-9822(03)00602-X.View ArticlePubMedGoogle Scholar
- Wolf YI, Rogozin IB, Koonin EV: Coelomata and not Ecdysozoa: evidence from genome-wide phylogenetic analysis. Genome Res. 2004, 14: 29-36. 10.1101/gr.1347404.PubMed CentralView ArticlePubMedGoogle Scholar
- Rogozin IB, Wolf YI, Carmel L, Koonin EV: Ecdysozoan clade rejected by genome-wide analysis of rare amino acid replacements. Mol Biol Evol. 2007, 24 (4): 1080-1090. 10.1093/molbev/msm029.View ArticlePubMedGoogle Scholar
- Aguinaldo AM, Turbeville JM, Linford LS, Rivera MC, Garey JR, Raff RA, Lake JA: Evidence for a clade of nematodes, arthropods and other moulting animals. Nature. 1997, 387 (6632): 489-493. 10.1038/387489a0.View ArticlePubMedGoogle Scholar
- Philippe H, Lartillot N, Brinkmann H: Multigene analyses of bilaterian animals corroborate the monophyly of Ecdysozoa, Lophotrochozoa, and Protostomia. Mol Biol Evol. 2005, 22 (5): 1246-1253. 10.1093/molbev/msi111.View ArticlePubMedGoogle Scholar
- Felsenstein J: Evolutionary trees from DNA sequences: a maximum likelihood approach. J Mol Evol. 1981, 17 (6): 368-376. 10.1007/BF01734359.View ArticlePubMedGoogle Scholar
- Friedman N, Ninio M, Pe'er I, Pupko T: A structural EM algorithm for phylogenetic inference. J Comput Biol. 2002, 9 (2): 331-353. 10.1089/10665270252935494.View ArticlePubMedGoogle Scholar
- Siepel A, Haussler D: Phylogenetic estimation of context-dependent substitution rates by maximum likelihood. Mol Biol Evol. 2004, 21 (3): 468-488. 10.1093/molbev/msh039.View ArticlePubMedGoogle Scholar
- Castillo E, Gutierrez JM, Hadi AS: Expert systems and probabilistic network models (Monographs in Computer Science). 1996, New York , SpringerGoogle Scholar
- Press WH, Flannery BP, Teukolsky SA, Vetterling WT: Numerical recipes in C: The art of scientific computing. 1992, New York , Cambridge University Press, 2ndGoogle Scholar

## Copyright

This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.