\section{Introduction}

Gene fragmentation and dispersal of fragments are found in various organisms
ranging from Euglenozoa (found in kinetoplastids, euglenoids, and
diplonemids) to Alveolates (found in apicomplexa, ciliates, and
dinoflagellates). See (Benne \emph{et al.}, 1986; Gillespie \emph{et al.},
1999; Marande and Burger, 2007; Kamikawa \emph{et al.}, 2007; Nowacki \emph{%
et al.}, 2008; Walker, 2007; Spencer and Gray, 2011). Gene fragmentation
occurs when the RNA is not made from a single precursor, derived directly
from the genome, but it is reconstituted from different RNA sources. This
can also happen without dispersal, \emph{i.e. }spreading the genes for these
RNA parts over the genome, intermingled with other genes.

RNA editing, first used to refer to RNA alteration processes in
trypanosomatids, entailing extreme mitochondrial (mt) gene fragment
dissemination, was discovered in 1986 (Benne \emph{et al.}, 1986).
Trypanosomatids are parasitic unicellular organisms belonging to the
kinetoplastid order. Kinetoplastids are characterized by a 'kinetoplast',
the strongly staining location of \emph{all} the mt DNA molecules within
their \emph{single} mitochondrion. Kinetoplastid RNA editing refers to
post-transcriptional sequence alteration via insertion and deletion of
uridylate residues at specific sites of mt RNAs.

Trypanosomal mt RNAs are encoded by two types of mt DNA molecules: several
thousand small minicircles and a few dozen maxicircles (Scott, 1995; Simpson
and Thiemann, 1995; Simpson \emph{et al.}, 2000; Madison-Antenucci \emph{et
al.}, 2002). Maxicircles contain ribosomal RNA genes and genes mostly
encoding subunits of respiratory chain complexes. Many subunit genes need
editing upon transcription. Both kind of circles encode guide (g) RNAs
necessary for editing of maxicircle encoded RNAs (Blum \emph{et al.}, 1990).

Very extensive editing is illustrated by the \emph{Trypanosoma brucei} cox3
gene. Cox3 is easily found (\emph{e.g.} by looking at homologies) in \emph{%
Crithidia fasciculata} and \emph{Leishmania tarentolae}: the transcript
undergoes limited editing (\emph{e.g.} in \emph{C. fasciculata} only 32 Us
are inserted and 2 Us deleted at 14 sites).

The \emph{T. brucei} cox3 'gene' was identified later: its transcript has to
be edited overall (Feagin \emph{et al.}, 1988). While the G, A and C
nucleotide sequence is maxicircle derived, the U nucleotide sequence is
generated by editing, using minicircle encoded gRNAs, with $547$ Us inserted
and $41$ Us deleted at $223$ sites. The transcript derived from the 'GAC
sequence' (called the 'cryptogene' (Simpson and Shaw, 1989)) has less than
half the edited size.

In \emph{T. brucei}, most of the mt encoded proteins come from such
cryptogenes ($9$ out of $17$, with $5$ unedited and $3$ partially edited
transcripts remaining). The extensive editing of cryptogene transcripts is
known as pan-editing (Simpson and Shaw, 1989; Sturm and Simpson, 1990). The
precise molecular mechanisms of editing, using sequential basepairing
between gRNAs and mRNA, are described in (Schnaufer \emph{et al.}, 2003) and
references therein (Madison-Antenucci \emph{et al.}, 2002; Blum \emph{et al.}%
, 1990; Feagin \emph{et al.}, 1988; Simpson and Shaw, 1989). The consequence
of pan-editing is that information necessary for production of (some)
proteins is spread and mixed over the entire mt DNA.

Here, we propose a mathematical model to support the concept (Speijer, 2006;
Speijer, 2007) that such gene fragment dispersal can function as a
biological adaptation preventing gene loss caused by intense intraspecific
competition, due to selection on growth rate in limited spaces (\emph{e.g.}
a host). The process lowers the chance that rapidly growing deletion mutants
replace \emph{all} the more complex, ecologically versatile, organisms, as
large deletions will contain linked 'active' segments. It thus works to
'preserve by association'. To put it another way: mt gene scattering
prevents the translation of short-term advantages into long-term disaster.

