A note on the population genetic consequences of delayed larval development in insects

Observations by Dobzhansky’s group in the 1940s suggesting that the presence of recessive genotypes could account for lower larval developmental rates in Drosophila melanogaster were not confirmed at the time and all subsequent investigations on this subject focused on the analysis of ecological models based on competition among pre-adult individuals. However, a paper published in this journal in 1991 eventually confirmed the finding made by Dobzhansky and his co-workers. In this report, we provide a theoretical analysis of the population genetic effects of a delay in the rate of larval development produced by such a genetic mechanism.


Introduction
Based on an analysis of samples from natural populations of Drosophila pseudoobscura collected from three localities on Mount San Jacinto (California), Dobzhansky et al. (1942) showed that homozygotes for genes located on the second and third chromosomes in this species had a lower developmental rate than expected. These authors suggested that, under controlled culture conditions, the viability and developmental rates were dependent on the demographic density of the larvae. In contrast to the genetic finding, the latter observation was subsequently confirmed by several authors (Bakker, 1959(Bakker, , 1969De Witt, 1960;Robertson, 1963Robertson, , 1964Barker and Podger, 1970;Huang et al., 1971;Tosic and Ayala, 1981;Mather and Caligari, 1981;Ménsua and Moya, 1983;Adell et al., 1988), not only for drosophilid populations but for other insects as well. As a result, all subsequent research on this subject has centered on the analysis of ecological models based on competition among pre-adult individuals. Oliveira and Cordeiro (1981) detected a major gene effect on delayed development in Drosophila melanogaster that was later associated with a recessive gene located on the second chromosome of this fly (Oliveira et al., 1991).
Since the development of molecular genetics in the early 1990s many authors have examined the role of genes involved in the regulation of developmental timing and body size determination in Drosophila. Most of these stud-ies focused on genes related to apoptosis and to the production or release of ecdysone and other ecdysteroids (molting hormones). The results of some relevant articles (all dealing with this hormone in D. melanogaster) are summarized below. Other authors described the importance of additional gene-controlled systems, such as those involved in the regulation of adenosine deaminase-related growth factor ADGF-A, in developmental timing (Dolezal et al., 2005) and in the regulation of the activities of oxidative phosphorylation complexes through nuclear-encoded mitochondrial tyrosyl-tRNA synthetases (Meiklejohn et al., 2013). Oldham et al. (2000) and Quinn et al. (2012) provide appropriate reviews of this subject. Sliter and Gilbert (1992) showed that loss of function mutations at the dre4 gene caused ecdysteroid deficiency and developmental delay or arrest. Later, Bialecki et al. (2002) showed that isoform-specific null mutants for the orphan member E75A gene had a reduced ecdysteroid titer that resulted in developmental delay. McBrayer et al. (2007) showed that prothoracicotropic (PTTH) hormone production was not essential for molting or metamorphosis but that its lack resulted in significantly delayed larval development. Ghosh et al. (2010) showed that the gap gene giant (gt) regulated ecdysone production through specification of PTTH and that mutants for this gene exhibited delayed larval development (in addition to other effects).
Using elements from the theory of finite differences, difference and differential equations, and dynamic systems, we examined the problem based on a model with discrete generations and with the total number of individuals (adults plus larvae) kept constant. We also summarize the results obtained with alternative models involving an increase or decrease in the total population number. The case of continuous generations (to be presented in detail elsewhere) is also briefly discussed.
Description of the main model in which the total number of individuals is kept constant over discrete generations The central principle of this model is that the total number of individuals (adults + larvae) is kept constant in a population with discrete generations. For this, let A and a be a pair of segregating alleles at an autosomal locus. In the homozygous state the recessive gene a induces a delay in larval development (increase in the diapause phase), measured in an integer number of generations: crosses therefore take place randomly (panmixia) only among individuals belonging to the same generation. The population size, kept rigorously constant throughout generations, is assumed to be so large that the effects of random sampling genetic drift are negligible. Consequently, the mathematical treatment is deterministic and gene frequencies are interpreted as probabilities. The effects of differential migrations, selection pressure and mutation are also considered negligible.
In the equations describing the model we will use the following parameters: N 0 = number of adults in the initial population at generation 0; this number equals the total number of individuals (adults plus larvae) at any time or generation K ³ 0; q 0 = initial frequency of the recessive allele a [q 0 = f 0 (a)]; p 0 = initial frequency of the normal allele A [p 0 = f 0 (A) = 1-q 0 ]; T = delay in larval development determined by the genotype aa, measured in an integer number of generations; K = index corresponding to an integer number of generations; N K = number of adult individuals in the population at the K-th generation; N' K = number of larvae in the population at the K-th generation; q K = frequency of the recessive gene a among adults at the K-th generation; p K = frequency of the normal allele, also among adults belonging to the K-th generation; Q K = frequency of the recessive gene in the total population (adults plus larvae) at the K-th generation; P K = frequency of the normal allele, also in the total population, at the K-th generation (P K = 1 -Q K ).

Derivation of pertinent equations
Based on the foregoing definitions: (a) the total number of adults at generation K is given by (b) the genotypic proportions among AA, Aa and aa adult individuals at generation K are, respectively: (N K-1 /N K )p K-1 2 , (N K-1 /N K )2p K-1 q K-1 , and 0 if K £ T or (N K-1 /N K )p K-1 2 , (N K-1 /N K )2p K-1 q K-1 , and (N K-1-T /N K )q K-1-T 2 if K > T; (c) the number of larvae aa present at generation K after the maturation of normal individuals is given by Using the foregoing equations, it is not difficult to show that N K+1 + N' K+1 = N K + N' K = ... = N 0 + 0 = N 0 and that, if P K = f(A) in the total population, then P K = (N K /N 0 )p K = (N K /N 0 )(N 0 p 0 /N K ) = p 0 , i.e., allele frequencies and the number of individuals remain constant when the whole population (adults plus larvae) is considered, as assumed by the model.
For K £ T, we have, among adults of the population, a mechanism similar to the classic model of total selection against homozygous recessive individuals aa, despite the absence of selection against any individual of the population (when adults and larvae are considered together), as shown by the equations derived in the previous paragraph. This results in q K+1 = q K /(1+q K ), a recurrence equation that has the general solution in simple analytical form q K = q 0 /(1+Kq 0 ).

Dynamic behavior of allele frequencies
Below, we analyze the behavior of allele frequencies p K and q K .
Analysis of the function f(x) = x + 1/x in the interval 0 < x < 1 shows that f'(x) = 1 -1/x 2 < 0; f(x) is therefore a monotonically decreasing function of x in this interval.
If p K > p K-T , then we have p K + 1/p K < p K-T + 1/p K-T = p K + 1/p K+1 and, therefore, p K+1 < p K ; similarly, p K < p K-T implies that p K+1 > p K , and p K = p K-T implies that p K+1 = p K .
In the interval 0 < K £ T, p K increases such that when K surpasses T, we have p K > p K-T ; p K then decreases to a new point where p K < p K-T is reached and we will then have another increase in p K , and so on. Thus, there will be an oscillation in gene frequencies among adults of the population. With the passing of generations, i.e., when K increases, the variations become smaller. In the limit case, when K tends to infinity, the gene frequencies reach an equilibrium (p inf = p e , q inf = q e ). Figures 1-3 show numerical examples of oscillatory convergence to equilibrium points (p e ) over 10 generations of random crosses, for populations starting with initial gene frequencies p 0 = 0.1, ... , 0.9 and T values of T = 1, 2, and 3. These graphs and others shown in this paper were prepared using routines and packages from Mathematica â v.8.0 (Wolfram Research).
The preliminary global analysis presented above already shows a point of real theoretical interest, namely, that despite intense selective pressure among adult individuals (but not among the total population of adults plus larvae), the recessive allele is never eliminated from the population, but rather is kept in equilibrium at a point q e > 0. This is true except when the number of generations involved in the developmental delay is very large or tends to infinity such that allele a is practically eliminated from the adult population and all larvae fail to develop.

Equilibrium gene frequencies
In a generic instant we have the following population composition (Figure 4).
If we consider an equilibrium situation where N K+1 = N K = N K-1 = N K-2 = ... = N K-T = N e q K+1 = q K = q K-1 = q K-2 = ... = q K-T = q e = 1-p e , from N K + N' K = N 0 it follows that N e + TN e q e 2 = N 0 . Since N e = N 0 p 0 /p e , it follows that N 0 p 0 (1+Tq e 2 )/(1-q e ) = N 0 , and, hence, p 0 + Tp 0 q e 2 = 1-q e and, finally,    As expected, q e = q 0 if T = 0 (simple case of panmixia without selection) and q e = 0 if T = ¥ (special case of complete selection against recessive individuals aa).
Considering that: p K = N K-1 p K-1 /N K = N 0 p 0 /N K for any K, N K = N K-1 (1-q K-1 2 ) for any K £ T, and N K = N K-1 (1-q K-1 2 ) + N K-T-1 q K-T-1 2 for any K > T, it is not difficult to show that the recurrence equation 1/p K+1 = p K-T + 1/p K-T -p K can be rewritten in the more convenient form

1/p K = 1/p 0 -S [j = 1, T] {(1-p K-j ) 2 /p K-j }.
This recurrence equation allows straightforward derivation of the equilibrium point p e . Indeed, when K tends to infinity, the limit for the right-hand side of the above equation is 1/p 0 -T(1-p e ) 2 /p e and it follows that p e = p 0 (1+Tq e 2 ) and p e = 1 -q e = 1 -[(1+4Tp 0 q 0 ) 1/2 -1]/(2Tp 0 ). Table 1 shows the sets of equilibrium points q e for several values of T, as a function of the initial values q 0 = 0.0, 0.1, ..., 0.9, 1.0.  For the case T = 1, the recurrence equation above can be rewritten as
Since |r| = |p e 2 -1| < 1 for any value of p e in the interval 0 < p e < 1, the equilibrium point p e is asymptotically stable and convergence to this point occurs at a geometric rate. Since r is always smaller than zero in this interval, convergence to equilibrium is oscillatory. Figure 6 shows the values of r and q e as functions of the initial gene frequency q 0 for the case T = 1.
For a generic value of T the method that follows is used for determining convergence rates as well as equilibrium properties.   and q e (0 < q e < 1) as functions of the initial gene frequency q 0 (abscissa axis).

Equilibrium dynamics
As shown above, for any K > T,
Based on the function f(r) = r T+1 -(a+1)r T + a, r ¹ 1, -1 < a < 0, it can be shown that f(0) = a < 0, f(1) = 0 and df(r)/dr = f'(r) = (T+1)r T -T(a+1)r T-1 . Also, as r tends to +¥, f(r) tends to +inf and as r tends to -¥, f(r) tends to +¥ if T is odd, and to -¥ if T is even. When T is even, f'(r) is an increasing function of r in the open interval (-¥, 0), decreasing in the open interval (0, (a+1)T/(T+1)) and again increasing for r > 0. Therefore, r = 0 is a relative maximum and r = (a+1)T/(T+1) a relative minimum. When T is odd, f'(r) is a decreasing function of r in the open interval (-inf, (a+1)T/(T+1)) and an increasing function of r for r > (a+1)T/(T+1). Therefore, r = 0 in the case when T is odd is an inflection point, whereas r = (a+1)T/(T+1) is a relative minimum (as in the case when T is even). Hence, we conclude that when T is even, the equation r T+1 -(a+1)r T + a = 0, r ¹ 1, -1 < a < 0 does not have any real root since f(r) tends to -¥ as r tends to -¥; when T is odd, the above equation has just one real negative root since f(r) tends to +¥ as r tends to -¥.
For the cases T = 1 to 4 it is possible to obtain explicit solutions in simple analytical form for the values of the eigenvalues r; for T > 4 the eigenvalues can be determined using numerical methods. Analytical procedures applied in the case T < 5 and extensive numerical analysis of the equation for any values of T > 4 have shown that (a) when T is odd, there is always one real eigenvalue, with |ri| < 1 and -1 < ri < 0, and (T-1)/2 pairs of complex conjugates r j,k = A ± iB, with (A 2 +B 2 ) 1/2 < 1 for any initial gene frequency p 0 ; (b) when T is even, there are always only T/2 pairs of complex conjugates with the same properties as stated before. We therefore conclude that the set of equilibrium points p e = f(p 0 , T) is a stable one, with convergence always occurring to it for any initial value p 0 in the open interval (0, 1) and for any integer value of T.

Description of the model in which there are fluctuations in the total number of individuals
In order to study the effect of an increase or decrease in population size, it is enough to introduce the additional parameter c: N 1 p 1 = cN 0 p 0 into equation N 1 p 1 = N 0 p 0 of the previous section. Consideration of this parameter leads straightforwardly the equation p K+1 = 1/[2-p K +(1-p K-T ) 2 / (c T p K-T )], K > T. This is the same equation obtained in the analysis of the constant population size model if c = 1.
At equilibrium, p K+1 = p K = p K-T = p e and the above equation reduces to p e (1-p e ) 2 (1-c T ) = 0, c ¹ 1. When c ¹ 1, the only possible solutions for the foregoing equation are p e = 0 and p e = 1, i.e., the normal allele is completely lost or fixed. Intuitively, it is easy to show that p e = 0 when c < 1 (decreasing population size). Thus, for any K > T, the proportion of recessive genes introduced by larvae ending their diapause period will become increasingly greater as generations go by, in relation to the normal allele present in the adult population, the size of which is decreasing. The in-verse, i.e., when c > 1 (population increasing in size), is also obviously true since proportionally fewer recessive genes will be introduced into the adult population by emerging larvae. Clearly, then, p e tends to 1 as generations go by.

Brief remarks on a model with continuous generations
Let P(t), R(t), Q(t) and L(t) be the normalized proportions {P(t) + R(t) + Q(t) + L(t) = 1} of adults AA, Aa and aa, and larvae aa in generation t. Among adults, the normalized genotypic proportions are given by: In the total population (adults + larvae), the allelic frequencies are given by p(t) = P + R/2 and q(t) = Q + L + R/2, p(t) + q(t) = 1.
After a time interval Dt, the frequency of allele A among adults of the population is given by
The expression for p(t+Dt) can be used for quantitative analysis (to be presented elsewhere) in the case of continuous generations. One expects that equilibrium values for any fraction 0 < e < 1 in T < T+e < T+1 would interpolate between the equilibrium curves that represent the set of equilibrium points for the integer cases of T and T+1 generations in delayed larval development.