Toward Reconciling Inferences Concerning Genetic Variation in Senescence in Drosophila melanogaster

Corresponding author: Frank H. Shaw, Department of Ecology, Evolution and Behavior, University of Minnesota, 1987 Upper Buford Circle, St. Paul, MN 55108., fshaw{at}superb.ecology.umn.edu (E-mail)

Communicating editor: P. D. KEIGHTLEY

	ABSTRACT

TOP ABSTRACT Empirical studies of variance... Statistical challenges for the... THE MODEL AND STATISTICAL... RESULTS DISCUSSION APPENDIX 1 LITERATURE CITED

Standard models for senescence predict an increase in the additive genetic variance for log mortality rate late in the life cycle. Variance component analysis of age-specific mortality rates of related cohorts is problematic. The actual mortality rates are not observable and can be estimated only crudely at early ages when few individuals are dying and at late ages when most are dead. Therefore, standard quantitative genetic analysis techniques cannot be applied with confidence. We present a novel and rigorous analysis that treats the mortality rates as missing data following two different parametric senescence models. Two recent studies of Drosophila melanogaster, the original analyses of which reached different conclusions, are reanalyzed here. The two-parameter Gompertz model assumes that mortality rates increase exponentially with age. A related but more complex three-parameter logistic model allows for subsequent leveling off in mortality rates at late ages. We find that while additive variance for mortality rates increases for late ages under the Gompertz model, it declines under the logistic model. The results from the two studies are similar, with differences attributable to differences between the experiments.

WHY do organisms age? Ultimately, senescence must arise because the strength of selection declines with age (HAMILTON 1966 ; CHARLESWORTH 1994 ). We define senescence as the persistent decline in age-specific fitness components (i.e., reproduction and survival) due to internal physiological decline (ROSE 1991 ). However, the theory of age dependence of selection sensitivity does not describe the genetic factors, expressions, or interaction that cause the senescence phenotype. This more mechanistic problem requires specific genetic models for aging, and biologists commonly cite two genetic theories—"mutation accumulation" (MEDAWAR 1952 ) and "antagonistic pleiotropy" (WILLIAMS 1957 )—to account for the causes and genetic architecture of senescence.

Consider a mutation that leads to an increase in mortality rate (µ_t) at time t. [The mortality rate µ_t is the limiting value of age-specific mortality, q_t, as the age interval becomes infinitesimally small (ELANDT-JOHNSON and JOHNSON 1980 ). It is often estimated as the central death rate (LEE 1992 )]. If that gene is expressed early in life (for example, before reproductive maturity), natural selection tends to reduce the frequency of the mutation in the population, whereas, if the effects of the gene are not expressed until late in life, natural selection will have a much weaker impact on the frequency of a similar mutation. Thus, late-acting deleterious genes are more likely to accumulate over evolutionary time than are early acting deleterious genes. MEDAWAR 1952 argued that this age-specific increase in the expression of deleterious mutations leads, in turn, to an increase in age-specific mortality rates. WILLIAMS 1957 suggested that the same late-acting deleterious genes that natural selection fails to remove may actually be selected for if they have beneficial (i.e., "antagonistic pleiotropic") effects early in life.

To test the adequacy of these models, biologists have sought to make specific predictions and to conduct critical experiments. Williams noted that negative associations between late and early fitness traits are expected under antagonistic pleiotropy. If this trade-off arose as a result of optimality selection (PARTRIDGE and BARTON 1993 ), there would be no standing variation for the traits involved, and consequently the traits would not respond to selection. But if these traits are controlled by alternative alleles maintained at a polymorphic equilibrium, then selection upon one trait is expected to produce negatively correlated changes in a second trait (ROSE 1984 ; CHARLESWORTH 1994 ). Studies of age-specific selection upon longevity and age-of-reproduction with Drosophila have generally revealed such correlated selection responses (ROSE and CHARLESWORTH 1980 ; LUCKINBILL et al. 1984 ; ROSE 1984 ; ZWAAN et al. 1995 ). These are consistent with some expectations of the antagonistic pleiotropy model but may be explained by other genetic mechanisms (CLARK 1987 ; CURTSINGER et al. 1995 ) or elements of the selection design (PARTRIDGE and BARTON 1993 ; LEROI et al. 1994 ; PROMISLOW and TATAR 1998 ). Direct measures of age-specific genetic correlation structure are also expected to reveal negative covariances (ROSE and CHARLESWORTH 1980 ). In Drosophila melanogaster negative covariances have been reported between early and late age-specific fecundity (ROSE and CHARLESWORTH 1981 ; ENGSTROM et al. 1989 ), and TATAR et al. 1996 found a negative relationship between fecundity at the earliest and mortality at the latest ages, although all other age-specific pairs of covariances within and between these traits were positive or not different from zero.

Compared to work on antagonistic pleiotropy, there have been fewer predictions and tests specifically aimed at addressing the relevance of mutation accumulation to senescence, although several novel expectations have recently been proposed (TANAKA 1993 ; CHARLESWORTH and HUGHES 1996 ). The most commonly discussed prediction is that additive genetic variance for fitness traits should increase with age. In addition, CHARLESWORTH and HUGHES 1996 have recently pointed out that both inbreeding load and dominance variance for survival (not mortality) rates should increase with age under mutation accumulation but not under antagonistic pleiotropy. To date, the first prediction, that addi-tive variance should increase with age, has been most widely explored. Although there was something of a hiatus after initial tests of this prediction (ROSE and CHARLESWORTH 1980 , ROSE and CHARLESWORTH 1981 ; KOSUDA 1985 ; ENGSTROM et al. 1989 ), in recent years the prediction has once again become the focus of a concerted research ef-fort. At least one test was made of the prediction that inbreeding loads and dominance variance should increase with age under the mutation accumulation model (CHARLESWORTH and HUGHES 1996 ). They found that inbreeding load increased significantly with age for male mating success and for age-specific survival estimates. Dominance variance increased significantly with age for age-specific survival but not for mating success.