\section{Modelling the effects of mitochondrial gene fragment dispersal in
changing environments}

Our model describes the `life cycle' of the mt genetic material of a
trypanosomatid. In the absence of functional selection we do not distinguish
between competition at the levels of individual mt genomes (\emph{e.g.}
maxicircles) and whole organisms. At the molecular level small circles
outcompete larger circles by faster replication while at the level of the
organism, organisms carrying less mt DNA can also grow faster. \emph{At both
levels mt deletions will be selected then: they even reinforce each other.}
We need to include the following parameters to determine the relationship
between some of them:

\begin{enumerate}
\item Population size ($s_{p}$), $10^{10}$ individuals;

\item Deletion size ($l$);

\item Probability of a deletion ($p$), $10^{-6}$;

\item Level of fragmentation/editing ($k$);

\item Number of generations in each environment (in the case of parasites
`host');

\item Replication advantage of smaller DNA ($r(x)$);

\item Confidence parameter ($d$) that governs the level of approximation in
our second simulation.
\end{enumerate}

As mitochondrial DNA polymerases are much more error-prone than their
nuclear counterparts (see (Larsson, 2010) and references therein), a value
of $10^{-6}$ for the deletion probability is not unrealistic. We first model
the process of replication and deletion in its most general form, and then
simplify it, to enable us to simulate for realistic values of our
parameters. In all our simulations, we assume there are two environments, $A$
and $B$, \emph{e.g.} $A$, tsetse fly and $B$, human host. A sequence '$a$'
of base pairs of length $n_{a}$ is necessary for survival in $A$ and a
sequence '$b$' of length $n_{b}$ is needed in $B$, but not in $A$. For
simplicity we excluded sequences needed in both (violet segments in Figure
1). Total length of mt DNA will therefore be $n=n_{a}+n_{b}$ base pairs. We
chose $n_{a}$ to represent $65\%$ and $n_{b}$ $35\%$ of total mt DNA ($\sim
23000$ basepairs). We assume the parasite to be in environment $A$.

The level of mt DNA fragmentation for an individual parasite is simply
modelled as follows. The fragmentation and spreading that pan-editing of RNA
entails, intertwines subsequences of $a$ and $b$ modelled by the sequence $%
a_{1}b_{1}a_{2}b_{2}...a_{k}b_{k}$. We assume the DNA and hence this
sequence to be circular, so that $b_{k}$ connects with $a_{1}$, and we call $%
k-1$ the level of fragmentation. Initially, length of subsequence $a_{i}$ is
$n_{a}/k$ and length of subsequence $b_{i}$ is $n_{b}/k$. We model
replication rounds, in which for every round, with probability $p$, deletion
of DNA with random length $l$, at a random position, may occur. Any deletion
in an '$a$' sequence is lethal. Surviving individuals replicate such, that
the smaller their total DNA, the faster they replicate. This advantage
(given as an increased representation in the next round) is modelled by a
function $r(x)$. For example, a complete deletion of $35\%$ ($n_{b}$) of the
DNA will translate in a $3/2$ stronger representation in the next generation.

