Model of invasion of a population by transposable elements presenting an asymmetric effect in gametes

Dynamics of population invasion by transposable elements is analyzed and simulated, using a model with a very large number of transposition sites. The properties of the model are determined in the framework of a conflict between transposition capabilities of the elements and their harmful effects on the host genome. Equations are developed for the mean and the variance of the number of elements at equilibrium. We use simulations


INTRODUCTION
About 15% of the eukaryote genome consists of a family of repeated and dispersed DNA sequences. Many of these sequences have been described before, and some of them have been found capable of mobility (review in Berg and Howe, 1989).
Several models have been proposed to characterize the distribution laws of these transposable elements in populations as a function of different variables such as their transposition and excision rate, and the selective values given to carrier individuals (reviewed in Charlesworth, 1985, and Brookfield, 1986.
Generally speaking, in all sexed organisms, these models have shown that a family of elements could be kept in stable equilibrium by the opposed effects of replicative transposition and selection against the harmful carriers.
However, much experimental research has proved the existence of self-regulation mechanisms by which the probability of transposition of an element decreases as a function of the number of elements of the same family present in the host genome (reviewed in Berg and Howe, 1989). Different models have shown that such selfregulation could also lead to a state of stable equilibrium for the distribution law of a given family of elements (Charlesworth and Charlesworth, 1983: Langley et al, 1983, Charlesworth and Langley, 1986Langley et al, 1988;Rio, 1990).
In Drosophila melanogaster, the research on hybrid dysgenesis induced by families of I, P and hobo elements (reviewed in Berg and Howe, 1989;and in Berg and Spradling, 1991) has generated a set of data by which more specific basic models can be conceived. Such models have been proposed for describing the evolution of such systems, in consideration of some of their characteristics, but either by dealing only with the case of a single transposition site (Ginzburg et al, 1984;Uyenoyama, 1985) or with an infinite number of sites, and analyzing the selective values at the individual level (Brookfield, 1991).
In the present article, we analyze a model for a family of transposable elements whose transposition and excision rates are functions of the copy number, and whose dysgenic effects depend on the type of crossing. The invasion conditions of these elements in a population are determined analytically. When a state of internal equilibrium exists, the mean and the variance of the distributions are found.
Simulations are used to verify the equations, and as a basis for discussing invasion rates of this type of element in populations.

DESCRIPTION OF THE MODEL
The model developed here is based on the opposing actions of a transposase and of a repressor, whose reciprocal concentrations, and thereby the effets, depend on the element copy number. The equilibrium or disequilibrium existing between these 2 components depends on the direction of crossing, and is established in the zygote. Adults have a selective value linked to the possible dysgenetic effects affecting the zygote they come from.
In the model the number of sites (T) is supposed large enough to allow any transposition into an empty site. The gametes are characterized by the number (ranging from 0 to T) of active elements they contain. For the ovum, this number is taken as an index of the concentration of repressor, as we assume the rate of transposition is simply controlled by the repressor present in this gamete. Considering T as very large, the frequency of occupied sites does not appear as a pertinent parameter and we address here only the distribution of copy number per gamete. The element copy number distributions in the spermatozoa and in the ova of generation t are considered to be identical. We define them by (pt(0), pt(1), ... pt(T)), with: The zygote obtained by crossing an ovum containing i elements and a spermatozoon containing j, is denoted (i, j). We must distinguish between 3 types of crossing (table I).

First case
The spermatozoa contain no elements, so that j = 0. We suppose that the transposable elements have no effect on the (i, 0) zygote because the equilibrium between the repressor concentration and the number of elements in the egg is not disturbed.
The selective value of these zygotes is taken as reference, and is therefore set equal to unity, w(i, 0) = 1. Finally, we let (1 -G(i)) be the frequency of the gametes without elements in the set of gametes produced by the (i, 0) type zygotes.
When the sites of the elements are on a single pair of chromosomes, and when there is no recombination possible, we have G(i) = 1/2 for i > 0.

Second case
The ovum has no repressor (i = 0) and the spermatozoon has elements (j > 0). In this configuration there is a high level of element activity, and we let A(j) be the mean increase in the number of elements for the type (0, j). A(j) is the mean number of elements created, less the mean number of elements lost. W(j) = 1-S(j) represents the selective value of these zygotes. [1 -B(j)] is the frequency of gametes without elements in the set of gametes produced from type (O J ) zygotes.
Third case i > 0 and j > 0.
We define a(j), w(j) = 1 &mdash; s(j) and b(i,j). The values a and w are supposed to depend only on j, which induces the disequilibrium between the number of elements and the repressor concentration. As in the first case, the repressor concentration of the ovum is assumed to balance its element copy number. In the zygote, the disequilibrium therefore depends only on the number of elements j introduced by the spermatozoon. The values a(j) and w(j) = 1 &mdash; s(j) correspond to the mean increase in the number of elements and to the selective value of the zygotes (i,j). In these zygotes, the presence of repressor limits the activity of the j elements introduced, which means a(j) < A(j) and w(j) > W(j). Finally, and as before, we define (1 &mdash; b(i,j)) to be the frequency of gametes without elements resulting from these zygotes. Table II  In order for the element frequency to be able to increase, the function pt+1 (!) -T p t(O) = pt+1 (0) -(1 &mdash; L P t(i)) must be strictly negative in a neighborhood of i=i 1 (Pt(1) = 0,... ,pt(T) = 0) .
As p t+i (0) is differentiable in (0, ... , 0), we have to compute the partial derivatives in (0, ... , 0). We get for 1 ! k $ T (Appendix 1): A sufficient condition for an increase of elements starting from a small initial number of gametes possessing elements is therefore G(k) + W(k)B(k) > 1 for 1 ! k ! T. When the gametes introduced have few elements, it is sufficient for the first inequalities alone to be satisfied..
These inequalities are easy to interpret because G(k) [respectively W(k)B(k)] is the probability that a (k, 0) type zygote [resp (0, k)] that would be viable and nonsterile will produce a gamete containing at least one element.
It will be seen that this model generalizes that of Ginzburg et al (1984), which assumes a single insertion site, or that the number of elements has no effect, and on a single pair of chromosomes. In this case, the notation is G(i) = 1/2 for 1 ! i ! T; S(j) = S; s(j) = 0; B(j) = !3+ 1/2(1 -/3) = 1/2 + /3/2 for 1 ! j ! T; b(i,j) = 1 for 1 ! i ! T and 1 ! j ! T (!3 being the probability that the maternal genome be contaminated by transposable elements in a (0, j ) type mating. The T inequalities are identical with G(k) = 1/2, W(k) = 1 -S and B(k) = 1/2 + /3/2. So we once again find the necessary and sufficient condition of expansion which they reach in their special case, ie fl > S/(1 -S).
In the present model, B(k) is not fixed, but rather depends on the transposition process. In the case of only one chromosomal pair and assuming that the k elements of the paternal chromosome are not excised, that the increase in the number of elements is A(k) for any (0, k) zygote, and that any new element is inserted randomly in 1 of the 2 chromosomes, we get B(k) = 1 -(1/2)!!)+!.
The kth inequality is then written Each inequality yields a relation between S(k) and A(k) for determining the conditions under which the element copy nomber will increase.
The hatched area of figure 1 corresponds to the values of S(k) and A(k) verifying this inequality. The harmful effect of the transposable elements can increase as the increase in the number of elements created itself becomes greater. ' Analysis of the positions of equilibrium Analysis of the mean Pt = ( Pt (0) , ... , ,pt(T)) is an equilibrium point if p t+ i(i) = p t (i) for 0 ! i ! T. The modelling described here cannot be used to determine the p t+1 values as a function of P t for i > 0. This is possible only if we know, for each type of zygote, the distribution of the gametes produced as a function of the number of elements they contain. Each of these distributions requires T parameters (the sum of them being less than or equal to unity) in order to be defined. It should be possible to reduce this excess of parameters by adopting assumptions concerning the mode of action of the transposable element. This problem is not addressed in the present paper. Instead, we attempt to obtain the equations for the mean and the variance of the distributions of elements at equilibrium, when it exists. Such equilibria have been found for the corresponding model of Ginzburg et al (1984).
The mean and variance depend on the parameters previously defined. This way we get Pt+ i (0) = Do!D (see Analysis of initial element propagation conditions) and the mean E(X t+1 ) = E(Y t+1 ) for the variables X t+i (resp Y t+ ,), number of elements in the ova (resp in the spermatozoa), of the (t + 1)th generation (see Appendix 2). At a point of equilibrium we have: When pt(0) ! 1 we can use the variable X', which follows the law of X t conditioned by the gametes containing elements. We get E(X t ) _ (1 -p)E(X'), with p = p t (0). For a point of equilibrium (p, pl, ... , pT) other than (1, 0, ... , 0), equation [1] can be written: In the case considered before of a species having only one pair of chromosomes, and supposing that the elements are not excised, the (i, j) type yields no gametes without elements; so we have b(i, j) = 1 for i > 0 and j > 0, and therefore e = 0, which corresponds to an equilibrium where all of the gametes would possess elements. p = 0 is then a solution of the equation [1]. It is the only equilibrium possible if d > 0, because c > 0, (w(j) > W(j)). This situation occurs in particular when the inequalities related to the element copy number growth conditions are verified.
When there is more than 1 pair of chromosomes, or when the element can be excised, b(i,j) may be other than unity; but it approaches it very quickly as i and j increase. It is therefore not surprising to find populations in equilibrium in which p(0) can be considered zero. If it is, equation [2] is reduced too, and the mean number E(X') of elements per gamete satisfies: If a(.) and w(.) are linear functions (a(j) = a.j.w(j) = i -(s.j), this equation is written: in which Var(X') designates the variance in the number of elements per gamete. E(X') therefore does not depend only on the mean increase and the selective value, but through the variance of X' it also depends on the dispersion of the insertion-excision process. This variance therefore deserves being analyzed.

Analysis of the variance
Let Var(i, j) be the variance of the number of elements of gametes produced by type (i, j) zygotes. Even with a deterministic model of the number of transpositions and excisions in these zygotes, Var(i, j) is not zero.
Var(i, j) is analyzed in Appendix 3 for the case of a single pair of chromosomes. These new parameters are introduced in order to calculate the variances Var(X t+d = Var(Y t+i ) of the number of elements in the gametes of generation (t + 1).
At a point of equilibrium we have Var(X t+1 ) = Var(X t ), which provides a third condition [3] of equilibrium (see Appendix !,). This condition depends on the third moment of X'. Even when p(0) = 0 and the functions a(.), w(.) and Var(.,.) are simple, the simulations (see below) have shown the importance of the third moment, which has in no case been found to be close to zero. Equations [2] and [3] cannot therefore be used to find E(X') and Var(X').
As the invasion dynamics of the elements are just as interesting as their mean and variance at equilibrium, we chose to simulate the process rather than simplify the equations by approximation. However, the mean increase in the number element per type of zygote is not sufficient and we must take into account the way a (i, j) zygote produces new elements.

SIMULATION AND NUMERICAL ANALYSES Evolution simulation program
To reduce the simulation program run time and have a first approach to the process, the program computes the case of a single pair of chromosomes without recombination. This has its effect on the numerical results by way of G(.), B(.), b(.,.) and Var(., .), but does not change the mean value. Moreover it will allow an introduction to the general features of the phenomena.
The user has to define the element copy number distribution in the gametes of the original generation, as well as the functions A(.), a(.), W(.) and w(.). Table III summarizes the list of parameters used in these simulations.
The functions allowed are of the form: The mean increases are therefore the result of a (U, u) transpose and a (V, v) excision process (Charlesworth and Charlesworth, 1983).
To obtain the p t (n) frequencies at the tth generation, we have to determine for each (i,j) type zygote the gametes it will produce. However, the knowledge of A(j) and a(j) are not sufficient, and the distribution of the transposed and excised elements around these means is necessary. The program allows the user to choose between a distribution ranging between the two integers to either side of the mean, or a Poisson distribution. From such a distribution the final composition of gametic types is determined, giving to each chromosome produced its number of new elements.
The simulation stops at the tth generation when the frequencies p t (n) and p t -i (n) are within 10-6 of each other (0 fi n fi T).
The stability of the mean and the variance of the tth generation is verified by computing the mean and the variance of the (t+1)th generation from the formulae that led to equations [2] and [3].

Examinations of a few special cases
The examples considered here are based on linear functions (D = E = d = e = F = f = 1), and the increase of elements is distributed over 2 consecutive integers.
Using the notation of the model, we get for the average increases in the number of elements: and for the selective values of the zygotes: Using the available experimental data (from Bingham et al, 1982;Engels, 1988; Berg and Spradling, 1991) for the P-M system of Drosophila melanogaster, orders of magnitude were defined along with rates of insertion, excision and selective values.
A first series of simulations, carried out as a check, shows as expected that, when there is no deleterious effect (S = s = 0) the mean of the number of elements increases indefinitely, at a rate that depends on the insertion and excision rates. In a second series of simulations, the relations between the harmful effects of the elements and their regulation capacities were examined.

Variations with counterselection and self-regulation of elements
When the mobility of the elements causes harmful side effects, the variations depend on the ratios between the various parameters. For a mean increase of the order of 0.25 per element in dysgenic mating ((0, j) zygote), and considering an excision rate 100 times smaller than the insertion rate (Engels, 1988), we get U = 0.252 and V = 0.002. If the self-regulation phenomena did not exist in the (i, j) zygotes, the parameters u, v and s would be equal to U, V and S, respectively. Under these conditions (fig 2), the simulations show that the transposable elements may invade the population rather quickly (250 generations) when the deleterious effect is not too great (S = 0.05), but can only set in once S reaches a threshold value, which is S = 0.11 here. For intermediate values of S(S = 0.08), the invasion time will be greater (500 generations) and only a part of the gametes would have transposable elements.
The self-regulation effect in the (i, j) zygotes can then be analyzed by assigning values 10 times smaller to u and v, or u = U/10 and v = V/10, and choosing a low deleterious effect of s = 0.012. The simulation is initialized with 1 gamete among 1000 carrying a single element, or p o (1) = 10-3 . The curves obtained for different values of S (0.05 and 0.08) are given in figure 3. Here, the elements can totally invade the population, while the selection coefficient S against the dysgenic zygotes is < 0.11. The invasion is slower than before and the population reaches a stable equilibrium in 1 500 to 2 000 generations, but with a higher mean number of elements (13.4 on the average, and with a SD of 4.8). On the other hand, the frequency of the gametes without elements rapidly diminishes and becomes practically zero in < 350 generations.
It will be noted that the invasion is not necessarily assisted if the rare founder gametes carry a large number of elements. This is due to the fact that the majority of the zygotes in the first generations will be dysgenic and very highly counterselected.
When the spermatozoa introduced carry more than 1/S elements, the (0, j) zygotes provided will be all dysgenetic (W(j) = sup(0,1 -S.j) = 0) and there will be no invasion. However, the invasion might occur for lesser values. Such a case is illustrated in figure 3 with p o (11) = 10-3 and S = 0.08. We then observe that the evolution proceeds in a manner similar to the previous cases, but a little less rapidly.
Effects of variations in the insertion-excision and self regulation rates As the excision rates under dysgenic conditions are known only approximately (to within a factor of 10; see Engels, 1988), a set of simulations were run to get an idea of its role in the case of our model. Using the conditions of the preceding example, of a mean increase of 0.25, we chose to increase the excision rate 10-fold, or: and s = 0.012.
The results obtained (fig 4) show that the evolutions are quite similar to those obtained previously, with equilibrium achieved in 1500 to 2 000 generations with a mean number of elements of 12.6 and a SD of 4.7. In fact, when the insertion rate is much greater (10-100-fold) than the excision rate, variations in the latter have little influence on the evolution dynamics.
In the same way, the role of element mobility regulation and the dysgenic effects observed in (i, j) type zygotes can be evaluated by modifying the values of the corresponding parameters. We examined the case of relatively ineffective regulation in this way, by increasing the values of u and v, and generally the simulations show that the invasion time is shorter.
In the examples given in figure 5 (with U = 0.252, V = 0.002, u = U/2, v = V/2 and s = S/2), the invasion process is possible only for S < 0.11, and the equilibria are achieved at between 250 and 500 generations. In the limiting case, this situation tends toward the one presented before (fig 2), in which self-regulation did not exist in the (i, j ) zygotes.
Finally, a last set of simulations was undertaken, to estimate the effect of very high transposition rates that could lead to a mean increase A of one element. The values chosen (fig 6) of U = 1.5, V = 0.5, u = U/10, v = V/10 and s = 0.03 when the invasion takes place (S E [0, 0.24)) showed that the invasion is faster (from 200 to 800 generations for the 3 examples considered) but leads to equilibrium values very similar to the previous ones (average of 13.7 and a SD of 4.9).

Examination of an invasion in a sequence of stages
The model proposed can also be used to study an invasion occurring in successive waves. An original population A is invaded by transposable elements and then, after a certain number of generations, a part of its individuals emigrate into a fresh population B having none of these elements. Under these conditions, the distri-bution of elements in the gametes introduced into population B will be more or less close to the equilibrium distribution of population A, but usually very different from the initial distribution introduced into this same population A.
The simulations were carried out to create mixed populations from 90% gametes without elements and 10% gametes originating from a parent population in equilibrium. In all cases, the mixed populations evolved toward the equilibrium state of the parent population while the parameters remained unchanged.
However, the analysis of several families of transposable elements has revealed the formation of deleted elements in the course of the generations, which might play a role in the dynamic regulation of the invasion (Black et al, 1987;Jackson et al, 1988;P6riquet et al, 1990;Raymond et al, 1991).
The parameters then have to be modified, as the regulatory process decreases the mobility of the elements and thereby their dysgenic effects. Considering the importance of the insertion rate, a series of simulations was made with only this parameter modified, in order to represent the invasion of a series of 3 populations. The parent population is defined as U = 0.27, V = 0.02, S = 0.05, u = U/10, v = V110 and s = 0.012. ' At equilibrium, 10% of the gametes are introduced into a population with no elements, which then evolves with U = 0.22. The same process will be repeated in a third population for which U = 0.15.
The results given in figure 7 show that each population then reaches a specific state of equilibrium, with lower and lower average numbers of elements: 12.6, 6.5 and then 3.0. It is also seen that these states of equilibrium are achieved more and more quickly. Such a process of successive colonizations can therefore lead to the establishment of a geographical differentiation in the distribution of the number of elements. It should be pointed out, though, that considering the formation of deleted elements, each population will contain a mixture of different elements. The values obtained by our simulations concern only complete and active elements.

DISCUSSION
The model and dynamic simulations of transposable elements presented here are based on a genetic approach to the phenomena of hybrid dysgenesis, described mainly for D melanogaster. The model leads to a generalization of the model of Ginzburg et al (1984) for a large number of transposition sites.
The equations for the mean and the variance of the number of elements at equilibrium depend on the third moment, however, which cannot be neglected. The validity of the model has been confirmed by simulations, mainly by examining the population invasion dynamics. The main results of these simulations can now be discussed and compared with knowledge acquired on P and hobo elements in D rnelanogaster, although the model is an oversimplified version for the complex mechanisms of regulation known for these elements.
The parameter values used in our simulations come from the insertion and excision rates observed in dysgenic matings of the P-M system. Although these values are only orders of magnitude, they do show the impact of the various model parameters. The element invasion conditions is given by a system of inequalities stating the relation between the selective coefficient of the dysgenic zygotes and the net increase in the element copy number A. When it occurs, the invasion is slower with increasing S. With the values used, the critical threshold of S is of the order of 0.10. The conflict between the deleterious effect of the elements and their invasion can be decided either by increasing the transposition rate or by decreasing their deleterious effect, or by acquiring a self-regulation process.
In the case of P and hobo dysgenic systems, the harmful effect due to sterility is smaller at the temperatures usually encountered by the individuals, and a complex self-regulation system was developed. The present model does not account for cytoplasmic type repression effects, but clearly shows the role of chromosomal regulation which, while diminishing the harmful effects when the egg has elements (s < S), favors the survival of the carriers and allows the effective invasion of the population by a large number of elements.
It is particularly interesting to note that, for a species like D melanogaster, the values used lead to a state of equilibrium in 1 500 generations, for a population that had one gamete to begin with out of 1000 carrying a single element. Considering the species, such an invasion process would be of the order of some 100 years, which is compatible with the hypothesis of recent invasion by P and hobo elements (Kidwell, 1983;Anxolab6h!re et al, 1988 ;Pascual and P6riquet, 1991).
So the importance of setting up self-regulation mechanims seems to be primordial in the evolution of these systems. We can then understand that mutant elements capable of participating in the global regulation process can be retained by selection.
This would be the case of the KP elements and other deleterious elements in the P-M system, as well as the Th element in the hobo system. A model has been developed recently that includes the effect of this type of element (Brookfield, 1991). The author shows that deleterious elements can indeed be favored by selection if they are a favored substrate for transposition, and that the populations then consist of a combination of complete and deleted elements. In wild P strains, the total number of element is ! 50. For those analyzed, some one-third are complete, or ! 15 of these elements, which corresponds to the values obtained in our simulations.
Finally, the impact of the deleted elements on the dynamics of the whole also appears to be important in the case of invasions in series. In our simulations, we saw that invasion was not facilitated simply when the gametes of the original population contained more than one element. This suggests that the diffusion of elements starting from individuals that are part of a population that has already been invaded would be based rather on those individuals containing the fewest complete elements.
Moreover, the presence of incomplete elements, limiting the activity of active ones, leads the receiving population to a state of equilibrium at a lower number of complete elements. When the process is repeated from population to population with a low elements transfer flow, it can lead to the onset of a decreasing gradient of the element copy number in all of the populations. This process is thus consistent with the invasion of European strains by American P ones and their dilution into M cytotype, leading to the gradual variation currently observed. Here again, the values obtained in our simulations, from 12.6 to 3.0 complete elements, are compatible with data observed in wild populations: from 35 elements in all for the French populations, 7 or 8 for those of central Asia (P6riquet et al, 1989). It will be noted that, if the element transfer flow continued in the course of generations, all of the populations would then tend to homogenize, which is the case of the I and Hobo elements, for which there is no evidence of geographical differentiation (reviewed in Bregliano and Kidwell, 1983;Pascual and P6riquet, 1991). For the P-M system, for which the invasion in D melanogaster appears to be the most recent, the homogenization process would be under way and recent observations on the French population do in fact show a trend in this direction over the last decade (Fleuriet et al, 1992).
As the model presented here suggests, this process could extend to all European populations, but this would require a relatively long time. We will be able to answer these questions by analyzing wild populations and comparing the results with an extended version of the present model, introducing recombination and segregation between more than one pair of chromosomes and taking into account some aspects of the mechanisms of regulation of these elements. Such a program is under investigation. '