Two issues have emerged with respect to the variance trajectory prediction. First, formal models reveal that both antagonistic pleiotropy and mutation accumulation should lead to an increase in additive genetic variance (CHARLESWORTH 1990 , CHARLESWORTH 1994 ; PARTRIDGE and BARTON 1993 ; CHARLESWORTH and HUGHES 1996 ). Additive genetic variance at late relative to early ages approaches infinity as the age-specific function of selection sensitivity approaches zero (CHARLESWORTH and HUGHES 1996 ). Thus, this prediction is a test of the validity of current models for both explanations rather than a differential diagnostic. Second, concerted empirical studies were launched to evaluate the prediction in terms of age-specific reproductive traits and mortality. In the first study of its kind HUGHES and CHARLESWORTH 1994 ; HUGHES 1995 ; CHARLESWORTH and HUGHES 1996 concluded that additive genetic variance for fitness traits increased monotonically with adult age. PROMISLOW et al. 1995 soon after confirmed that additive variance for mortality and female fecundity increased initially with age but discovered that the age-specific additive variance component for mortality eventually declined. It is our purpose here to understand and resolve the discrepancy in these reports concerning mortality. We begin with a short description of the two studies to highlight the demographic, genetic, and statistical issues. We then present a novel statistical technique that solves fundamental problems in the original analyses.

	Empirical studies of variance for age-specific mortality

TOP ABSTRACT Empirical studies of variance... Statistical challenges for the... THE MODEL AND STATISTICAL... RESULTS DISCUSSION APPENDIX 1 LITERATURE CITED

HUGHES and CHARLESWORTH 1994 used quantitative genetic methods to estimate the genetic variance components for age-specific mortality at three different ages in a population of male D. melanogaster. The experimental design made use of 40 genetic lines each isogenic for a different wild-type third chromosome. These lines were divided into five independent blocks. Crosses of the lines were performed within each block, yielding 80 sib groups, which were observed in a total of 304 subcohorts of 20 males each. Approximately 6080 deaths were observed at weekly intervals over 14 wk. For mortality between eclosion and 3 wk and from 5 to 7 wk, genetic variance components were not significantly different from zero. At 9–11 wk posteclosion, in each of five independent blocks, Hughes and Charlesworth found substantial additive and dominance variance for age-specific mortality (reported in CHARLESWORTH and HUGHES 1996 ). Initially these data were argued to be consistent with the mutation accumulation theory for the evolution of senescence (HUGHES and CHARLESWORTH 1994 ). Further analysis based on inbreeding load was used to support the original interpretation of senescence evolving via mutation accumulation (CHARLESWORTH and HUGHES 1996 ).

PROMISLOW et al. 1995 conducted an experiment similar to that of Hughes and Charlesworth but involving substantially more flies. In the Promislow et al. experiment, a single-block design with 10 genetic lines crossed to create 25 heterozygous genotypes was employed. Each of these genotypes was observed in 4 cages with ~320 flies of each sex per cage, for a total of 100 cages containing 65,134 flies. Deaths were observed every day over 65 days; mortality rates were calculated on a weekly basis. Promislow et al. also found that the additive genetic variance for log mortality rate, V_A(ln(µ_t)), markedly increased through reproductive maturity for females and significantly increased but less so for males. But contrary to the results of Hughes and Charlesworth, at late ages, V_A(ln(µ_t)) declined for both males and females.

The pattern of declining variance observed by Promislow et al. runs counter to any expectation derived by existing mathematical models of either antagonistic pleiotropy or mutation accumulation. On the basis of current assumptions, both models predict an increase in additive variance for fitness traits with age. This increase is expected to reach a maximum at the age of last reproduction when the sensitivity of selection becomes zero, and the variance should remain at this level, or potentially increase, through all postreproductive ages where the sensitivity remains at a minimum (CHARLESWORTH and HUGHES 1996 , Table 1 and Equation 6). Promislow et al. proposed five possible explanations for their results: (1) temporal changes in the laboratory environment, i.e., period effects; (2) effects of reproduction and its associated costs; (3) heterogeneity of individual frailty; (4) bias at advanced ages when sample size is diminished; and (5) individual age-dependent genetic processes specific to somatic senescence. Our present purpose is to sort out the differing observations of Hughes and Charlesworth and Promislow et al. and their rather different sets of associated implications.

	Statistical challenges for the study of mortality variance

TOP ABSTRACT Empirical studies of variance... Statistical challenges for the... THE MODEL AND STATISTICAL... RESULTS DISCUSSION APPENDIX 1 LITERATURE CITED

Although there are many differences between the studies of Hughes and Charlesworth and Promislow et al., the basic idea behind them is the same. Cohorts of Drosophila lines among which the genetic relationships are known are observed from eclosion until mortality. Mortality is observed at frequent intervals for each subcohort, and variance component analysis is performed on discrete time mortality rates estimated from these observations. Although this analytical approach to the data is pragmatic, there are several problems, all having to do with its treating the mortality rate, µ_t, as having been directly observed when only an estimate of µ_t is obtained. If n flies in a subcohort are alive at the beginning of a certain time interval and k of them die during the interval, then k is a binomial random variable with sample size n and success probability p_t = exp(-µ_t). So the estimate of ln(µ_t) is ln(-ln(k/n)). At early ages when µ_t is small, there is high probability that k = 0, in which case the estimate of ln(µ_t) is undefined, and at late ages when n is small there are only a few possible values of k and the estimate is very crude. At very late ages, when some genotypes have died out altogether (so n = 0), there is no estimate of ln(µ_t) for those genotypes.