The process is modelled using a Markov chain approach, with simplifying
assumptions described in the appendix. We study the value $t_{\max }(k)$
monitoring the expected number of generations until no individuals are left
with '$b$' of length $n_{b}$. Any individual, missing part of $b$, can not
survive in the `next' environment ($B$) anymore. The increase of $t_{\max
}(k)$ gives us the added protection that additional fragmentation ($k$)
brings. At level $k$ we have $(n_{k}/k+1)k$ many states in our Markov chain
(see appendix). Realistically, $n_{b}$ is set around $8000$ base pairs and $%
k $ can be over $200$, making the Markov chain too big to handle.

By exploiting symmetries and simplifying our model we approximate our Markov
chain, obtaining a manageable model that we can simulate for realistic
values of $n_{b}$ and $k$. An additional feature of this approximation is
parameter $d$, the confidence level, which allows us to smoothly interpolate
between the simplified and the original Markov chain. Varying $d$ suggests
that already for small values we obtain accurate approximations of the
original process. See appendix and Figure 2.

\section{Results obtained by modelling the effects of increasing
mitochondrial gene fragment dispersal}

In a first simulation (see Figure 1B) we decided to compare our highly
simplified genome with and without a fragmented gene ($k=2$) under partially
relaxed selection (when only part of the genome is needed) with concomitant
strong (intraspecific) competition. We simulated a modest (\emph{e.g.} in
the bloodstream much less of the mt genome is needed) situation in which
most ($2/3$) of the mt DNA is still necessary for survival (environment $A$:
'tsetse'). Starting with $10^{10}$ individuals, we look at the number of
generations needed to reduce the amount of individual organisms that have
retained 'ecological flexibility' (\emph{i.e.} they still contain all
genetic information to be able to make the switch to another part of the
life cycle) to $1$. Without any fragmentation: after $108$ generations no
organism has the possibility to switch to the next stage of the life cycle.
With $1$ split: no such individual is left after $156$ generations (values
obtained starting with $108$ individuals: $96$ and $146$ respectively; with $%
10^{12}$ individuals: $119$ and $169$). See the table below:

\bigskip

\begin{equation*}
\begin{tabular}{|c|ccc|}
\hline
&  & Population size &  \\ \hline
$k$ & $10^{8}$ & \multicolumn{1}{|c}{$10^{10}$} & \multicolumn{1}{|c|}{$%
10^{12}$} \\ \hline
1 (no split) & $96$ & $108$ & $119$ \\
2 (one split) & $146$ & $156$ & $169$ \\ \hline
\end{tabular}%
\end{equation*}

\bigskip
\bigskip

The allowed generation time under periods of partially relaxed selective
pressure has been extended by more than $40\%$. We next altered the 'simple'
simulation by using a fixed deletion size for gene fragments (instead of
allowing the full range of sizes of our previous model) while at the same
time making it more informative by varying the amount of gene fragmentation (%
$k$, see Figure 1C). This (full Markov chain) simulation gave a surprising
result: we obtained a direct, \emph{quasi-linear} correlation between the
degree of dispersal/fragmentation ($k$) and the number of generations after
which a complete population loses ecological competence (see Figure 2). In
other words, the function $t_{\max }(k)$ has an almost linear growth rate.

\section{Comparing the model to values observed}

A triatomine bug (the insect vector used by \emph{Trypanosoma cruzi}) will
contain on average $105$ parasites (Kollien and Schaub, 1998); while a bite
with the tsetse fly will transmit $0-40,000$ (mean, $\sim \emph{3},000$)
\emph{T. brucei} infective metacyclics, giving rise to over $10^{10}$
parasites in a Vervet monkey (\emph{Chlorocebus pygerythrus}) at peak
infection ((Thuita \emph{et al.}, 2008) and references therein). We chose
populations of $10^{8}$, $10^{10}$ and $10^{12}$ parasites in our simple
modelling (see the table) and $10^{10}$ for our full Markov chain approach.
Under the parameter values of our model the size of the population does not
seem of major importance (cf the table).

We chose to model $\sim 65\%$ of the total mt DNA in use, and $\sim 35\%$
free from functional selection (for trypanosomatids in mammalian hosts, much
less mt DNA seems essential). Doubling time in mammals can be $4.5$ hours if
unimpeded by immune responses, while infections can be sustained for months.
Thus, our model should look at effects over many generations. In our basic
approach we compare only two instances ('no split' and 'only one split'),
but allow all possible deletion sizes. In the full Markov chain approach we
use an approximation (depending on $d$) allowing us to model random
deletions of varying sizes.

Because we have modelled optimal dissemination (always strictly
intermingling constitutively used and conditionally used regions), we
tentatively infer that the fragmentation values of Figure 2 are relevant
below $150$. Recall that the \emph{complete} \emph{T. brucei} required set
of gRNAs for \emph{all} lifecycles is about $150$ (Hong and Simpson, 2003).

Surprisingly, only $77$ rapid passages (at least $600$ generations) in mice
of tsetse fly infective \emph{T. brucei} already gave rise to a homogeneous
population of parasites that could no longer infect the insect and did not
display mt (oligomycin sensitive) ATPase activity anymore (Hajduk and
Vickerman, 1981), demonstrating rapid takeover by mt DNA mutants. These
mutants in the end indeed \emph{divided more rapidly} than the wildtype in
blood but did not contain large scale mt deletions. This could of course be
explained by the fact that certain parts of the present-day highly dispersed
mt genome are still necessary and thus protect against large scale deletions
under these circumstances, as described by our model of selection for
linkage.

Takeover by deletion mutants could thus have been even faster, as mt
replication times are possibly limiting when doubling every $4.5$ hrs. This
also means (cf. Figure 2) that at a fragmentation level of between $\sim 10$
and $20$, deletion mutants would contribute about as much to depletion of
the wt population as all other inactivating mutations ($>600$ generations).
For \emph{T. brucei} we would expect a fragmentation level of $>20$ and $%
<<150$, which seems realistic. However, all \emph{T. brucei} parameter
values are very difficult to estimate in real life, because:

\begin{enumerate}
\item It is not known precisely how little mt DNA is essential in the
mammalian host,

\item Population sizes and generation time estimates are complicated by
waves of parasitemia, reflecting the immune response to alternating variant
surface glycoprotein expression. Very limited mt function and repeating
population bottlenecks could explain why \emph{T. brucei} has the most
extensive editing observed (see below).
\end{enumerate}

A further conclusion, shedding light on another point of contention
(Landweber, 2007; Speijer, 2008): physical linkage of all mt genes is not
essential for the model to work. This aspect is important when considering
the cox1 gene in the parasitic diplonemid \emph{Diplonema papillatum}
apparently split up in $\sim 250$-bp fragments, located on individual
unlinked DNA circles (Marande \emph{et al.}, 2005) and the existence of
minicircles with only one gRNA in certain trypanosomatids (though often
physically linked).

As the modelling makes clear, it is only reducing the chance of large
deletions giving replication advantages which is essential. Losing a
substantial fraction of -one gRNA encoding- minicircles (a 'large' deletion)
will still compromise viability directly, regardless of \emph{physical}
linkage. Networks themselves did probably evolve to make takeover by
(deletion) mutants even less likely (Borst, 1991).

\section{Evolution of RNA editing: gene fragment dispersal to counter gene
loss?}

We model fragment dissemination as advantageous to parasites (Speijer, 2006;
Speijer, 2007), but others defend `neutral' models (such as 'Constructive
Neutral Evolution' (Covello and Gray, 1993; Gray \emph{et al.}, 2010)); see
(Lukes \emph{et al.}, 2009; Speijer, 2010; Lukes \emph{et al.}, 2011;
Flegontov \emph{et al.}, 2011; Speijer, 2011). Parasitism is not essential
for our model to function: any free-living ancestor in periodically changing
environments adapting its mt function could favour mt gene fragmentation
when encountering strong intraspecific competition.

