\section{Introduction}

   Cyg X-3 is an X-ray binary that was first discovered in X-rays \citep{Gia67}.
   A short orbital period of 4.8 hr \citep{Par72,Bec73} and
a massive Wolf$-$Rayet donor star \citep{van96,Koch02} make Cyg X-3
a tight high-mass X-ray binary system.
   The nature of the compact object is
still a matter of debate, but there has been growing evidence
of a black hole \citep[e.g.,][and references therein]{Shrader10}.
   Cyg X-3 is the brightest X-ray binary in the radio and
has displayed frequent flaring activity
since the discovery of a giant flare in
1972 \citep{Gre72,Waltman94,Waltman95,Waltman96}.
   \citet{Waltman94} identified three types of variability in
the 1988$-$1992 light curves at 2.25 and 8.3 GHz:
frequent minor flaring activity of $<$1 Jy,
intermediate flares ($>$1 Jy), and major flares ($>$10 Jy).
   In addition, there are pre-flare quenched states of tens of mJy
prior to the flares, and flickering activity during quiescence.
   In very long baseline interferometry (VLBI) observations,
clearly resolved jet structures on tens of milliarcsecond (mas) scales
have been observed  \citep{Newell98,Mio01,Mil04,Tudose07,Tudose10}.
   Even in the observations with unresolved structure, model
fits indicate expanding jet-like structures
\citep[e.g.,][]{Gel83,spen86,Mol88,Sch95,Sch98}.
    However, the VLBI observations have mostly been carried out during major flares,
in particular after the peak, and rarely performed in the rise to the peak
mainly due to the slow reaction time of the VLBI facilities
after the triggering of flares.
    Therefore, it is important to probe the rise phase to understand
the evolution of flares and associated initiation of jets.

   Cyg X-3 exhibits the canonical X-ray states such as low/hard state,
very high state, intermediate state, high/soft state, ultrasoft state,
and quiescence.
   However, the X-ray spectra in Cyg X-3 are complex and significantly
different from other X-ray
binaries \citep[e.g.,][]{MR00,Fender04,S08,Hja09,kol10}.
   In Cyg X-3, the X-ray hardness intensity diagram (HID),
a plot of X-ray flux versus X-ray color (hardness),
does not show a typical hysteresis unlike other black hole X-ray
binaries \citep[][and references therein]{kol10}.
   \citet{Fender04} predict that a relativistic ejection is expected
to occur in a transition to a harder X-ray state across the ``jet line",
while soft states do not seem to produce jets.
   Based on further analysis of the HID, \citet{kol10} conclude that
Cyg X-3 crosses the jet line and a major radio flare occurs
when there is a transition from a ``hypersoft state"
(or canonical state of the ultrasoft state) to a flaring, harder X-ray state.

   On 2007 June 1, Cyg X-3 was found in a radio flare at 1$-$30 GHz with
the RATAN-600 radio telescope \citep{Tru07}.
   The radio flare reached $\sim$1.6 Jy at 1 GHz and
3.9 Jy at 30 GHz the next day, and declined to $<$1 Jy at 1$-$30 GHz
within the next two to three days (Figure 1).
   A few days prior to the radio observation,
the hard X-rays (15$-$50 keV) rapidly increased with the hardening
in the soft X-rays (see the third and bottom panels in Figure 1).
   After a soft state observed since 2007 April \citep{Krimm07}, an extreme soft,