These problems were handled in different ways in the original analyses of the Hughes and Charlesworth and Promislow et al. data sets. Hughes and Charlesworth, with subcohort size 20, grouped several age classes together to obtain an average mortality rate for their variance component analysis and, in their analysis, which fitted the Gompertz model to their data, used ln(µ_t + 1) in place of ln(µ_t) so that this variable would be defined when µ_t = 0. This latter transformation changes the statistical model and necessarily reduces V_A(ln(µ_t + 1)) at early ages. At very late ages, Hughes and Charlesworth did not estimate V_A(ln(µ_t)) because of loss of lines. The Promislow et al. experiment, designed in part to reduce the sample size problem apparent in the earlier Hughes and Charlesworth experiment, increased subcohort size to >300. To increase the demographic power, Promislow et al. sacrificed quantitative genetic power; they observed fewer genotypes. Mortality rates of zero, which were few, were treated as missing data, which also introduces some bias.

The Hughes and Charlesworth and Promislow et al. results were obtained by variance component analysis of estimated ln(µ_t + 1) or ln(µ_t) at discrete ages. Such analyses assume the phenotypic variable (in this case log mortality rate) is measured accurately. Because it was not, and in fact was very crudely estimated at both early and late ages, the assumptions of the standard quantitative genetics model are badly violated. According to statistical theory the correct approach is to use a model that properly accounts for the fact that mortality rates are only estimates. In principle, this is easily done. The log mortality rates ln(µ_t) are considered an unobserved random vector obeying a standard quantitative genetics model. The conditional distribution of deaths given µ_t is an independent binomial in each cohort-age class. The parameters of the model can be estimated by maximum likelihood. In practice, such an analysis is very difficult. Although such models are beginning to attract some research interest among statisticians, the estimation methods proposed have only been applied to much simpler data than we have here, and some proposals are only crude approximations to maximum likelihood. We used Markov chain Monte Carlo (MCMC) maximum likelihood (THOMPSON and GUO 1991 ; GEYER 1994 ) to maximize the likelihood for this very complex model by brute-force computation.

Because the variable of interest ln(µ_t) is not observed and only crudely estimated at early and late ages, treating ln(µ_t) nonparametrically with a different latent variable for each cohort-age class results in too many covariance component parameters to be estimated with precision. Thus we adopted a parametric model for the mortality curve. This parameterization makes the analysis feasible but necessarily will present substantial limits on our inferences, as we shall see.

Two such models were employed in the original analyses of the Hughes and Charlesworth and Promislow et al. data, which we reconsider here. These are the GOMPERTZ 1825 model and the more recent logistic model (VAUPEL and YASHIN 1985 ; VAUPEL 1990 ). The Gompertz model has µ_t, the mortality rate at age t, given by

(1)

rewritten

(2)

where B = ln(C). We call B the Gompertz intercept parameter and A the Gompertz slope parameter. The logistic model was developed to take account of the differing levels of frailty exhibited by individuals who are assumed to have the same underlying Gompertz mortality rate. An individual's mortality risk is then zµ_t, where z is its frailty. We assume individual frailty within each subcohort has a gamma distribution with mean 1 and variance at birth. If we subsume the frailty z into µ_t we arrive at the expression

(3)

which has three parameters, the two underlying Gompertz parameters C and A and a variance of individual frailty parameter . The special case = 0 gives the Gompertz model (1), but for > 0, the mortality rate levels off at late ages, a phenomenon observed in large-scale studies (CAREY et al. 1992 ; CURTSINGER et al. 1992 ).

In much recent work (FUKUI et al. 1993 ; PLETCHER and CURTSINGER 1998 ; PROMISLOW and TATAR 1998 ) the logistic model (3) has been treated as a three-parameter model capable of fitting logistic mortality curves. In this work, the parameter is referred to as the "deceleration" parameter because its value controls the extent to which mortality-rate curves level off. The original intent of the model has recently been reemphasized (SERVICE et al. 1999 ).

Using regression analysis of mortality rate on age (Hughes and Charlesworth) or maximum-likelihood methods (Promislow et al.), both studies found significant genetic variance for the Gompertz slope parameter. Only in Promislow et al. was significant genetic variance for the intercept parameter detected, and then only for males. In both experiments, estimates for the covariance between slope and intercept were reported to be negative but also not signficantly different from zero (HUGHES 1995 ; PROMISLOW et al. 1995 ). The logistic model was fitted using maximum likelihood to the Promislow et al. data and was found to have a significantly better fit than the Gompertz; with this model significant additive variance was detected for the slope, intercept, and parameter for both sexes. All of these analyses, however, had the flaw of assuming that the µ_t were observed when they were actually only estimated.

In a new analysis reported here, the data from the Hughes and Charlesworth and Promislow et al. experiments are modeled to take account of the chance variability in the observed deaths. Means and covariance component matrices are estimated for the parameters of the logistic model. Cohort size, though it does affect the power of the analysis, does not seriously bias the estimates for early and late life V_A(ln(µ_t)). In this new analysis, the datasets yield similar results in that an increase in additive variance is observed. The data of Promislow et al. continue to reveal a decline in additive variance from mid- to advanced ages but the data of Hughes and Charlesworth lack adequate power to resolve this phenomenon. The variance between these experiments may be attributable to both biological and technical differences.

	THE MODEL AND STATISTICAL METHOD

TOP ABSTRACT Empirical studies of variance... Statistical challenges for the... THE MODEL AND STATISTICAL... RESULTS DISCUSSION APPENDIX 1 LITERATURE CITED

The Gompertz model of mortality (1) is derived from the observed tendency for mortality rates to increase exponentially with age in many organisms. We can then define a mortality rate trait y = (A, B), where A and B are the Gompertz slope and intercept parameters, and assume that y has a multivariate normal distribution with mean ß (a 2-vector) and variance V (a 2 x 2 matrix). We are interested in the genetic and phenotypic variation to be found in y in a population. As usual in quantitative genetics, we partition y into the independent genetic and residual components y = g + r. The covariance of breeding values g for different genotypes is proportional to the numerator relationship matrix, and the covariance of residual values r for different genotypes is zero.