Quickly changing environments could also be at the basis of the glycosome
(Hannaert \emph{et al.}, 2003; Gualdron-Lopez \emph{et al.}, 2012). Rapidly
changing oxygen levels would lead to either aerobiosis with mt respiration
or anaerobiosis with glycosomal activity for energy generation in
non-parasitic ancestors. For these ancestors, the lack of selective pressure
on the mt genome has been invoked as the raison d'etre for RNA editing
before (Cavalier-Smith, 1997).

The fact that novel kinetoplastids and diplonemids have been found in anoxic
deep-sea basins recently, make such an ancestor more likely (Lara \emph{et
al.}, 2009; Edgcomb \emph{et al.}, 2011). Indeed, despite corrections to the
kinetoplastid phylogenetic tree (Moreira \emph{et al.}, 2004; Katz \emph{et
al.}, 2012), \emph{pan-editing is still seen as ancestral}, originating 500
to 700 million years ago (Lukes \emph{et al.}, 1994; Fernandes \emph{et al.}%
, 1993). Kinetoplastids are part of the Euglenozoa (with diplonemids and
euglenoids). Interestingly, another form of mt gene fragmentation occurs in
the \emph{diplonemid D. papillatum} (Marande \emph{et al.}, 2005; Marande
and Burger, 2007), as described above.