or ultrasoft state was observed with $INTEGRAL$ during
May 21$-$26 (JD 2,454,242$-$2,454,247; \citep{Soldi07,Bec07}.
   During May 27$-$29 (JD 2,454,248$-$2,454,250),
the onset of an X-ray state transition was observed with $INTEGRAL$
and six days later during
June 5$-$12, it reached a hard state \citep{Bec07}.
   Our VLBI observation was carried out
after $INTEGRAL$ reported Cyg X-3 in its ultrasoft state at 3$-$20 keV
\citep{Soldi07}.
   In this paper, we present a VLBI observation of
the 2007 May-June flare of Cyg X-3 in the rise phase
during a state transition from ultrasoft to hard state.

\section{VERA Observation, Data Analysis and Results}

   Cyg X-3 was observed at 22 GHz for $\sim$9 hr,
during ${\sim}$15$-$24 UT on 2007 May 29,
with VLBI Exploration of Radio Astrometry (VERA).
   VERA consists of four 20 m antennas located at
Mizusawa, Iriki, Ogasawara and Ishigaki-jima.
    The baseline length of VERA ranges
from 1019 km (Iriki$-$Ishigaki) to 2270 km (Mizusawa$-$Ishigaki),
and the expected angular resolution at 22 GHz is $\sim$1.2 mas,
depending on the projected baseline in the $uv$-plane.

   A left-handed circular polarization (LCP) was received
because the LCP signal is only available in VERA.
   The received signals were sampled with 2$-$bit quantization
and filtered to 16 MHz bandwidth using VERA digital filter bank
(Iguchi et al. 2005).
   Since the recording rate of VERA DIR2000 recorder is 1024 Mbit s$^{-1}$,
we can record $16~{\rm MHz~bandwidth} \times 16~{\rm IFs}$ in Nyquist sampling
($16~{\rm MHz} \times 2~{\rm samples~Hz^{-1}} \times 2~{\rm bit~sample^{-1}} \times 16~{\rm IFs}
= 1024~{\rm Mbit~sec^{-1}}$).
   Of 16 IFs, 15 were assigned to Cyg X-3
and one IF was assigned to another target,
W75N \citep{Kim2013},
separated by $2.^{\circ}03$ from Cyg X-3.
   All VERA antennas have a facility of dual-beam receiving system
\citep{Kawaguchi00} for simultaneously observing a target and
reference source within $2.^{\circ}2$.
Since a short-duration flare of Cyg X-3 which we present in this paper was
readily detected in the fringe search (see below),
we did not carry out the phase-referencing VLBI technique.
   The data correlation processing was carried out on the
Mitaka FX Correlator \citep{Chi91} at the Mitaka campus of the
National Astronomical Observatory of Japan (NAOJ).

   The VLBI data analysis was performed using
the Astronomical Image Processing System
\citep[AIPS;][]{Greisen2003}.
   For amplitude calibration,
the flux density scale was set using system temperatures measured with
the chopper-wheel method \citep{Ulich76} and aperture
efficiency\footnote{http://veraserver.mtk.nao.ac.jp/restricted/CFP2010/status10.pdf}
depending on the elevation.
    In our VERA observation, we measured the power of a hot load
and a normal sky, corresponding to the ambient and sky temperature.
For the bandpass calibration,
we alternately observed BL Lac for 10 minutes and Cyg X-3
for 70 minutes in a cycle of 80 minutes,
which corresponds to the recording time of a tape of VERA's recorder.
   For the intrinsic variability of BL Lac, the intraday variability of
amplitude is a few percent \citep[e.g.,][]{Quirrenbach92}.
   Therefore, the calibration is free from the intrinsic variability of BL Lac.
   Fringe-fitting was carried out
separately for Cyg X-3 and the calibrator source, BL Lac.

   We detected fringes of Cyg X-3 with a threshold of the
signal-to-noise ratio of 4 in the AIPS task FRING.
   The visibility amplitude from the detected fringes for
each baseline is shown by the black crosses in Figure 2.
   All VERA stations had no systematic and technical problems
during our observation.
   In the Ishigaki and Ogasawara stations,
the system temperature highly varied up to
$\sim$1,500 and $\sim$3,000 K due to bad weather before 20:00 UT,
but became stabilized down to $\sim$500 and $\sim$200 K
after $\sim$20:00$-$21:00 UT, respectively.
   In the Iriki station, $\sim$200$-$250 K was maintained
except the rapid variation before 16:00 UT and $\sim$17:30$-$20:00 UT.
   The system temperature in the Mizusawa station remained fairly
constant, $\sim$150$-$200 K, during the observation.
   We excluded bad data from the data analysis for
system temperatures higher than $\sim$800 K or whenever the
system temperature rapidly changed by 50\% or more within tens of minutes.
   Therefore, Cyg X-3 was detected mainly at the Mizusawa$-$Iriki baseline
before 21:00 UT, while no detection was made at all baselines during
$\sim$17:30$-$20:00 UT.
   There is a short-duration flare (hereafter, “mini-flare”)
for 21:00$-$23:37 UT.
   The system temperature in all stations
became stabilized after 21:00 UT.
   Therefore, the observed mini-flare is independent of systematic effects
such as the antenna system temperature and source elevation
that began to decrease after 19:00 UT.
   Note that the mini-flare is detected on all baselines,
indicating that the mini-flare is
an intrinsic phenomenon of the target source.
   The peak flux density is greater than 1 Jy on
shorter baselines such as
the Iriki$-$Ogasawara and Iriki$-$Ishigaki baselines,
while it is less than 1 Jy on longer baselines
such as the Ogasawara$-$Ishigaki and Mizusawa$-$Ishigaki baselines.

    From the variation of the visibility amplitudes presented in Figure 2,
we infer that the detected fringes are not from a point source.
   The change of visibility amplitudes during the mini-flare
in Figure 2 is intrinsic to the source Cyg X-3,
and reminiscent of, for example,
few hour-scale small flares reported by \citet{Newell98}.
    For the time-variation in the visibility amplitude during the flares,
they suggested that
structure variation was seen.
    They indeed produced a sequence of images from the 20-minute snapshots,
implying a 0.3c jet event
that clearly confirmed the change in the source structure during the flares.
   To investigate whether our observations of Cyg X-3
indicate temporal structural variation, we performed the following two tests.
   First, we tested an image by performing a hybrid mapping method
for all the detected fringes assuming no structural time variation, and
compared the corresponding model visibility amplitudes,
shown by the blue dashed lines
in Figure 2, with the observed visibility amplitudes.
   In the Mizusawa$-$Ishigaki baseline,
the model visibility amplitude over the observing time
(the blue dashed line in Figure 2)
is
larger than the observation
(the black crosses), indicating
that the model is too compact compared to the real size of Cyg X-3.
   Therefore, such a simple model without time variation
does not seem to well represent the observation.

   The second test we performed was an elliptical Gaussian fit
to the observed visibility amplitude.
   Due to the limitation in $uv$ coverage with
only six baselines,
we cannot directly produce a hybrid mapping image for a time-scale short
enough to trace a change of source structure along with the variation of
flux density.
   We checked for phase closure during observations
of Cyg X-3 and obtained a value near zero within ${\pm}10^{\circ}$,
which justifies adopting a single Gaussian fitting
to the visibilities.
   We fitted the visibility amplitude directly
in the $uv$-plane without making images.
   The model fitting was carried out using
Caltech's Difmap program.
   We adopted a single elliptical Gaussian component with a time bin of
10 minutes during the mini-flare, starting from 21:14 UT.
We obtained a flux density, major axis, minor axis, and
position angle of
the
major axis for each time bin from the model fit.
   The orange solid lines in Figure 2 are the visibilities
produced based on the model-fit results.
   In Figure 3, we present the
model-fit results showing
total flux density,
intrinsic source sizes of major and minor axes and position angle.
   The model flux density for each baseline,
shown by the orange solid lines for each time bin,
is in good agreement with the observation (Figure 2).
   Error analysis for the total flux density, the sizes of major and minor
axes and the position angle were estimated using Difwrap \citep{Lovell00}.
   The reduced ${\chi}^2$ for the model fits with time bins of
10 and 20 minutes are ${\lesssim}0.82$.
   The model fit with a time bin of 30 minutes
has a reduced ${\chi}^2 {\lesssim} 1.1$,
approximately equivalent to $1 {\sigma}$ \citep[e.g.][]{Dodson03}.
   The time-evolving features of axes and position
angles for the time bins of 10, 20 and 30 minutes
are consistent with each other.
   In the following discussion,
we exclude the first data point in Figure 3
because it is a marginal detection,
and the last data point
because the time bin is very short ($<8$ minutes).
   The major axis kept increasing since 22.5 UT,
${\sim}1$ hr before the peak flux density.
   In the two mini-flares observed by \citet{Newell98},
the extended structures (jet expansion) appeared
${\sim}1$ hr prior to the peak of each flare,
proportional to the increase in flux density
\citep[see Figures 4 and 5 in][]{Newell98}.
   Similarly,
the increase of the major axis size in Figure 3 may indicate
a structural change, or jet ejection, during the mini-flare
due to synchrotron emission.
   The position angle of the major axis is relatively stable
with time,
particularly after ${\sim}21.75$ UT (the bottom panel in Figure 3).
   This position angle of ${\sim}50{\pm}10^{\circ}$ is
in agreement with the observations of \citet{Newell98} and \citet{Sch95}.
   In the following section,
we discuss the possible causes of the structural change.

\section{Discussion and Conclusion}

   In this paper,
we present a short-lived ($\lesssim$3 hr) flaring event, or mini-flare
during an X-ray state transition from ultrasoft to hard state
in the early phase of the 2007 May$-$June flare.
   No radio observation was reported prior to the peak of a few Jy
except our VLBI observation.
   To make matters worse, no hard X-ray observation was made
during most of the 2007 May-June flare.
   However, it is clear that
a few days before our observation
the hard X-ray flux of $Swift$/BAT rapidly increased
after its $\sim$12.5 day low flux quenched state  (Figure 1).
   At the same time, the soft X-ray flux ($RXTE$/ASM)
decreased and a few days later
the radio flux was observed near the peak.
   It is well known that the strong correlation of
the radio and soft X-ray emission disappears immediately before
the flaring periods \citep{watanabe94,Chu02}.
   In flaring states, there is a correlation between
the radio and hard X-ray emission,
while the soft X-ray is anti-correlated with the two wavebands \citep{Mc99}.
   Therefore, at the time of our observation, the 2007 May$-$June flare was
probably already in the early rise.
   To confirm this,
we checked the multiwavelength light curves
of a few previous flares.
   In most of the major ($>$10 Jy) and intermediate ($<$10 Jy) flares,
the radio fluxes at GHz wavebands increase rapidly
from a few hundred mJy to a few Jy
when the $Swift$/BAT flux increases
from $<$0.005 to $\sim$0.02 counts s$^{-1}$ cm$^2$
\citep[$Swift$/BAT Web site; RATAN Web site\footnote{http://www.sao.ru/cats/$\sim$satr/XB/CygX-3/};
e.g.,][]{Tru08}.
   Therefore, at the time of our observation,
the 2007 May$-$June flare was most likely in its early rise phase
since the $Swift$/BAT flux was $\sim$0.02 counts s$^{-1}$ cm$^2$ (Figure 1).
   In Cyg X-3, small-scale flares on a time-scale of a few hours or less
have been observed during low-level activity of $<$600 mJy at 5$-$43 GHz
\citep[e.g.,][]{Newell98,Marti00,Ogley01,Mil09}, and
associated jet ejections were confirmed in some cases
\citep{Newell98,Marti00}.
   However, a flare on a similar time scale prior to
major or intermediate flares, such as the mini-flare of $\sim$1 Jy
that we have presented in this paper, has not been reported.

   Due to the absence of published X-ray spectra
during the radio flare of 2007 May-June
($\sim$JD 2,454,250$-$2,454,256; see Figure 1),
it is not clear whether Cyg X-3 underwent an intermediate
or very high state prior to the observed hard state
\citep{Bec07}.
   \citet{kol10} find that there are two types of flares:
flares with hard X-rays and without hard X-rays,
and the latter does not reveal
any jet-like structure \citep[also see][]{Fender04}.
   They point out that the first two $INTEGRAL$ observations
in 2007 May (JD 2,454,241 and JD 2,454,246) during a quenched state
correspond to the hypersoft or ultrasoft state (Figure 1).
   Therefore, a further transition to a harder state in
2007 May$-$June would result in crossing the jet line in HID.
   From Figure 1, it is clear that Cyg X-3 already changed
its state from the ultrasoft to a harder state
with rapidly increasing hard X-rays of $Swift$/BAT
during our observation.
   Therefore, the structural change in our observation (Figures 2 and 3)
indicates that the mini-flare is plausibly produced
by a jet ejection in the early phase of the 2007 May$-$June flare.
   Furthermore, \citet{Tru07} point out that,
 prior to the major flare in 2007 May$-$June,
the observations of $Swift$/BAT in mid- to late 2007 April
\citep{Krimm07}
and $INTEGRAL$ in 2007 May 21
\citep{Soldi07}
clearly indicate a  possible jet ejection event.

 What can be a cause of the mini-flare?
   The most appropriate model to account for such a mini-flare
would be the so-called shock-in-jet model \citep{MG85}.
   The model is an analytic, parameterized shock model describing
the synchrotron emission associated with relativistic jets
to account for the radio and near-infrared light curves
of flaring activities in quasars and microquasars
\citep{Turler00,Turler11}.
   The model can reproduce the variable jet emission in microquasars,
very similar to that of extra-galactic jets in quasars,
suggesting that the physical characteristics of the relativistic jets
are independent of the black hole masses
\citep{Turler07}.
   As for small-scale flares in other microquasars, for example,
\citet{Turler04} reproduced a series of flaring events of
a few hours observed at 2.2 ${\rm {\mu}m}$ and 2$-$6 cm
during a plateau state in the microquasar GRS 1915+105.
   In the case of Cyg X-3,
the model accounts for a series of small-scale, tens of minutes to an hour,
flaring activities in the decay phase of a half-day long radio flare
(${\lesssim}500$ mJy) at 1.4$-$43.3 GHz observed in 2002 January 25
\citep{Mil09}.
    The model can also reproduce the major and intermediate flares of Cyg X-3
by adopting a different time-scale, peak flux density
and frequency of each flaring activity
at 2.2 ${\rm {\mu}m}$ and 2.25$-$15 GHz
\citep[e.g.,][]{Turler11},
suggesting that the same phenomena are at work.
    The difference is that the brighter flares peak at lower frequencies
and have longer time-scales than fainter flares, which is
consistent with the formation of shocks further downstream in the jet
\citep[for further discussion, see][]{Mil09,Turler11}.

   In our observation, the model fits with
a single elliptical Gaussian component
are in good agreement with the observed visibility amplitudes.
   The model fits indicate the variation of axis size
during the mini-flare, similar to
the previous observations of the two mini-flares
at different flux scales, both of which included jet imaging
\citep{Newell98}.
   A possibility of jet ejection is also supported by
other multiwavelength observations in the radio, soft and hard X-rays
during the 2007 May$-$June flare
\citep{Tru07,Krimm07,Soldi07,Bec07,kol10}.
   In particular, it has been suggested that the jet ejection is highly
expected during a state transition from soft to harder X-ray state
\citep[e.g.,][]{kol10, Fender04}.
   The shock-in-jet model
for microquasars suggests
that various flares with different time scales can be described by
a synchrotron flare
\citep{Turler00,Turler04,Turler11,Mil09}.
     Therefore, together with such observations and models
for flaring activities in Cyg X-3 and other microquasars,
our model fits suggest
a high possibility of synchrotron flares and
an
associated jet ejection
during the X-ray state transition from ultrasoft to hard state
in the 2007 May$-$June flare.
    The observation of jet activity in
this particular state
transition has been a long-awaited issue.
   We present
observational evidence of such an expected
jet ejection event in the state transition for the first time.
The next challenge is to image a variety of jet ejection events with the
combined array of VERA and Korean VLBI
Network\footnote{http://kvn-web.kasi.re.kr/en/en{\_}index.php}
which will provide better UV coverage and
allow us to make snap shot images of
rapidly variable sources.

\acknowledgments
We thank the anonymous referee
for constructive suggestions and a thoughtful review of the manuscript
that has helped to substantially improve the contents of this paper.
   We thank Katsunori Shibata at NAOJ for arranging the dynamic scheduling.
   We thank Jim Lovell and Kotaro Niiuma for providing Difwrap and for useful comments on it.
   We also thank Makoto Inoue and Kiyoaki Wajima for comments on Difwrap.
   S.J.K. acknowledges partial supports from
the Korea Science and Engineering Foundation (R01-2008-000-20002-0),
a WCU Grant (No. R31-10016),
and a grant from the Kyung-Hee University in 2013.


\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.