This notion extends naturally to more complicated models. For the logistic model (3) we define y = (A, B, ) with B = ln(C) as in the Gompertz model, and as before we assume that y has a multivariate normal distribution with mean ß (now a 3-vector) and variance V (now a 3 x 3 matrix).

Because (A, B, ) are not observable, it is impossible to assess the validity of the assumption that they should be multivariate normally distributed. Promislow et al. did examine the normality of their estimates of ln µ_t for discrete age classes and found these to be normally distributed.

Change in additive genetic variance with age:
Although the primary interest in the data we are analyzing is whether the additive genetic variance for the mortality trait V_A(ln(µ_t)) increases or decreases with age, neither the Gompertz nor the logistic model has this change of variance with time as an explicit parameter. The variance V of the "phenotypic" vector y does not change with time. We have to do a calculation to see what the detailed form of each model implies about the variance of ln(µ_t) as a function of age.

In the Gompertz model with ln(µ_t) = A + Bt, the phenotypic variance for mortality rate can be expressed as a function of age:

(4)

This is a quadratic function of t, and the coefficient of t² is positive because it is a variance. V(t) is the variance matrix of y = g + r at age t. Because g and r are independent, we have V(t) = V_A(t) + R(t), where V_A(t) is the variance of g and R(t) the variance of r. Both V_A(t) and R(t) are given by equations similar to (4) and so are quadratic in t with a positive coefficient of t².

Hence V(t) and V_A(t) and R(t) increase with t for large t. It is important to understand that this is a built-in feature of the Gompertz model. The Gompertz model cannot model a decrease in V_A(t) with age at late ages unless it decreases for all ages of interest. Thus an increase in V_A(t) at late ages is only the expected "built-in" behavior of the Gompertz model.

An analogous analysis for the logistic model is not possible because of the complexity of the functional form (3) of the specification of µ_t. However, we can get some idea about the behavior of the variance of ln(µ_t) for large t by noting that ln(µ_t) converges to an asymptote ln(A) - ln() for large t. Hence ln(µ_t) is asymptotically constant in t and so is its variance. Whether the variance increases or decreases to this constant is not obvious. If A is highly correlated with , a decrease in variance would result. We observed such a correlation and decrease of additive genetic variance with age at late ages in the Promislow et al. data. It seems likely that this is also a built-in behavior of the model.

The behavior of both models for "large t" may not be relevant because this may refer to times beyond which all flies have died. The fact that mortality variance under the Gompertz model cannot increase and then decrease, however, is a significant constraint.

Recent theoretical work (CHARLESWORTH and HUGHES 1996 ) has emphasized the importance of dominance variance for log mortality rate in distinguishing between the mutation accumulation and antagonistic pleiotropy hypotheses of senescence. The addition of a third covariance matrix of parameters to those of the additive and residual components was beyond our computational capabilities because shortcuts to likelihood calculation (see Appendix 1) would then not have been possible. The presence of dominance variance in these data, documented by both Hughes and Charlesworth and Promislow et al., will bias upward both the additive genetic (which should nevertheless not be taken as general genetic variance) and the residual variance components. Because dominance variance as a proportion of total phenotypic variance was found in earlier analyses to be small, we feel that its absence from our current model is of little consequence to the questions we are attempting to answer; e.g., Is the additive genetic variance for the log of late life mortality risk in these datasets increasing or decreasing?

Monte Carlo maximum likelihood:
The probability of survival at age t is written

for the Gompertz model. Note that this depends on the parameters ß and V of the quantitative genetics model only through y = (A, B). For the logistic model we still have p_t = exp(-µ_t), but plugging in (3) for µ_t gives a more complicated function involving y = (A, B, ). It is still true that the distribution of the observed data n given y depends only on y and does not involve the parameters of the quantitative genetics model.

Let n be a vector of census numbers from k ages, so that n_t is the number of individuals surviving at age t. Because the number of individuals surviving over a given period beginning with t - 1 and ending with t is binomially distributed, we can write the likelihood of n_t-1 - n_t individuals not surviving over the period as being proportional to

The conditional density for n given y is

and the complete data likelihood function is

(5)

where the distribution of y is multivariate normal

(6)

where is the parameter vector (ß and V are functions of ) and neither C₁ nor C₂ contains . If there is no mortality during an interval in a given subcohort (n_t-1 = n_t), we have no contribution to f(n, y) due to that subcohort in that interval. This eliminates problems in early age classes when mortality rates are low. Similarly, at late ages, f(n, y) is unaffected by subcohorts that have completely died out (n_t = 0).

To get the likelihood of a given vector n we must integrate out y from (5):

Although this cannot be done analytically or by numerical integration, the integral can be done using Monte Carlo integration and the resulting likelihood function can be maximized (THOMPSON and GUO 1991 ; GEYER 1994 ).

Fix a point ₀ in the parameter space. For reasons of computational convenience we use the likelihood ratio as the likelihood for the problem. Then

The last expression is a conditional expectation. Thus if we can generate random variates Y_i having the conditional density f₀(y|n), we can evaluate the likelihood ratio as

(7)

For large K this is close to what we would have obtained if we could have done the integrals exactly. Note that this approximates the likelihood for all parameter values using a sample from the conditional distribution of y, given the observed data n for one fixed parameter value ₀. The approximation will be good only in a neighborhood of ₀, so ₀ must be close to the maximum-likelihood estimate (MLE). This may require some preliminary trial and error (see GEYER and THOMPSON 1992 for more discussion of this issue).

The distribution of Y is multivariate normal, but the conditional distribution given the observed data n required here is very complicated, so complicated that it is not possible to obtain independent samples Y_i. However, we can easily obtain a Markov chain sample using the Metropolis algorithm (METROPOLIS et al. 1953 ). Let Y_k denote the vector of A, B, and values for each genotype. We generate a multivariate normal random vector Z_k with mean zero and fixed variance matrix V₀. Then the next point in the sample after Y_k is either Y_k+1 = Y_k + Z_k or Y_k+1 = Y_k, depending on the value of