Transcript editing has become less extensive in more recently evolved
species (\emph{e.g.} compare cox3 editing in \emph{L. tarentolae} and \emph{%
C. fasciculata} with its pan-editing in \emph{T. brucei}). Editing loss
results from reverse transcription of (almost completely) edited RNA
followed by homologous recombination of the cDNA with mt DNA. The 5' and 3'
homology requirements of the cDNA would sometimes result in a need for 5'
editing of the 'new' transcript encoded by the mt DNA upon recombination.
This is actually seen in cox3 and cytb in \emph{L. tarentolae} and \emph{C.
fasciculata} (Landweber, 1992).

Another nice example of recent editing domain length reduction during
cryptogene evolution is found in the ND8 gene of three related insect
trypanosomatids, \emph{again strongly suggestive of correlations between
life cycle complexity and editing extent} (see below; Gerasimov \emph{et al.}%
, 2012). Pan-editing entails a lot of extra complexity and energetic costs,
so losing it is indeed what one would indeed expect, as soon as reduction of
life cycle complexity allows it.

\section{Pan-editing counters gene loss}

In 1993 Covello and Gray (Covello and Gray, 1993) introduced a general model
for the evolution of different RNA editing forms. RNA editing activity is
first acquired by (a combination of) pre-existing enzymes. Mutations at
'editable' nucleotide positions in the genome occur next. Later on, editing
becomes essential after fixation by a chance process in which an altered
form replaces the original without a selective advantage. In the case of
kinetoplastid editing this model could in principle explain the emergence of
a few 'limited' editing instances, but it is very hard to envisage how it
explains the rapid acquisition of multiple instances of pan-editing with
hundreds of gRNAs.

This model does not identify selective pressure(s) responsible for an \emph{%
active} increase of editing potential. This pressure was postulated to
reside in the fact that gene scattering could protect against loss of
temporarily non expressed mt genes during periods of intense intraspecific
competition (Speijer, 2006), by warding off large deletion mutants
outcompeting wildtype kinetoplastids (see above and Figure 1A). Thus,
ecological flexibility is retained, allowing kinetoplastids to occupy highly
diverse (parasitic) niches and undergo extensive speciation (cf. the
repeated evolution of \emph{T. brucei} strains that can not infect their
insect vectors anymore (Lun \emph{et al.}, 2010)).