With probability min(1, R) we "accept" the proposed change and set Y_k+1 = Y_k + Z_k. Otherwise we set Y_k+1 = Y_k and we have two consecutive samples that are the same. This algorithm generates a Markov chain Y₁, Y₂, ... , having the distribution of interest f₀(y|n) as its equilibrium distribution. Verification of stronger convergence properties of this Markov chain, such as geometric ergodicity, was not done because such calculations are very difficult and highly technical except for very simple models. However, we can expect that for large-enough K, the Monte Carlo approximation to the likelihood (7) will be reasonably accurate.

In all our Markov chain sampling, the matrix V₀ was taken to be diagonal so that independent normal changes were made to each parameter. A nondiagonal V₀ matrix can give more efficient sampling, but it is necessary to already have good samples to determine which V₀ is best, and at that point there is no reason to go back and get a better sample. The diagonal elements of V₀ were determined by trial and error to give about a 20% acceptance rate in the Metropolis algorithm, which some literature suggests is a good target (GELMAN et al. 1996 ).

If we use the fact that f(n|y) does not depend on the parameters we obtain

because the terms f(n|Y_k) in the numerator and denominator cancel. Using the specific form (6) of f(y), we get

where now denotes only the variance and covariance component parameters involved in V, where the mean parameters have now been written in regression form Xß, where X is a fixed known matrix, and we have switched to lowercase for the Markov chain sample y₁, y₂, ... Taking logs we get

(8)

for the log likelihood, where

This quantity can be calculated very rapidly for a given set of parameter values and ß by means of eigenvalue decomposition of the numerator relationship matrix (see Appendix 1) and used in the above expression to maximize the Monte Carlo approximation to the log likelihood (8).

The analysis thus proceeds in two stages. The first stage involves generating, via the Metropolis algorithm, points y₁, ... , y_K that are a representative sample from the conditional distribution of y given the observed data n and parameter value (₀, ß₀). With the sample recorded in a computer file, we proceed to the second stage of the analysis, the maximization of (8) to find MLEs, and . Maximization of (8) was carried out using the quasi-Newton maximizer in the Numerical Algorithms Group library, E04UCF (separately available under the name NPSOL from Stanford University).

One of the virtues of the Monte Carlo approach is that the observed Fisher information, which is the negative of the matrix of second partial derivatives of the log likelihood (the negative of its Hessian matrix), is well approximated by the negative of the Hessian of (8), which is calculated by the optimization software. Standard errors were calculated using the diagonal of the Fisher information matrix and confirmed using likelihood profiles wherever possible.

When the covariance matrix of estimates became infeasible (i.e., had negative eigenvalues), the objective function (8) could not be calculated. Consequently, in the Hughes and Charlesworth data, the parameters making up the covariance matrix were transformed to the elements of the triangular matrix of the square root of the covariance matrix (MEYER and SMITH 1996 ), where constraint to feasibility (diagonal elements positive) was easy to enforce.

The Hughes and Charlesworth data consisted of five independent blocks with weekly census numbers for 60–68 subcohorts per block over 16 wk. The same Markov chain sample was used for all the analyses of the Hughes and Charlesworth data (different parameters for each block and a pooled analysis with the same parameters for each block). A Monte Carlo sample size of K = 10,000 was used. A spacing of 2000 iterations was used because we expected slow mixing of the Markov chain, but the results were much the same when a spacing of only 10 iterations was tried.

The Promislow et al. data resembled a single block of the Hughes and Charlesworth data but with census numbers every 5 days through day 65. The Markov chain sampling was much easier for the Promislow et al. data because the peak in the likelihood function was much more pronounced than in the Hughes and Charlesworth data. We therefore used spacing of one iteration; i.e., each iterate was recorded.

A multivariate normal distribution of missing data (Y) is assumed for the likelihood function, as is common in quantitative genetics. Because these data are unobserved, there is no way to tell if this assumption is violated or not. The discrete time estimates for mortality rate used in both the Hughes and Charlesworth and Promislow et al. original analyses were normally distributed.

	RESULTS

TOP ABSTRACT Empirical studies of variance... Statistical challenges for the... THE MODEL AND STATISTICAL... RESULTS DISCUSSION APPENDIX 1 LITERATURE CITED

The results of these analyses are in the form of estimates of mean vectors and covariance component matrices for the parameters of the logistic model. To determine the functional form of V_A(t) for the logistic model, we again used Monte Carlo integration. The Monte Carlo sample contains A_k, B_k, and _k values for the genetic component of these quantities for each subcohort-age class with k running over the Monte Carlo sample. Variances for log mortality rate were calculated for each day through the lifespan of the longest-lived flies for each simulated dataset. The mean values for these variances are plotted in Figure 1 Figure 2 Figure 3 Figure 4 Figure 5. To assess the distribution of the mean genetic variance for mortality rate as a function of time, the above process was repeated using 1000 samples from the distribution of the estimated parameters

where is the MLE of parameter means and covariance matrices and I() is the Fisher information matrix calculated from the Hessian from the optimization software. This process was not practical for the Hughes and Charlesworth results because most of the samples from the distribution of the parameters were not feasible; in particular the sampled variance for Gompertz slope was negative for half the samples.

View larger version (46K): In this window In a new window Download PPT slide   Figure 1. Additive genetic variance trajectories over time for the five independent blocks of the Hughes and Charlesworth experiment. Variance at a given time was determined by simulating 1000 blocks of cohorts having the additive genetic covariance matrices estimated in the analysis. This provided, for each simulated block, a distribution of genetic effects for log mortality rate at each age class. The mean variance of these, over the 1000 blocks, is plotted.

View larger version (15K): In this window In a new window Download PPT slide   Figure 2. Additive genetic variance trajectories as a function of percent mortality for the five independent blocks of the Hughes and Charlesworth experiment. Values for VA (ln(µt)) were calculated as in Figure 1.

View larger version (9K): In this window In a new window Download PPT slide   Figure 3. Additive genetic variance trajectories over time for males and females in the Promislow et al. experiment. Values for VA(ln(µt)) were calculated as in Figure 1.

View larger version (9K): In this window In a new window Download PPT slide   Figure 4. Additive genetic variance trajectories as a function of percent mortality for males and females in the Promislow et al. experiment. Values for VA(ln(µt)) were calculated as in Figure 1.

View larger version (10K): In this window In a new window Download PPT slide   Figure 5. Additive genetic variance trajectory as a function of time for females in the Promislow et al. experiment plotted with 95% confidence limits. Variance trajectories were calculated as in Figure 1 for 1000 different covariance matrices sampled from the asymptotic distribution of .

Hughes and Charlesworth data:
The Hughes and Charlesworth data for 7-day counts of mortality for the 304 genotypes divided into five blocks were analyzed under the logistic model. The results are shown for the separate blocks with 12 variance components and three block means each (Table 1 and Table 2, Figure 1 and Figure 2).

Under the logistic model, V_A(ln(µ_t)) climbs or is constant over the lifespan of the flies (Figure 1 and Figure 2). Because the estimates under the logistic model for the additive genetic variance components pertaining to were vanishingly small, feasibility of this covariance matrix was a problem in both the pooled and the separate block analyses. When this covariance matrix was singular, it was impossible to back-transform the Hessian matrix developed by the optimization software from the square root values (see above) to calculate standard errors of the estimates. We therefore removed these components from the final model as reported in Table 2. The removal of these components from the model resulted in negligible likelihood change. Note that the residual covariances and, more significantly, the block means for remained in the model. Though the variance of is not detectable, the mean estimate is significantly greater than zero, indicating that the logistic model fits these data better than the Gompertz model where = 0.

The separate block analyses, for which 45 covariance components were estimated, had a combined maximum-likelihood ratio 42 points higher than the pooled analysis (not shown), which has 9 covariance components. This demonstrates that there are significant differences among the covariance component matrices (P < 0.001). This test included residual covariance matrices, which appear to be mostly responsible for the differences.

Genetic variance for the intercept parameter is the only significant estimate in most of the genetic covariance matrices. Both the covariance between slope and intercept and the variance for slope itself are always less than a standard error from zero. From this it is evident that, although the estimates point to a slow increase in genetic variance, a much more rapid increase or even a decline cannot be ruled out.

We can also plot the trajectories for genetic variance as a function of percentage mortality (Figure 2). This demonstrates that much of the steep decline at very late ages seen in Figure 1 is taking place after the vast majority of the flies in a given block have died. The steep declines are extrapolations following the logistic model using parameters estimated from earlier age classes when mortality rates could be estimated with more precision. Because V_A() is not detectable and is here set to zero, V_A(ln(µ_t)) has a constant asymptotic value, V_A(A)/ with the block mean . The additive genetic variance for Gompertz slope, V_A(A), is itself very small and not significant for these data, so the late-age asymptotic values for V_A(ln(µ_t)) are small in comparison to the additive genetic variance for Gompertz intercept, V_A(B), which is significant and dominates V_A(ln(µ_t)) at early ages.

Promislow et al. data:
The data for 5-day counts of mortality on 65,134 male and female flies in 100 cages were analyzed under the logistic model. The first two observations (days 5 and 10) were discarded because there is a high degree of variability in these observations and interest was in the middle and late stages of life. Means and covariance matrices given are for model parameters per 5 days beginning at age 10 for males and females (Table 3).

As was the case in the Hughes and Charlesworth data, the standard error for the mean indicates that it contributes significantly to the model. In the Promislow et al. data, we have significant levels of other variance components including V_A() as well. The genetic correlation between Gompertz slope and is positive and significant for both males (r = 0.77) and females (r = 0.96). This leads to a decline in variance at late ages, because the mortality curves are asymptotic to A/, which, if A and are highly correlated, varies little among cages. Interestingly, the residual correlations between A and (0.58 for males, 0.54 for females) are lower and are not even significantly different from zero for the females.

Figure 3 and Figure 4 show the trajectory of genetic variance under the logistic model as a function of time and percent mortality. The peak variance, at day 38, coincides with the time that egg laying ceased and also with the time when about half the flies had died. By simulating the mortality parameters from the sampling distribution and plotting a large number of them, we were able to show that variance V_A(ln(µ_t)), which is a complicated function of these parameters, has a significant decline at late ages. A 95% confidence region for the female V_A(ln(µ_t)) is shown in Figure 5. Of the 1000 trajectories plotted, 989 of them were declining after day 45.

	DISCUSSION

TOP ABSTRACT Empirical studies of variance... Statistical challenges for the... THE MODEL AND STATISTICAL... RESULTS DISCUSSION APPENDIX 1 LITERATURE CITED

Each of the datasets reanalyzed here was generated to gain insight into the pattern of genetic variance for mortality rate at different ages. The original analyses reached contradictory conclusions. On the basis of discrete analyses and fitting of the Gompertz model to the data, Hughes and Charlesworth concluded that genetic variance for log mortality rate increases late in life. On the basis of discrete analysis of a larger dataset, Promislow et al. concluded that the genetic variance decreases at late ages. Their conclusion calls into question the primary assumptions underpinning both the mutation accumulation and antagonistic pleiotropy models of senescence. The analysis we present here relies upon parameteric mortality models to estimate simultaneously demographic and genetic parameters. In this context we largely substantiate the earlier findings, but we document a crucial difference in power between the two experiments. Genetic variance for log mortality rate increases through most of the lifespan of the flies in the Hughes and Charlesworth experiment, certainly well beyond the point where the sensitivity of selection with respect to fitness, S(w), is zero, but the increase is not statistically significant. In the Promislow et al. experiment, V_A(ln(µ_t)) increases until the end of reproduction and then declines; both the increase and the decline are significant.