Pan-editing seems an effective way of making large deletions improbable. It
necessitates presence of \emph{all} cognate small gRNAs (containing
information for the U sequence) to express a cryptogene (encoding the 'GAC
sequence'). These `gene fragments' colonize mt DNA so that every random
large scale deletion will now contain some information still under selective
pressure (Speijer, 2006; Speijer, 2007). The following observations fit our
model.

\begin{enumerate}
\item Especially transcripts encoding gene products which are crucial in all
life cycle stages should be prone to pan-editing (such as ATP6 and RPS12,
encoding components of the \emph{T. brucei} mt F1Fo ATPase and ribosome
respectively). It is exactly because they are always required that
fragmenting their genes with concomitant spread of cognate gRNA genes in mt
DNA leads to efficient 'mt DNA integrity checkpoints'.

\item The role of gRNAs as such mt DNA integrity checkpoints can also be
performed by evenly distributed tRNA genes. Compare \emph{e.g.} tRNA genes
in human mt DNA (Anderson \emph{et al.}, 1981) and \emph{E. coli}'s multiple
rRNA operons distributed in such fashion (Blattner \emph{et al.}, 1997).
However, exceptionally, all trypanosomatid mt tRNA genes are absent, mt
tRNAs coming from the cytoplasm (Hancock and Hajduk, 1990).

\item There seems to be a clear correlation between life cycle complexity
and amount of pan-editing, \emph{i.e.} gene fragment dissemination (see
above and \emph{e.g.} (Gerasimov \emph{et al.}, 2012). \emph{T. brucei},
with its highly complex life cycle, still encodes a COX III cryptogene, much
reduced in length (see above), so its loss in the mammalian host would be
less advantageous. Much more importantly, its dispersed gRNA segments
protect parts (encoding NADH dehydrogenase subunits) of the genome not used
in the fly (Speijer, 2006; Speijer, 2007). Thus, gene fragmentation implicit
in mt RNA editing can be seen as an instance of a general tendency to
'evenly' distribute genetic information over the genome. In this way large
scale deletion of pieces released from selection becomes impossible, no
longer threatening ecological flexibility.

\item Our modelling demonstrates that the number of generations a population
retains individuals with a full mt genome complement increases
quasi-linearly with increases in RNA editing (\emph{i.e.} with its
concomittant gene scrambling). It thus allows \emph{positive selection of
small, incremental increases of editing potential} under the appropriate
circumstances (\emph{e.g.} less frequent environment-host exchange).
\end{enumerate}

Constructive Neutral Evolution, ignoring 'evolvability' aspects tries to
explain pan-editing as resulting from genetic drift and population
bottlenecks only, but our model explains both its rapid spread and current
patterns. These patterns correlate not only with a history of population
bottlenecks, but \emph{specifically with alternating selection pressures on
mt DNA} (Speijer, 2006; Gray \emph{et al.}, 2010; Lukes \emph{et al.}, 2011;
Flegontov \emph{et al.}, 2011; Speijer, 2011).

\section{Kinetoplastid mitochondrial genome evolution}

Every time kinetoplastids readjust to different environments, founder
effects can occur. The observed rapid speciation and development of
unexpected, 'weird' biochemical properties could then be a natural result.
Such founder effects might explain acquisition of (limited) RNA mt
editing/gene fragmentation in Euglenozoa (Walker, 2007) in the first place.
The kind of selection described here would then quickly give pan-editing in
lineages with alternating mitochondrial demands.

Present-day kinetoplastids show an enormous diversity in mt DNA: maxicircles
of up to 200 kb, gRNA encoding minicircles of 1-3 kb with only one to
several gRNA genes, and 200 kb circles possibly encoding hundreds (Simpson
\emph{et al.}, 2000; Maslov and Simpson, 1994; Blom \emph{et al.}, 2000).
Minicircles clearly evolved later, possibly in response to coding capacity
demands of pan-editing. This increased mt DNA size and led to networks
combating minicircle loss (Borst, 1991). As mentioned, over time editing
becomes less extensive in some kinetoplastids (Landweber, 1992;Gerasimov
\emph{et al.}, 2012), reflecting the complexity and overall burden of the
elaborate editing system.

\section{Conclusion}

We have modelled mt gene scrambling, matching it with observations regarding
the frequency of pan-editing and the ecology and population sizes of
trypanosomatids in which it occurs. Important modelling parameters include:
the population size, the deletion size (as a percentage of total genome
size), the probability of a deletion in each new generation and the
replication advantage value for a deletion. The outcome for every simulation
was the number of generations after which no viable individual was left. We
took the \emph{T. brucei} lifecycle to deduce fixed values, only varying
population size and amount of gene fragmentation. Population size did not
seem a major determinant. Our output correlated extent of
fragmentation/dispersal and increase in number of generations retaining
ecological competence, which turned out to be a quasi-linear function. The
'fragmentation value' of \emph{T. brucei} mt DNA seems to be in the range
suppressing large scale deletion population takeover. To prevent such
takeover, condemning parasites to die with their host, a gene architecture
of tremendous complexity evolved. It links survival of genes released from
selection pressure to survival of genes still under such control. Thus,
while fragmentation reduces efficiency it seems to enhance long term
ecological success. We are currently looking at gene fragmentation of \emph{%
Plasmodium falciparum} (Feagin \emph{et al.}, 2012) and ciliates (Nowacki
\emph{et al.}, 2011) in light of this mechanism, with only the malaria
parasite possibly fitting the mould.

\bigskip

\emph{Acknowledgements. }We thank Nick Dekker for his expertise in
preparing Figure 1, Ronald de Wolf for helpful discussions and Andrew Moore
for valuable comments and suggestions to improve the manuscript. This
research has been partially funded by a Dutch BSIK/BRICKS grant and by Vici
grant 639-023-302 from the Netherlands Organization for Scientific Research
(NWO). SS is supported by the Royal Society. Part of this work was performed
while SS was visiting the CWI; its financial support is gratefully
acknowledged.