Although the Hughes and Charlesworth data do not exhibit a significant peak and decline of variance, standard errors on the parameter estimates are large because of the small subcohort sizes, and a pattern like that seen in the Promislow et al. data cannot be ruled out. Because the genetic variance of could not be detected, it is not included in the trajectories of Figure 1 and Figure 2. The absence of in the model may account for the absence of decline in overall genetic variance until extrapolation beyond the final deaths. In Hughes and Charlesworth's original discrete analysis, variance estimates were not made for very late ages (>12 wk) because loss of whole subcohorts might have biased the results. Thus, estimates for late ages where mortality curves that might have yielded declining genetic variance were not available. The discrete age estimates of the genetic variance at early ages provided little information because few deaths were available from which to estimate the mortality rate. The constraints of demographic power in the discrete analysis extended as well to the original parametric analysis. The Gompertz model based on regression was fit to mortality curves beginning at early ages when there was no detectable variance, and this was bound to lead to the conclusion that genetic variance for log mortality rate is increasing at a highly significant rate. In the current analysis, significant variance is found at the intercept but variance for slope is not detectable. Therefore, in our estimated plots of genetic variance for mortality (Figure 1 and Figure 2), the observed increases in V_A(ln(µ_t)) are not significant. While genetic variance for was also not detectable, the mean of was significantly different from zero, and the dependence of overall genetic variance on this mean value in the logistic model, along with the small estimate for slope variance, results in an extrapolated decline in V_A(ln(µ_t)) but at ages well beyond those of any interest. In sum, we observe an increase in genetic variance for mortality in three of five blocks, but the power is insufficient to infer whether these increases are statistically significant, and whether subsequently variance declines at advanced ages.

In the Promislow et al. dataset, both the early increase in variance and the decline in variance for late ages under the logistic model are significant. Significant genetic variances are detected for all parameters of the logistic model. Evidently, original discrete analysis, aided by large subcohort sizes, showed the patterns of the current analysis quite well. We may gain some insight into the positive genetic correlation between the Gompertz slope parameter and from this agreement. Because in the logistic model the mortality curves are asymptotic to A/, our analysis is expected to detect a decline in variance at late ages if A and are highly correlated across genotypes. The presence of the decline in the discrete analysis of Promislow et al. suggests that its occurrence in the present analysis is a real biological feature of aging, as is the asymptotic plateau of mortality rates. These features of mortality trajectories are subject to numerous explanations.

As noted by Promislow et al., two technical causes for the declining genetic variance can be imagined. First, in the Promislow et al. study the cage density declined with cohort age. KHAZAELI et al. 1996 experimentally ruled out density effects upon overall mortality deceleration, yet interaction of density and genotype has not been explored and may yet contribute to the genetic variance decline. In Hughes and Charlesworth, density was controlled by adding marked replacement males, though these replacements were generally younger than those remaining. Second, environmental changes through the course of the Promislow et al. experiment could have produced demographic period effects, but independent observations suggest that is not a likely explanation. In numerous trials with cohorts of various genotypes, mortality rates in D. melanogaster are consistently observed to plateau at a level between 0.25 and 0.35 (PLETCHER et al. 1998 ). At the same time, mortality rates at earlier ages are observed to vary among such cohorts. In these trials the early-age genetic variance for mortality must degenerate because little phenotypic variance for the trait is observed at older ages. Thus, the result of Promislow et al. is likely to be robust across experimental periods, but this should be addressed explicitly in the future.

The constancy of the mortality asymptote is also a feature of both male and female lifetables when sexes are held separately. This is relevant to the possibility that the genetic variance for mortality declines in the mixed-sex cages of Promislow et al. because of changes in reproductive behavior and its associated mortality costs (PARTRIDGE 1987 ; CHAPMAN et al. 1995 ; TATAR and CAREY 1995 ). In the Promislow et al. dataset the decline in mortality genetic variance coincides with the final episodes of egg laying (TATAR et al. 1996 ). If variation among genotypes in mortality is due to variation in reproduction, then at late ages, when few genotypes are reproductively active, variation in mortality may decline. The male-only cohorts of Hughes and Charlesworth, on the other hand, experienced considerably less reproductive stress and would not be expected, under this explanation, to present a change in mortality variance at the end of reproductive life. However, as noted, the convergence of mortality plateaus to moderate levels is not a unique feature of mixed-sex cages of the sort used in Promislow et al. Therefore, the late-age decline in mortality variance may well be a feature of large-scale demographic studies where reproductive costs are held at a minimum. Furthermore, in the study of Promislow et al., per capita fecundity decreases from the first week onward (TATAR et al. 1996 ); this spans the period of marked increase of mortality variance. It seems unlikely that variance in reproductive activity could account for this complex pattern of mortality variance. As above, this explanation could be tested by repeating the Promislow et al. design with single-sex cages.

Two further biological explanations for the decline in genetic variance for mortality are related to the biodemographic problem of why mortality rates level off. Mortality deceleration was first recognized as a potential outcome of compositional heterogeneity within cohorts (VAUPEL and YASHIN 1985 ; VAUPEL 1990 ). When individuals within a cohort vary in frailty (the intrinsic baseline multiplier of age-dependent mortality risk, ), the sum trajectory of mortality across the cohort can decelerate even when the mortality risk of each individual increases exponentially with age. The total rate decreases less rapidly with age as the higher risk individuals drop out. In the logistic-frailty model of VAUPEL and YASHIN 1985 , describes the frailty distribution and thereby affects the pattern of mortality deceleration. Conversely, when a mortality trajectory is empirically fitted by the logistic model, the estimate of may describe otherwise unobservable heterogeneity for frailty. Of course, just because we observe mortality deceleration, it does not mean there is necessarily heterogeneity among individuals in frailty. It could also be due to intrinsic slowing of the aging process within individual flies.

Compositional heterogeneity for frailty could play a strong role in the decline of mortality variance observed by Promislow et al. In these data, both and V_A() are significantly different from zero. If mortality deceleration is caused by frailty heterogeneity, these observations imply that the distribution of individual frailty is correlated genetically among subcohorts; a subcohort of a particular genotype at the second chromosome has not only a particular Gompertz slope and intercept but a particular variance in individual frailty as well. This might be explained by interactions of the second chromosome genes, which are uniform in each subcohort, with heterogeneous aspects of their environment. Genetic variance for results if sensitivity to environmental variation varies among the second chromosome genotypes. These conditions are plausible for D. melanogaster. Additive variance for reaction norms has been documented for chromosomal isolates with respect to competitive ability (GURGANUS et al. 1998 ). Furthermore, both specific and general sources of variation were present in Promislow et al.

In the Promislow et al. design, each second chromosome line was associated with an uncontrolled array of balancer stock genes originating from the creation of the line. Thus, a potentially different array of third and X chromosome genotypes was associated with each second chromosome genotype. This could affect both and V_A() in two ways. The distribution of frailty within a cohort may be due to the chance inheritance of various first and third chromosomes among individuals of the group. In this case, frailty would covary with second chromosomes, but not in a meaningful manner. On the other hand, the arrays of background chromosomes may specifically interact with the genotypes at the second chromosome. With gene-by-background interaction, covariance of second chromosomes with frailty heterogeneity will exist when the extracted chromosomes vary in their norms of reaction.

In contrast to Promislow et al., Hughes and Charlesworth placed each third chromosome (which was derived from the Ives population) on two independent Ives genetic backgrounds to avoid confounding third chromosome effects with effects due to other chromosomes and to eliminate the possibility of hybrid dysgenesis. Subcohort values of might therefore be lower than those in the Promislow et al. experiment and genetically uncorrelated with each other. The mean value of in the Hughes and Charlesworth results is an order of magnitude less than that in the Promislow et al. data for both males and females. However, the observed value of in Promislow et al. is typical of values estimated for highly isogenic, inbred lines studied in the same laboratory. In this context, the absence of detectable V_A() in the Hughes and Charlesworth data, otherwise attributable to subcohort sample size, may actually be due to the highly isogenic background. Clearly, it will be useful to combine the strengths of these studies to estimate V_A() with both large-scale cohorts and relative homogeneity of background.

The primary alternative to the hypothesis that the observed bending of mortality curves results from heterogeneity of frailty is that the rate of senescence slows with age at the level of individual experience. In this case, the empirically fit parameter coincidentally reflects an internal epigenetic process. This is the least understood explanation for mortality deceleration and genetic variance decline. Promislow et al. provisionally suggested that these patterns could result "if late-acting deleterious mutations are expressed at all ages subsequent to the age of onset" (PROMISLOW et al. 1995 , p. 846). Mutations may affect the probability or the timing of state changes in regulatory or signal transduction systems. When a limited number of such systems affect the performance of an integrated unit, reliability theory predicts that failure rates will reach an asymptote (GAVRILOV and GAVRILOVA 1991 ). Further, with a limited number of subsystems, the progressive accumulation of mutations will have redundant effects on the failure outcome. In this way there may be little genetic variance for late-life mortality. Consistent with this model, the first detailed examination of age-specific mutation effects upon mortality revealed mutational covariance among age classes and a striking lack of mutational variance for mortality at late ages (PLETCHER et al. 1998 ). In contrast, current assumptions of models for both mutation accumulation and antagonistic pleiotropy consider that mutational and allelic effects are expressed at only a single age. The effect of among-age covariance and age-specific mutational variance upon the evolution of mortality and its standing genetic variance has yet to be incorporated in either evolutionary model of senescence.

	ACKNOWLEDGMENTS

F.H.S. thanks S. C. Stearns, T. Kawecki, J. W. Curtsinger, S. D. Pletcher, R. Shaw, and P. M. Service for many helpful discussions of this work. The experimental work of D.E.L.P. and M.T. was supported in part by the American Federation of Aging Research.

Manuscript received November 12, 1998; Accepted for publication March 1, 1999.

	APPENDIX 1

TOP ABSTRACT Empirical studies of variance... Statistical challenges for the... THE MODEL AND STATISTICAL... RESULTS DISCUSSION APPENDIX 1 LITERATURE CITED

The Monte Carlo likelihood-ratio function requires the calculation of

which would be very time consuming if the matrix V, of order 300 in the case of the Promislow et al. data and of order 912 in the Hughes and Charlesworth data, were not easy to invert for a given point () in the 12-dimensional parameter subspace it occupies. The following method was applied to accomplish the inversion and multiplication (THOMPSON 1976 ).

The phenotypic covariance matrix V is the sum of the additive genetic and residual covariance matrices

where denotes an outer product, A is the additive genetic relationship matrix, I is the identity matrix, and _A and _R are the additive genetic and residual covariance component matrices, respectively. If T is the matrix whose rows are eigenvectors of corresponding to vector of eigenvalues , we have TAT' = diag(_i). In the case of the data at hand, most of the eigenvalues are zero and there are only a few distinct eigenvalues of varying multiplicity in all. Transforming y_k to Ty_k, we have

which, given that TIT' = I, is a block diagonal matrix where the blocks (one for each A matrix eigenvalue) are of order 3. Conveniently, with new values of , the eigenvalues and eigenvectors do not change, and inverting this matrix is an easy matter of inverting k 3 x 3 matrices

where k is the number of distinct eigenvalues _i in A (three in the case of the Promislow et al. data and three to eight depending on the block in the Hughes and Charlesworth data).

The product

becomes

where T_i,l is the lth eigenvector (of n_i) and V^-1_i is the block of V^-1, corresponding to the ith distinct eigenvalue.

It is noted that the structure of V as the sum of constant multiples of A and I is crucial in the derivation of the above. If dominance is included in the genetic model and

the eigenvalues and eigenvectors of V change unpredictably with each value of so that a full inversion is required at each point in the parameter space.

	LITER