
==== Front
Genetics
Genetics
genetics
Genetics
0016-6731
1943-2631
Oxford University Press US

37450606
10.1093/genetics/iyad133
iyad133
Investigation
Population and Evolutionary Genetics
AcademicSubjects/SCI01180
AcademicSubjects/SCI01140
Featured
The infinitesimal model with dominance
https://orcid.org/0000-0002-8548-5240
Barton Nicholas H Institute of Science and Technology, Am Campus I, A-3400 Klosterneuberg, Austria

Etheridge Alison M Department of Statistics, University of Oxford, 24–29 St Giles, OX1 3LB Oxford, UK

Véber Amandine MAP5, Université Paris Cité, CNRS, 45 rue des Saints-Pères, 75006 Paris, France

Martin G Editor
Corresponding author: MAP5, Université Paris Cité, 45 rue des Saints-Pères, 75006 Paris, France. Email: amandine.veber@parisdescartes.fr
Conflicts of interest The author(s) declare no conflict of interest.

10 2023
14 7 2023
14 7 2023
225 2 iyad13302 11 2022
23 6 2023
28 8 2023
© The Author(s) 2023. Published by Oxford University Press on behalf of The Genetics Society of America.
2023
https://creativecommons.org/licenses/by/4.0/ This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

The classical infinitesimal model is a simple and robust model for the inheritance of quantitative traits. In this model, a quantitative trait is expressed as the sum of a genetic and an environmental component, and the genetic component of offspring traits within a family follows a normal distribution around the average of the parents’ trait values, and has a variance that is independent of the parental traits. In previous work, we showed that when trait values are determined by the sum of a large number of additive Mendelian factors, each of small effect, one can justify the infinitesimal model as a limit of Mendelian inheritance. In this paper, we show that this result extends to include dominance. We define the model in terms of classical quantities of quantitative genetics, before justifying it as a limit of Mendelian inheritance as the number, M, of underlying loci tends to infinity. As in the additive case, the multivariate normal distribution of trait values across the pedigree can be expressed in terms of variance components in an ancestral population and probabilities of identity by descent determined by the pedigree. Now, with just first-order dominance effects, we require two-, three-, and four-way identities. We also show that, even if we condition on parental trait values, the “shared” and “residual” components of trait values within each family will be asymptotically normally distributed as the number of loci tends to infinity, with an error of order 1/M. We illustrate our results with some numerical examples.

The classical infinitesimal model is a simple and robust model for the inheritance of quantitative traits. In previous work, Barton, Etheridge, and Véber showed that—when trait values are determined by the sum of a large number of Mendelian factors, each of small effect—one can justify the infinitesimal model as limit of Mendelian inheritance. Here, the authors extend these results to include a nonlinear dominance effect.

infinitesimal model
dominance
ERC 10.13039/100010663 250152 101055327
==== Body
pmcIntroduction

In the classical infinitesimal model, a quantitative trait is expressed as the sum of a genetic and a nongenetic (environmental) component, and the genetic component of offspring traits within a family follows a normal distribution around the average of the parents’ trait values, and has a variance that is independent of the trait values of the parents. With inbreeding, the variance decreases in proportion to relatedness. When trait values are determined by the sum of a large number of Mendelian factors, each of small effect, as we show in Barton et al. (2017), one can justify the infinitesimal model as a limit of Mendelian inheritance. Crucially, the results of Barton et al. (2017) show that the evolutionary forces such as random drift and population structure are captured by the pedigree; conditioning on that pedigree, and trait values in the population in all generations before the present, the within-family distributions in the present generation will be given by a multivariate normal, with variance determined by that in the ancestral population and probabilities of identity by descent that can be deduced from the pedigree. If some traits in the pedigree are unknown, then averaging with respect to the ancestral distribution, the multivariate normality is preserved. It was also shown that under some forms of epistasis, trait values within a family are still normally distributed, although the mean will no longer be a simple function of the traits in the parents (as there are epistatic components which cannot be observed directly).

We emphasize that as a result of selection, population structure, and so on, the trait distribution across the population can be far from normal; the infinitesimal model as we define it only asserts that the within-family distributions of the genetic component of the trait are Gaussian, with a variance–covariance matrix that is determined entirely by that in an ancestral population and the probabilities of identity determined by the pedigree. Moreover, as a result of the multivariate normality, conditioning on some of the trait values within that pedigree has predictable effects on the mean and variance within and between families. In other words, knowing the trait values for some individuals in the population does not distort the multivariate normality of the distribution of the unobserved traits, and the mean and covariances of these traits may be derived explicitly (albeit after rather tedious calculations).

In this paper, we show that this extraordinary robustness of the infinitesimal model extends to include dominance. The distribution of the genetic part of the trait will once again be a multivariate normal distribution whose mean and variance is expressed in terms of the variance components in an ancestral population and probabilities of identity by descent determined by the pedigree, but now, with just first-order dominance effects, the identities required will involve up to four genes. As with the case of epistasis, the mean is not a simple function of the trait values in the parents, and there is nontrivial covariance between families. One can think of the genetic component of the trait values within a family as consisting of two parts. Both are normally distributed. In the additive case, the first reduces to the mean of the trait values of the parents; with dominance it will be random (even if we condition on knowing the parental traits), but the same for all individuals in the family. What is at first sight surprising is that even if we condition on knowing the trait values of the parents, this shared quantity is normally distributed. Assuming there is no mutation to ease the presentation [the effect of mutation was studied in Barton et al. (2017)], our first contribution is to show how to calculate its mean and variance from knowledge of variance components in the ancestral population and the pedigree, both with and without knowledge of the trait values of the parents. Knowing the trait values of the parents shifts the mean in a predictable way, the variance is independent of the parental trait values. The second part of the trait value, which is independent for each offspring in the family, is independent of the first; it encodes the randomness of Mendelian inheritance. It is a draw from a normally distributed random variable with mean zero and variance again determined by the pedigree and variance components in the ancestral population. It is not affected by conditioning on parental trait values. This segregation of the trait into a shared part and a residual part that is independent for each member of a family is not the classical subdivision into additive and dominance components, but it arises naturally both in the formulation of the infinitesimal model and in its derivation as a limit of Mendelian inheritance for a large number of loci each of small effect. We give a more mathematical description of it in Equation (2).

Our work can be seen as an extension of that of Abney et al. (2000), who establish sufficient conditions for a Central Limit Theorem to be applied to the vector of trait values in the presence of dominance and inbreeding. Our second contribution in this work is to establish the magnitude of the error in that normal approximation, verify that in conditioning on the trait values of the parents of an individual we are not (unless those traits are very extreme or the pedigree is very inbred) leaving the domain where the normal approximation is valid, and write down the effect of knowing those parental trait values on the distribution of the individual’s own trait. A careful statement of our results can be found in Theorems 4 and 5. The notation we shall need is rather involved, but in a nutshell, we shall write the trait Z~i of a given diploid individual i in generation t as the sum over M loci of per-locus allelic effects that are functions of the allelic states χl1,χl2 of the two genes of individual i at locus l, plus an environmental contribution Ei (that we shall assume to be Gaussian):

(1) Z~i=z¯0+∑l=1M1M(ηl(χl1)+ηl(χl2)+ϕl(χl1,χl2))+Ei.

Here, z¯0 is the average trait value in the ancestral population (itself a sum of average allelic effects) and the sum encodes the contribution of all loci to the deviation from this average [each per-locus deviation being of order 1/M, see Barton et al. (2017) and the third section below for a justification]. In this sum, the term ηl(χl1)+ηl(χl2) models the additive part of the contribution of locus l and ϕl(χl1,χl2) models the part due to dominance. Assuming Mendelian inheritance and no linkage between the M loci, at each locus the allelic state χl1 is a copy of the allelic state of one of the two genes in the “first” parent of i, chosen at random, and χl2 is a copy of the allelic state of one of the two genes in the “second” parent of i, again chosen uniformly at random. Writing χli[1],1,χli[1],2 for the alleles at locus l in the first parent and χli[2],1,χli[2],2 for the alleles in the second parent, we can then write the sum over all loci in Equation (1) as the sum of an average parental contribution (shared by all offspring of these parents), and a residual term of mean zero that encodes the stochasticity of Mendelian inheritance (the actual genetic contribution of the parents minus their average contribution). To avoid introducing even more notation, here we simply write RAi and RDi for the parts of the residual due to the additive terms and to the dominance terms respectively. Explicit formulae are given in Equations (24)–(27). Doing so, we obtain

(2) Z~i=z¯0+1M∑l=1M{ηl(χli[1],1)+ηl(χli[1],2)2+ηl(χli[2],1)+ηl(χli[2],2)2}+1M∑l=1M{ϕl(χli[1],1,χli[2],1)+ϕl(χli[1],2,χli[2],1)4+ϕl(χli[1],1,χli[2],2)+ϕl(χli[1],2,χli[2],2)4}+RAi+RDi+Ei=:z¯0+Ai+Di+RAi+RDi+Ei.

The genetic component of the trait can thus be seen either as the sum of an additive part (Ai+RAi) and a dominance part (Di+RDi), or as the sum of a shared part (Ai+Di) and a residual part (RAi+RDi). Following the same strategy as in Barton et al. (2017), in Theorem 4 we show that even conditionally on (i.e. knowing) the parental traits Z~i[1] and Z~i[2], as M tends to infinity the residual part converges in distribution to a Gaussian distribution with mean 0 and a variance depending only on variance components in the ancestral population and on the probability of identity by descent between two parental genes (which is fully determined by the pedigree). Crucially, the limiting variance does not depend on the parental traits. This convergence happens at a rate proportional to 1/M. Turning to the shared part, we use a different approach to prove that conditional on Z~i[1] and Z~i[2], Ai+Di also converges to a Gaussian distribution as M tends to infinity. Again, the nonzero mean and the variance of the limiting normal distribution can be fully described, the variance is independent of the parental traits and the convergence happens at a rate proportional to 1/M. This is the content of Theorem 5, in the special (and most difficult) case when individual i was produced by selfing. For both the shared and the residual parts, the rate of convergence deteriorates when the pedigree is too inbred (leading to probabilities of identity by descent close to 1 between some pairs of parental genes), or when some traits in the population are too extreme (as knowing the trait value then gives us too much information about the unobserved underlying allelic states).

Our derivation of the infinitesimal model as the limit of a finite-locus model has two interesting corollaries. First, as mentioned above, we obtain that the error made by approximating the trait distribution within a family by a Gaussian distribution increases by a quantity of order 1/M in each generation. Consequently, for very large M, we expect the infinitesimal model with dominance to be valid for a time of the order of M generations, provided the population is not too inbred and no too extreme traits appear in the meantime. Second, the set of technical lemmas that are key to the proofs of these results, presented in Appendix E, show that the infinitesimal model leaves essentially no signature on the allele frequencies at any given locus: even knowing the ancestral traits, the distribution of the allelic state at a single locus in a given individual is barely distorted by selection acting on the trait and the result is that, at the population level, the allelic distribution evolves in an essentially neutral way. In particular, its variance only depends on the variance of the allele distribution in the ancestral population and on identities by descent, that are not changed by knowledge of the trait values.

The rest of the paper is organized as follows. In the next section, we define the identity coefficients (that is, probabilities of identity by descent) that we shall need to formulate the model precisely. We show how to compute them knowing the population pedigree in Appendix A and provide the corresponding Mathematica code in supplementary material (Barton 2023). Next, we spell out the model in terms of quantities that are familiar from classical quantitative genetics, and we explore its accuracy numerically in a devoted section. Finally, we derive this extension of the infinitesimal model as a limit of a model of Mendelian inheritance on the pedigree. The calculations are somewhat involved, and almost all will be relegated to the appendices. We must modify the strategy of Barton et al. (2017), which, although valid for the part of the trait value which is independent for each individual within the family, does not suffice for proving normality of the part of the trait value that is shared by all individuals within a family. To prove that this is normally distributed requires a new approach, based on an extension of Stein’s method of exchangeable pairs. To keep the expressions in our calculations manageable, we satisfy ourselves with presenting the details only in the case in which we condition on knowing the trait values of the parents of an individual, in contrast to the additive case of Barton et al. (2017), in which we conditioned on knowing all the trait values in the pedigree right back to the ancestral generation. Our approach could readily be extended to conditioning on knowledge of more trait values, which amounts to conditioning a multivariate normal on some of its marginals. In Appendix H, we present the new ideas that are required to control the way in which errors in the infinitesimal approximation accumulate from knowledge of trait values of more distant relatives in the presence of dominance.

Just as in the additive case, the key will be to show that because many different combinations of allelic states are consistent with the same trait value, knowledge of the pedigree, and the trait values of the parents of an individual in that pedigree, actually gives very little information about the allelic state at a particular locus in that individual, or about correlations between two specific loci. An important consequence of this is that, in practice, it is going to be hard to observe signals of polygenic adaptation, because even a large shift in a trait caused by strong selection does not yield a prediction about alleles at a particular locus.

Identity coefficients

In the case of an additive trait, the infinitesimal model can be expressed in terms of the variance in the ancestral population (that is, the base population which we shall call generation zero) and two-way identity coefficients at a single locus. Recall that two genes at a given locus are identical by descent if their allelic states are identical and were inherited from a common ancestor. Since we assume that individuals are diploid, we need to specify which genes we consider when defining the identity coefficients.

For two distinct individuals i and j in the same generation, we define Fij to be the probability of identity by descent between two genes (at a given locus), one taken uniformly at random among the two genes of individual i and one taken at random among the two genes of individual j. When i=j, Fii is defined to be the probability of identity by descent of the two distinct genes in the diploid individual i.

The definition naturally extends to subsets of three or four genes taken from two distinct individuals (again, at a given locus), for which we shall talk about three- and four-way identities. These quantities will be required to state our results below.

We use F122 for the probability that the two genes in individual 2 are identical by descent and they are identical by descent with a gene chosen at random from individual 1. We write F1122 for the probability that all four genes across individuals 1 and 2 are identical by descent; this corresponds to the quantity δ in Walsh and Lynch (2018, Chapter 11). We need an expression for the probability that each gene in individual 1 is identical by descent with a different gene in individual 2 and all four are not identical. We shall denote this by F~1212. This is denoted by (Δ−δ) in Walsh and Lynch (2018). Finally, we need the probability that the two genes in individual 1 are identical, as are the two genes in individual 2, but the four genes are not all identical, which we shall denote by F~1122. We illustrate the three- and four-way identities in Fig. 1. During the course of our mathematical derivations, it will be convenient to express all two-, three-, and four-way identities in terms of the nine possible four-way identities (Walsh and Lynch 2018, Fig. 11.5). This is illustrated in Fig. B1.

Fig. 1. Three- and four-way identities. Lines indicate identity by descent between genes. See the main text for further explanation.

In Appendix A, we discuss how to compute these identity coefficients given a pedigree. From now on we simply write “identity” instead of “identity by descent.”

The infinitesimal model with dominance

For ease of exposition, in this section we leave aside the environmental component of the trait value and we focus on its genetic component, which we denote by Z [so that in the notation of Equation (1), Z~=Z+E]. We first introduce the different quantities that are involved in this component of the trait value in a rigorous way, most of which were already hinted at in the Introduction, and then we compute the mean and variance of the shared and residual parts of Z with and without knowledge of the parental traits.

The population is diploid and trait values are determined by the allelic states at M unlinked loci. Each locus thus corresponds to a pair of genes. We assume that in generation zero (i.e. in the “ancestral” population), the individuals that found the pedigree are unrelated and sampled from an ancestral population in which all loci are in linkage equilibrium and are in Hardy–Weinberg equilibrium (that is, in the ancestral population the two allelic states at each locus in a given individual are sampled independently of each other and therefore the probability that an individual carries a given pair of alleles is given by the product of the probabilities of each allele being sampled).

In order to define the various quantities that enter into our model, we introduce notation to express the trait as a sum of effects over loci. However, we emphasize that once these components, all of which are familiar from classical quantitative genetics, have been calculated for the ancestral population, the model can be defined without reference to the effects of individual loci.

To adhere to the notation of Barton et al. (2017), we use χl1, χl2 for the allelic states of the two genes at locus l in a given individual in the pedigree. When we talk about the distribution of the allelic state of a single gene, we drop the superscript 1 or 2 and simply write χl. We write z¯0 for the mean trait value in the ancestral population and express the trait value of an individual as z¯0 plus a sum of allelic effects. The influence of each locus will scale as 1/M, where M is the total number of loci (assumed large). We write ηl(χl) to denote the (order one) scaled additive effect of the allele χl and ϕl(χl1,χl2) for the scaled dominance component (where ϕl is assumed to be a symmetric function of the two allelic states χl1 and χl2). That is, the total contribution of locus l to the trait value will be of the form

1M(ηl(χl1)+ηl(χl2))+1Mϕl(χl1,χl2).

We shall assume that both ηl and ϕl are uniformly bounded (i.e. they will all take their values in some finite interval [−B,B]). We also suppose that dominance effects are sufficiently “balanced” that inbreeding depression is finite at least in the ancestral population. More precisely, let χ^l denote an allele sampled at random from the distribution of alleles at locus l in the ancestral population, then ι defined by

(3) ι=1M∑l=1ME[ϕl(χ^l,χ^l)]

is bounded (as a function of M). This condition is crucial to our result. It is not obvious that it can hold, as the number of terms in the sum grows linearly with M while the scaling factor 1/M decreases much more slowly. Such a uniform bound is possible for instance if we consider a situation in which the contributions of the different loci compensate each other in a “random-walk-like” way, i.e. each expectation is either positive or negative (by the same amount, say), and the number of positive and negative expectations differ by at most O(M). An example is presented at the beginning of the section on numerics. Note however that the quantity ι may be bounded uniformly in M for many other reasons. For simplicity, we do not consider higher order dominance components (that is D×D—or more complex—components) here.

Remark 1 Note that χ^l is the random variable describing a draw from the distribution of allelic states at locus l in the ancestral population (generation 0), while we use χl to denote the allelic state at locus l in a given individual in the pedigree (living in generation t, say). A priori, the law of χl is a biased version of the law of χ^l, obtained after letting selection and drift act over t generations, but in Appendix E we shall show that, in effect, this distortion is very small for each given locus, and χ^l and χl have the same distribution up to a small error even if we condition on knowing the parental (or ancestral) trait values.

For an individual in the ancestral population, its allelic states at locus l, which we denote by χ^l1,χ^l2, are independent draws from a distribution ν^l on possible allelic states that we assume is known. It is convenient to normalize so that E[ηl(χ^l)]=0, E[ϕl(χ^l1,χ^l2)]=0, and for any value x′ of the allelic state at locus l, the conditional expectations E[ϕl(χ^l,x′)]=0=E[ϕl(x′,χ^l)]. We explain in the section on modeling Mendelian inheritance why these assumptions do not result in a loss of generality. The genetic component of the trait value takes the form [compare with Equation (1), the expression for the observed trait including environmental noise]

(4) Z=z¯0+1M∑l=1M(ηl(χl1)+ηl(χl2)+ϕl(χl1,χl2)).

Let us write i[1] and i[2] for the parents of the individual labeled i. As advertised in the Introduction, the genetic component of an offspring’s trait value has two contributions. The first one is shared by all its siblings, and is a random quantity which is characteristic of the family. The second contribution is unique to the individual and independent of the first one. In our proofs, we shall investigate these two parts separately. We shall use the notation Zi=(Ai+Di)+(RAi+RDi), where the shared part has been further subdivided into the contribution Ai from the additive component, and the contribution Di from the dominance component. The residuals RAi and RDi are determined by Mendelian inheritance and correspond to the contributions from the additive and dominance components respectively. Explicit expressions for these quantities are in Equations (24)–(29) below. In this notation, the additive part of the trait value is Ai+RAi and the dominance deviation is Di+RDi.

Trait values for a given pedigree

We now define the infinitesimal model in terms of classical quantities of quantitative genetics that can be expressed in terms of expectations in the ancestral population and identities determined by the pedigree. We use the notation of Walsh and Lynch (2018), which we recall in Table 1. Under the infinitesimal model, conditional on the pedigree, the components (Ai+Di) and (RAi+RDi) of the trait values of individuals in a family follow independent multivariate normal distributions. In Appendix B, the expressions presented in this section will be justified by taking the trait values determined by Equation (4) under a model of Mendelian inheritance. In writing down the infinitesimal model, we shall assume that as the number of loci tends to infinity, the quantities defined in the top part of Table 1 converge to well-defined limits.

Table 1. Coefficients of classical quantitative genetics (top) and elements of individual trait decomposition (bottom).

Additive variance	σA2=2M∑l=1ME[ηl(χ^l)2]	
Dominance variance	σD2=1M∑l=1ME[ϕl(χ^l1,χ^l2)2]	
Inbreeding depression	ι=1M∑l=1ME[ϕl(χ^l,χ^l)]	
Sum of squared locus-specific inbreeding depressions	ι*=1M∑l=1ME[ϕl(χ^l,χ^l)]2	
Variance of dominance effects in inbred individuals	σDI2=1M∑l=1M(E[ϕl(χ^l,χ^l)2]−E[ϕl(χ^l,χ^l)]2)	
Covariance of additive and dominance effects in inbred individuals	σADI=2M∑l=1ME[ηl(χ^l)ϕl(χ^l,χ^l)]	
Additive part of the shared component	Ai —defined in Equation (28)	
Dominance part of the shared component	Di —defined in Equation (29)	
Additive part of the residual	RAi —defined by Equations (24) + (25)	
Dominance part of the residual	RDi —defined by Equations (26) + (27)	
Genetic component of trait value	Zi=z¯0+Ai+Di+RAi+RDi	
Observed trait value	Z~i=Zi+Ei , Ei∼N(0,σE2)	
We use χ^l to denote an allelic state sampled from the distribution ν^l of possible allelic states at locus l in the ancestral population; χ^l1, χ^l2 are independent draws from the same distribution.

To simplify notation, we shall use 1 and 2 in place of i[1] and i[2] in our expressions for identity; thus, for example, F12≡Fi[1],i[2], and F11 will be the probability of identity by descent of the two genes in parent i[1]. The mean and variance of (Ai+Di) are then

(5) E[Ai+Di]=ιF12,

and

(6) Var(Ai+Di)=σA22(1+F11+F222+2F12)+σADI(F12+F112+F1222)+(σDI2+ι*)4(F12+F112+F122+F1122)+ι*4F~1212−ι*F122+σD24(1−F12+F22−F122+F11−F112+F~1122+12F~1212).

In this expression, the term proportional to σA2 is the variance of Ai, the term proportional to σADI is twice the covariance of Ai and Di and the remaining sum gives the variance of Di. Recall that we are assuming here that the ancestral population is in linkage equilibrium. With linkage disequilibrium there is an additional term, c.f. the remark below Equation (11). The components (A+D) are also correlated across families. For individuals labeled i and j, respectively,

(7) Cov((Ai+Di),(Aj+Dj))=2FijσA2+(Fijj+Fiij)σADI+F~ijijσD2+Fiijj(σDI2+ι*)−ι2FiiFjj+ι*F~iijj.

Note that, in contrast to our expression for the variance of Zi, in this expression, the subscripts i and j in the identities refer to the individuals themselves, not their parents; for example, the expression Fij is the probability of identity of two genes, one sampled at random from individual i and one sampled at random from individual j. We reserve letters for individuals in the current generation, and numbers for their parents.

If we combine the components RAi and RDi that segregate within families, we have that the sums (RAi+RDi) are independent of each other (due to the independence of the variables encoding Mendelian inheritance), mean zero, normally distributed random variables with variance

(8) Var(RAi+RDi)=(1−F11+F222)σA22+14(3F12−F1122−F112−F122)(σDI2+ι*)+14[3(1−F12)−(F11−F112)−(F22−F122)−F~112212−12F~1212]σD2+(F12−F112+F1222)σADI−ι*4F~1212.

Here again, the term proportional to σA2 is the variance of RAi, the term proportional to σADI is twice the covariance of RAi and RDi, and the remaining sum equals the variance of RDi. We calculate the mean, variance, and covariance of these different components in Appendix B. In order to recover the mean and variance of the trait values, we add the contributions of (Ai+Di) and (RAi+RDi) and observe that the identity F12 in our expressions for the variances of these quantities (which we recall was the probability of identity of one gene sampled at random from each of the parents i[1], i[2] of our individual) corresponds to Fii. This yields that, conditional on the pedigree,

(9) E[Zi]=z¯0+ιFii,

(10) Cov(Zi,Zj)=2FijσA2+(Fijj+Fiij)σADI+F~ijijσD2+Fiijj(σDI2+ι*)−ι2FiiFjj+ι*F~iijj,

and

(11) Var(Zi)=σA2(1+Fii)+σD2(1−Fii)+(σDI2+ι*)Fii+2σADIFii−ι*Fii2.

For a single individual, its trait value can only depend on the two alleles that it carries at each locus, so it is no surprise that this expression depends only on pairwise identities between those two genes. We remark that Equation (11) differs from the corresponding expression [Equation (11.6c) in Walsh and Lynch (2018)]. To recover exactly their expression, one must add (f~−Fii2)(ι2−ι*) to the right-hand side, where f~ is the probability of identity at two distinct loci in individual i. We see how to recover this term in Remark B1, but because we have assumed linkage equilibrium in our base population, for the period over which the infinitesimal model remains a good approximation, under our assumptions we have f~≈Fii2. This is not to say that there is not a significant contribution to the trait value from linkage disequilibrium; it is just that for any specific pair of loci it is negligible. We shall see a toy example that reinforces this point at the beginning of the section on modeling Mendelian inheritance.

We emphasize again that our partition of the trait values into a contribution that is shared by all individuals in a family and residuals differs from the conventional split into an additive part and a dominance deviation. The additive part of the trait is Ai=Ai+RAi and the dominance component is Di=Di+RDi. From our calculations in Appendix B, we can read off

(12) E[Ai]=0,E[Di]=ιFii,

(13) Var(Ai)=σA2(1+Fii),Cov(Ai,Di)=σADIFii,

and

(14) Var(Di)=σD2(1−Fii)+σDI2Fii+ι*(Fii−Fii2).

Remark 2 Notice that the purely additive case can be simply recovered by taking ϕl≡0, so that Di=0=RDi, and σA2 is the only nonzero variance coefficient. This yieldsE[Ai+Di]=0,Var(Ai+Di)=σA22(1+F11+F222+2F12),Cov((Ai+Di),(Aj+Dj))=2FijσA2,Var(RAi+RDi)=(1−F11+F222)σA22,

and finallyE[Zi]=z¯0,Var(Zi)=σA2(1+Fii),Cov(Zi,Zj)=2FijσA2.

Conditioning on trait values of parents

Under the infinitesimal model, the trait values of individuals across the pedigree are given by a multivariate normal. Therefore, standard results on conditioning multivariate normal random vectors on their marginal values, which for ease of reference we record in Appendix C, allow us to read off the effect on the distribution of Zi of conditioning on Zi[1] and Zi[2]. However, a little care is needed; we shall be justifying the normal distribution within families as an approximation as the number of loci tends to infinity, and we must be sure that asymptotic normality is preserved under this conditioning. We shall see that if, for example, parental trait values are too extreme, then the conditioning pushes us to a part of the probability space where the normal approximation breaks down. This is particularly evident in the toy example that we present in the section on modeling Mendelian inheritance. A justification for asymptotic normality even after conditioning is outlined in that section, and details are presented in the appendices.

Just as in the classical infinitesimal model, the mean and variance of the residuals RAi+RDi are unchanged by conditioning on the trait values of the parents [recall that these residuals encode the stochasticity due to Mendelian inheritance at each locus; expressions for RAi and RDi are given in Equations (24)–(27)]. For the shared components, the mean and variance will be distorted by quantities determined by the covariances between (Ai+Di) and Zi[1], Zi[2]. Let us write

(15) C(i,i[1]):=Cov((Ai+Di),Zi[1]),

with a corresponding definition for C(i,i[2]). Then, once again using 1 and 2 in place of i[1] and i[2] in our expressions for identities,

(16) C(i,i[1])=σA22(1+F11+2F12)+σADI2(F11+F12+2F112)+σD2(F12−F112)+(σDI2+ι*)F112−ι2F11F12,

with C(i,i[2]) given by the corresponding expression with the roles of the subscripts 1 and 2 interchanged. (A derivation of this expression is provided in Appendix B.) With this notation,

(17) E[(Ai+Di)|Zi[1],Zi[2]]=E[(Ai+Di)]+1Var(Zi[1])Var(Zi[2])−Cov(Zi[1],Zi[2])2×{(C(i,i[1])Var(Zi[2])−C(i,i[2])Cov(Zi[1],Zi[2]))×(Zi[1]−E[Zi[1]])+(C(i,i[2])Var(Zi[1])−C(i,i[1])Cov(Zi[1],Zi[2]))×(Zi[2]−E[Zi[2]])},

and

(18) Var((Ai+Di)|Zi[1],Zi[2])=Var(Ai+Di)−Var(Zi[1])C(i,i[2])2+Var(Zi[2])C(i,i[1])2Var(Zi[1])Var(Zi[2])−Cov(Zi[1],Zi[2])2+2Cov(Zi[1],Zi[2])C(i,i[1])C(i,i[2])Var(Zi[1])Var(Zi[2])−Cov(Zi[1],Zi[2])2.

(We have implicitly assumed that i[1]≠i[2]; in the case i[1]=i[2] the expression is simpler as we are then conditioning a bivariate normal on one of its marginals.)

Remark 3 In the purely additive case, things simplify greatly. From the expressions above, before conditioning, the mean of Ai+Di is zero (since ι=0), and the variance isσA22(1+(F11+F22)2+2F12).

Moreover,Var(Zi[1])=σA2(1+F11),Var(Zi[2])=σA2(1+F22),Cov(Zi[1],Zi[2])=2σA2F12,

andC(i,i[1])=12σA2(1+F11+2F12),C(i,i[2])=12σA2(1+F22+2F12).

Substituting into Equations (17) and (18), and observing that(1+F11)(1+F22+2F12)2+(1+F22)(1+F11+2F12)2−4F12(1+F11+2F12)(1+F22+2F12)=2((1+F11)(1+F22)−4F122)(1+F11+F222+2F12),

we find that conditional on the trait values of the parents, the mean and variance of Ai+Di reduce to (Zi[1]+Zi[2])/2 and zero, respectively, and we recover the classical infinitesimal model.

Although in the presence of dominance the expressions (17) and (18) are rather complicated, we emphasize that they are derived from knowledge of just the ancestral population and the pedigree, and are expressed in terms of familiar quantities from classical quantitative genetics.

Numerical examples

In this section, we present numerical examples to illustrate the accuracy of the predictions of the infinitesimal model, again disregarding the environmental component of the trait.

We first generated a pedigree for a population of constant size of N=30 diploid individuals over 50 discrete generations. Mating is random, but with no selfing. In order to facilitate comparison of different scenarios, the same pedigree was used for all subsequent simulations. In this way, the identity coefficients are held constant. As expected, the mean probability of identity between pairs of genes sampled from different individuals in generation t is close to 1−(1−1/2N)t.

We define a trait, Z, which depends on M=1,000 bi-allelic loci. There is no epistasis, so that the trait value is a sum across loci. In the examples here, we assume complete dominance, so that the effects of the three genotypes at each locus are either −α:−α:+α or −α:+α:+α. In order to ensure that the inbreeding depression ι is bounded, we need to have some “balance” and so we choose the effects at each locus according to an independent Bernoulli random variable with parameter H; that is, the probability that the effects across the three genotypes at locus l is −α:−α:+α is 1−H, independently for each locus. The effect size α is taken to be 1/M for all loci and H=12+2M. With these choices the additive and dominance variances will be O(1).

In the ancestral population, the allele frequencies were generated to mimic neutral allele frequencies with very low mutation rates, but conditioned to segregate at each locus. Thus, allele frequencies at every locus were sampled independently and according to a distribution with density proportional to (p(1−p))1−ϵ, with ϵ=0.001, but with those in [0,1/60] and [1−1/60,1] discarded (and the distribution renormalized). Then for each population replicate, these frequencies were used to endow each individual in the base population with an allelic type at every locus.

Variance components are defined with respect to this reference set of allele frequencies. For the population generated for the examples presented here, these values were σA2=0.269, σD2=0.073, and the inbreeding depression ι=−0.531. The additive and dominance components are uncorrelated in the base population (Cov(A,D)=0). In the numerical experiments that follow, each replicate population is started at time zero from a different collection of genotypes, sampled from this base distribution.

We first simulated a neutral model. Figure 2 illustrates how the different components of the trait values change over fifty generations of neutral evolution. Recall that we always use the same realization of the pedigree. For each replicate, we take an independent sample of allelic types at time zero. For each individual in the pedigree we evaluate the additive and dominance components A and D and then in each generation we calculate the mean and variance of these quantities across the 30 individuals in the population. This is only intended to give some feeling for the ways in which the components fluctuate through time. Of course the infinitesimal model is only providing a prediction for the distribution of trait values within families; a single realization will see substantial contributions to trait values from linkage disequilibrium (c.f. the toy example in the section on modeling Mendelian inheritance and Theorem 5). In the following figures, we compare these quantities to the detailed predictions of the infinitesimal model. The top row in Fig. 2 is a single replicate, while the bottom is the average over 300 replicates. On the left, we have the mean of the additive and dominance components and their sum; on the right, we have plotted the variance components. For a single replicate, there is indeed a substantial contribution from linkage disequilibrium. When we plot just the genic components (that is the sum over variances at each locus, ignoring the contribution from linkage disequilibrium), as expected, the picture is much smoother and we see that the predictions of the infinitesimal model are close to the values obtained by averaging over 300 replicates. Since linkage disequilibrium will dissipate rapidly, halving in each generation, it is the genic component that determines the long term evolution.

Fig. 2. Changes of the mean and variance of the additive part of the trait, the dominance part, and their sum over 50 generations of neutral evolution. The top row shows a single replicate, while the bottom row shows the average over 300 replicates using the same sequence of individuals spanning the 50 generations. The left column shows the means (G¯=A¯+D¯, A¯, D¯; black, or middle curve; blue, or bottom curve; red, or top curve), while the right column shows the variance components (VG=Var(G), VA=Var(A), VD=Var(D), VA,D=Cov(A,D); black, or top curve; blue, or middle top curve; red, or middle bottom curve; purple, or bottom curve). On the right, solid lines show the total variances and covariance, while the dashed lines show the genic component. These differ through the contribution of linkage disequilibrium, which generates substantial variation. The genic component changes smoothly, as expected with a large number (M=1,000) of loci. With M=1,000 loci, we expect the infinitesimal model to be accurate for about M∼30 generations. Simulations are made on a single pedigree with 30 individuals; variance components are measured relative to the ancestral population. The predicted values for these means and variances under the infinitesimal model are given in Equations (12)–(14) (note that the identity coefficients Fii increase through time due to genetic drift).

All components are measured relative to the base population. In practice, in natural populations, one does not have access to the ancestral population and so one measures components relative to the current population. This amounts to a change of reference Hill et al. (2006). We do not do this in our setting, as it would result in different variance components for every replicate.

In Fig. 3, we explore the relationship between the dominance deviation and inbreeding. Since we use the same pedigree for all our experiments, each individual is characterized by a single Fii (the probability of identity of the two genes at a given locus). For each of 1,000 replicates (that is independent samples of allelic types for the individuals in generation zero), we calculated the dominance deviation for each individual in the pedigree. The plot in Fig. 3 shows the dominance deviation averaged over those 1,000 replicates for each individual in the pedigree. Thus, there are 30 points in each generation, one for each individual in the population. As expected, the mean of the dominance component decreases in proportion to Fii, E[D]=−0.53Fii (recall that ι=−0.53 for our base population).

Fig. 3. The relation between the dominance deviation and the probability of identity of the two genes within an individual. There is one point for the average over 1,000 replicates for each of the 30 individuals in generations 5, 10, 20, 40 (black, or left-most group of points; blue, or second left-most group; purple, or second right-most group; red, or right-most group). (Recall that the pedigree is fixed, so identities are the same for each replicate.) The mean of D decreases as ιFii=−0.53Fii (solid line), in accordance with Equation (12).

Figure 4 shows how the (co)variance of A and D depends on identity Fii for pairs of individuals in the pedigree. As in Fig. 3, for each individual in the pedigree, A and D are calculated for each of the 1,000 replicates; Fig. 4 shows the variances and covariances of the resulting values for each of the 30 individuals in generations 5, 10, 20, and 40 and these are compared to the theoretical predictions. Note that since in the bi-allelic case σD2=ι*, the expression (14) for the variance of the dominance component reduces to

σD2(1−Fii2)+σDI2Fii.

Fig. 4. The variance and covariance of A and D versus identity Fii for individuals in the pedigree. As in Fig. 3, there are 30 points in each generation, one corresponding to each of the 30 individuals in the population. Generations 5, 10, 20, 40 (black, or left-most group of points; blue, or second left-most group; purple, or second right-most group; red, or right-most group). Here, again we use the shorter notation VA=Var(A), VD=Var(D), VA,D=Cov(A,D) and the theoretical predictions were derived in Equations (13) and (14).

Next, we consider the variances of the residuals RA and RD within families. One hundred pairs of parents were chosen at random from the population, and from each 1,000 offspring were generated. This was repeated for 10 replicates made with the same pedigree and the same set of parents; within-family variances were then averaged over replicates. In Fig. 5, in each plot there are 100 points, one for each pair of parents. The two lines correspond to least square regression (blue) and theoretical predictions (red) which can be read off from Equation (8). For readability, in the figure we use the notation VRA, VRD, and VRA,RD to denote the variance of RA, the variance of RD and the covariance between RA and RD, respectively. Using Equation (8) and the explanation below, together with the fact that σD2=ι* in our bi-allelic case, we have VRA=σA2(1−FW)/2, where FW=(Fi[1]i[1]+Fi[2]i[2])/2 is the within-individual identity averaged over parents 1 and 2;

VRD=σDI24(3F12−F1122−F112−F122)+σD24(3−F11−F22−F1122−F~1122−32F~1212);

and

VRA,RD=σADI2(F12−F(3)),

where F(3) is defined as follows:

F(3)=F112+F1222.

Fig. 5. The variance and covariance within families between the residual additive and dominance deviations RA and RD (VRA=Var(RA), VRD=Var(RD), VRA,RD=Cov(RA,RD)). One hundred pairs of parents were chosen at random from the ancestral population and from each one thousand offspring were generated. The within-family variances obtained in this way were averaged over 10 replicates (with the same pedigree and parents). Each of the 100 points in each plot corresponds to one pair of parents. The five outliers are families produced by selfing. The blue lines (or top lines) show a least-squares regression; the red lines (or bottom lines) are the theoretical predictions [see Equation (8)]. The two lines exactly coincide in the plot on the right.

The full force of our theoretical results is that even if we condition on the trait values of parents, the within-family distribution of their offspring will consist of two normally distributed components and, in particular, the variance components will be independent of the trait values of the parents. We test this by imposing strong truncation selection on the population. We retain the same pedigree relatedness, but working down the pedigree, each individual’s genotype is determined by generating two possible offspring from its parents and retaining the one with the larger trait value. In Fig. 6, we compare the results with simulations of the neutral population. Dashed lines are for the neutral simulations, solid ones for the simulation with selection. For the population under selection, we see an immediate drop in the total genetic variation, caused by the strong selection; there is significant negative linkage disequilibrium between individual loci, as predicted by Bulmer (1971). The blue is the additive component. We see that about one-third of the variance is dominance variance. The bottom row shows that the genic components are hardly affected by selection, as predicted by the infinitesimal model. With or without selection, the variance components change as a result of inbreeding.

Fig. 6. Comparison between a neutral population (dashed lines) and one subject to truncation selection (solid lines). Top row: change in means relative to the initial value (G¯=A+D¯, A¯,D¯; black, or top curves; blue, or middle curves; red, or bottom curves); middle: variances, including linkage disequilibria (top to bottom: VG=Var(A+D), VA=Var(A), VD=Var(D), VA,D=Cov(A,D); black, blue, red, purple). The bottom row is the changes to genic variances with time against predictions of the infinitesimal model. The values are averages over 300 replicates for the neutral case, 1,000 for the selected case, made with the same pedigree. There are M=1,000 loci, and thus we expect the infinitesimal model to be accurate for about M∼30 generations. Selection is made within families; for each offspring, two individuals are generated from the corresponding parents, and the one with the larger trait value retained.

Finally, Fig. 7 compares the variance components at 50 generations for neutral simulations with those with truncation selection as the number of loci increases from M=100 to M=104. Replicate simulations were generated as in Fig. 6. Under the infinitesimal model, these components should take the same values with and without selection. This is reflected in the simulations, with the covariance between the additive and dominance effects being the slowest to settle down to the infinitesimal limit.

Fig. 7. Convergence of the variance components at 50 generations, as the number of loci increases from M=100 to M=104 (same notation as in Fig. 6). Simulations with 50% truncation selection are compared with neutral simulations (solid, dashed lines). The replicate simulations were generated as in Fig. 6 (see main text). Regressions of the log absolute difference between selected and neutral variance components against ln(M) have slopes −0.62, −0.72, −0.70, −0.66 for VG, VA, VD, VA,D, respectively (see supplementary material for details). Thus, convergence is somewhat faster than M.

The infinitesimal model with dominance as a limit of Mendelian inheritance

In this section, we turn to the justification of our model as a limit of a model of Mendelian inheritance as the number M of loci tends to infinity. Although we shall focus on the distribution of the genetic components of the trait values in the pedigree, in this section we consider the general situation where the observed trait of an individual, Z~i, is the sum of a genetic component Zi and an environmental component Ei. Our mathematical assumptions on Ei are detailed in Main results below.

Our work is an extension of that of Abney et al. (2000), which in turn builds on Lange (1978). The distinctions here are that we explicitly model the component of the trait value that is shared by all individuals in a family separately from the part that segregates within that family; we identify the effect on each of these components of conditioning on knowing the trait values of the parents of the family; and we estimate the error that we are making in taking the normal approximation, thus providing information on when the infinitesimal approximation breaks down.

The fact that the genetic component of trait values within families is normally distributed is a consequence of the Central Limit Theorem. That this remains valid even when we condition on the trait values of the parents stems from the fact that knowing the trait value of an individual actually provides very little information about the allelic state at any particular locus. This in turn is because, typically, there are a large number of different genotypes that are consistent with a given phenotype. In Barton et al. (2017), this was illustrated through a simple example which can be found on p. 402 of Fisher (1918), which concerned an additive trait in a haploid population. Here we adapt that example to the model for which we performed our numerical experiments.

Suppose then that we have M bi-allelic loci. We denote the alleles at locus l by al and Al. The contributions to the trait of the three genotypes alal, alAl and AlAl are −α, −α, α respectively with probability 12−2M and they are −α, α, α with probability 12+2M. The effect size α=1/M. For simplicity, in contrast to our numerical experiments, we suppose that the probabilities of genotypes alal, alAl, AlAl are 1/4, 1/2, 1/4 respectively.

Now suppose that we observe the trait value to be k/M. What is the conditional probability that the allelic types at locus l, which we denote χl1χl2 are AlAl? For definiteness, we take M and k both to be even and l=1.

First consider the probability that the contribution to the trait value from locus 1 is +1/M. Let us write p+ for the (unconditional) probability that the contribution from locus 1 is 1/M, that is

p+=14+12(12+2M)=12(1+1M),

and p−=1−p+. Let us write Ψl/M for the contribution to the trait from locus l. We have

P[∑l=1MΨl=k|Ψ1=1]P[∑l=1MΨl=k]=P[∑l=2MΨl=k−1]P[∑l=1MΨl=k]=p+(M+k−2)/2p−(M−k)/2p+(M+k)/2p−(M−k)/2(M−1(M+k−2)/2)(M(M+k)/2)=(1+kM)12p+=(1+kM)1(1+1/M).

An application of Bayes’ rule then gives

P[χ11=A1,χ12=A1|∑l=1MΨlM=kM]=P[∑l=1MΨl=k|Ψ1=1]P[∑l=1MΨl=k]P[χ11=A1,χ12=A1]=(1+kM)1(1+1/M)P[χ11=A1,χ12=A1].

Similarly,

P[χ11=a1,χ12=a1|∑l=1MΨlM=kM]=(1−kM)1(1−1/M)P[χ11=a1,χ12=a1],

and

P[χ11=a1,χ12=A1|∑l=1MΨlM=kM]={(1+kM)(1/2+2/M)(1+1/M)+(1−kM)(1/2−2/M)(1−1/M)}×P[χ11=a1,χ12=A1].

In view of the Central Limit Theorem, we would expect a “typical” value of k to be on the order of M; conditioning has only perturbed the probability that Ψ1=1 by a factor k/M+O(1/M), which we expect to be of order 1/M. In the purely additive case, which corresponds to taking p+=p−=1/2, at the extremes of what is possible (k=±M), we recover complete information about the values of χ11, χ21; however, with dominance that is no longer true.

Notice that for the difference between the trait value of an individual and the mean over the population to be order one requires order M of the loci to be “nonrandom,” but observing the trait does not tell us which of the possible M loci these are. Similarly, performing the entirely analogous calculation for pairs of loci, and observing that

(M−2(M+k−4)/2)(M(M+k)/2)=14(1+kM)(1+k−1M−1),

we deduce that,

(19) P[χ11=A1,χ12=A1;χ21=A2,χ22=A2|∑l=1MΨlM=kM]=(1+kM)(1+k−1M−1)1(1+1/M)2×P[χ11=A1,χ12=A1;χ21=A2,χ22=A2]=P[χ11=A1,χ12=A1|∑l=1MΨlM=kM]×P[χ21=A2,χ22=A2|∑l=1MΨlM=kM]+P[χ11=A1,χ12=A1;χ21=A2,χ22=A2]×(1+kM)1(1+1/M)2(k−1M−1−kM).

For a “typical” trait value the last term in Equation (19) is order 1/M. When we sum over loci, this is enough to give a nontrivial contribution to the trait value coming from the linkage disequilibrium. However, although observing the trait of a typical individual tells us something about linkage disequilibria, it does not tell us enough to identify which of the order M2 pairs of loci are in linkage disequilibrium.

Essentially the same argument will apply to the much more general models that we develop below. In particular, for the infinitesimal model to be a good approximation, the observed parental trait values must not contain too much information about the allelic effect at any given locus, which requires that the parental traits must not be too extreme [corresponding to k in our toy model being O(M)].

In the additive case, it was enough to control the additional information that we gained about any particular locus from knowledge of the trait value in the parents. This is because, in that case, the variance of the shared contribution within a family is zero and independent Mendelian inheritance at each locus ensures that linkage disequilibria do not distort the variance of the residual component that segregates within families. With dominance, we must estimate the (nontrivial) variance of the shared component, and for this we shall see that we need to control the build up of linkage disequilibrium between pairs of loci. It will turn out that since all pairs of loci are in linkage equilibrium in the ancestral population, any given pair of loci will be approximately in linkage equilibrium for the order M generations for which the infinitesimal approximation is valid.

This does not mean that the linkage disequilibria do not affect the trait values, but because of the very many different combinations of alleles in an individual that are consistent with a given trait, observing the trait tells us very little about the allelic state at a particular locus. The allele at that locus can only ever contribute O(1/M) to the overall trait value.

As the population evolves, and we are able to observe more and more traits on the pedigree, we gain more and more information about the allele that an individual carries at a particular locus. In Barton et al. (2017), we considered an additive trait in a population of haploid individuals. In that setting, we showed that for a given individual, one does not gain any more information about the state at a given locus from looking at the trait values on the whole of the rest of the pedigree than one does from observing just the parents of that individual. In our model for diploid individuals with dominance, this is no longer the case; observing the trait values of any relatives, no matter how distant, provides some additional information about the allelic state at a locus. The difference arises from the fact that the contribution that a gene makes to the trait value of an individual depends not only on its own allelic state, but also on that of the other copy of the gene at that locus. As a result, we gain information about the allelic state in a focal individual by observing trait values in any other individuals in the pedigree with which it may be identical by descent at that locus. However, the amount of information gleaned about the allelic state of an individual from observing new individuals in the pedigree will decrease in proportion to the probability of identity, and so for distant relatives in the pedigree is very small; provided our pedigree is not too inbred, and trait values are not too extreme, we can still expect the infinitesimal model to be a good approximation for order M generations.

Environmental noise

Our derivations will depend on two approaches to proving asymptotic normality. The first, which we apply to the portion RAi+RDi of the trait values, uses a generalized Central Limit theorem (which allows for the summands to have different distributions), which provides control over the rate of convergence as M→∞. (It is this control that tells us for how many generations we can expect the infinitesimal model to be valid.) However, the Central Limit Theorem guarantees only the rate of convergence of the cumulative distribution function of the normalized sum of effects at different loci. Our proofs exploit convergence to the corresponding probability density function, which may not even be defined. To get around this, we can follow the approach of Barton et al. (2017) and make the (realistic) assumption that rather than observing the genetic component of a trait directly, the observed trait has an environmental component with a smooth density. This results in the trait distribution having a smooth density which is enough to guarantee the faster rate of convergence. In addition to the benefit in terms of regularity of the trait distribution, an environmental noise with a smooth distribution also reinforces the property that observing the trait value gives us very little information on the allelic state at a given locus: a continuum of combinations of genetic and environmental components may have led to the observed trait, in which each given locus contributes an infinitesimal amount. (To ensure sufficient regularity of the trait density, we could instead make the assumption that the distribution of allelic effects at every locus has a smooth probability density function.) The approach to proving asymptotic normality of the shared component uses an extension of Stein’s method of exchangeable pairs. Once again in the presence of environmental noise (to ensure that the trait distribution has a smooth density) we recover convergence with an error of order 1/M.

If the environmental component is taken to be normally distributed, then exactly as in Barton et al. (2017), we can adapt our application of Theorem C1 in Appendix C to write down the conditional distribution of the genetic components given observed traits; i.e. traits distorted by a small environmental noise, c.f. Remark F2.

Assumptions and notation

Recall that we assume that in generation zero, the individuals that found the pedigree are unrelated and sampled from an ancestral population in which all loci are assumed to be in linkage equilibrium. The allelic states at locus l on the two chromosomes drawn from the ancestral population will be denoted χ^l1,χ^l2. They are independent draws from a distribution on possible allelic states that we denote by ν^l(dx). Without loss of generality, by replacing ϕl(χ^l1,χ^l2) by

ϕl(χ^l1,χ^l2)−E[ϕl(χ^l1,χ^l2)|χ^l1]−E[ϕl(χ^l1,χ^l2)|χ^l2]+E[ϕl(χ^l1,χ^l2)],

and observing that the second and third terms on the right-hand side are functions of χ^l1 and χ^l2, respectively, which we may therefore subsume into ηl(χ^l), we may assume that for any value x′ of the allelic state at locus l, the conditional expectation

(20) E[ϕl(χ^l,x′)]=∫ϕl(x,x′)ν^l(dx)=0=E[ϕl(x′,χ^l)].

As a consequence, partitioning over the possible values of χ^l2, we have that the cross variation term

(21) E[ηl(χ^l1)ϕl(χ^l1,χ^l2)]=∫E[ηl(x′)ϕl(x′,χ^l2)]ν^l(dx′)=∫ηl(x′)E[ϕl(x′,χ^l2)]ν^l(dx′)=0.

With this modification of ϕl(x,x′),

(22) E[ϕl(χ^l1,χ^l2)]=0.

Moreover, still without loss of generality, by absorbing the mean into z¯0, we may assume that

(23) E[ηl(χ^l)]=∫ηl(x)ν^l(dx)=0.

In this notation, the genetic component of the trait of an individual in the ancestral population (which we denote by Z^ to make it clear that the following property is specific to individuals in generation 0) is

Z^=z¯0+1M∑l=1M(ηl(χ^l1)+ηl(χ^l2)+ϕl(χ^l1,χ^l2)),

and by Equations (22) and (23), we have E[Z^]=z¯0.

We assume that the scaled allelic effects ηl, ϕl are bounded; |ηl|, |ϕl|≤B, for all l. We also assume that all the quantities in the top part of Table 1 exist in the limit as M→∞.

Inheritance

We now need some notation for Mendelian inheritance. Recall that i[1] and i[2] are the labels of the parents of individual i in our pedigree, each of which contributes exactly one gene at each locus in a given offspring. Mendelian inheritance translates into the property that the gene passed on by parent i[1] was the one inherited from its own “first” parent (i[1])[1] with probability 1/2, or from its “second” parent (i[1])[2] with probability 1/2. Even though we do not distinguish between males and females, it is convenient to think of the chromosomes in individual i as being labeled 1 and 2, according to whether they are inherited from i[1] or i[2]. In particular, χli[1],1 and χli[1],2 will denote the allelic states of the two genes at locus l in parent i[1], respectively inherited from its own “first” and “second” parent. Again following the conventions of Barton et al. (2017), extended to account for the fact that we are now considering diploid individuals, we use independent Bernoulli(1/2) random variables, Xli, Yli to determine the inheritance of genes 1 and 2, respectively, at locus l in individual i. Thus, Xli=1 if the allelic state of gene 1 at locus l in individual i is inherited from gene 1 in i[1], and Xli=0 if it is inherited from gene 2 in i[1]. Likewise, Yli=1 if the allelic state of gene 2 at locus l in individual i is inherited from gene 1 in i[2], and Yli=0 if it is inherited from gene 2 in i[2].

In this notation, the trait of individual i in generation t is given by

(24) Zi=z¯0+Ai+Di+1M∑l=1M{(Xli−12)ηl(χli[1],1)+(12−Xli)ηl(χli[1],2)

(25) (25)+(Yi−12)ηl(χli[2],1)+(12−Yi)ηl(χli[2],2)}+1M∑l=1M{(XliYli−14)ϕl(χli[1],1,χli[2],1)(26)+(Xli(1−Yli)−14)ϕl(χli[1],1,χli[2],2)+((1−Xli)Yli−14)ϕl(χli[1],2,χli[2],1)

(27) +((1−Xli)(1−Yli)−14)ϕl(χli[1],2,χli[2],2)},

where

(28) Ai=12M∑l=1M(ηl(χli[1],1)+ηl(χli[1],2)+ηl(χli[2],1)+ηl(χli[2],2))

and

(29) Di=14M∑l=1M{ϕl(χli[1],1,χli[2],1)+ϕl(χli[1],1,χli[2],2)+ϕl(χli[1],2,χli[2],1)+ϕl(χli[1],2,χli[2],2)}.

The terms Ai and Di are shared by all descendants of the parents i[1] and i[2]. In the third section of this paper, we presented the mean and variance of their sum, conditional on the pedigree P(t). The sums (24)+(25) and (26)+(27) comprise what we previously called RAi and RDi, respectively; each has mean zero. They capture the randomness of Mendelian inheritance. They are uncorrelated with Ai+Di. Again, in a previous section we gave expressions for the variances and covariance of RAi and RDi in terms of the ancestral population and identities generated by the pedigree. These calculations allowed us to identify the mean and variance of the parts Ai+Di and RAi+RDi in terms of the classical quantities of quantitative genetics in Table 1. Since we are assuming unlinked loci, the asymptotic normality of these quantities when we condition on the pedigree, but not on the trait values within that pedigree, is an elementary application of Theorem D2 in Appendix D, a generalized Central Limit Theorem which allows for nonidentically distributed summands.

In Barton et al. (2017), we showed that in the purely additive case, the vector (RAi)i=1Nt, which determines the joint distribution of the trait values within families in generation t (recalling that in the additive case RDi=0), is asymptotically a multivariate normal, even when we condition not just on the pedigree relatedness of the individuals in generation t, but also on knowing the observed trait values of all individuals in the pedigree up to generation t−1, which we denote by Z~(t−1) (notice the difference between this notation and the notation Z~t for the observed trait of an individual living in generation t). Our main result extends this to include dominance, at least under the assumption that the ancestral population was in linkage equilibrium.

With dominance, the expression for the distribution of the mean and variance–covariance matrix of the multivariate normal Z1,…,ZNt conditioned on the pedigree up to generation t and some collection of the observed trait values of individuals in that pedigree up to generation t−1 is a sum of the quantities of classical quantitative genetics in Table 1, weighted by four-way identities and deviations of trait values from the mean. In principle, they can be read off from Theorem C1 in Appendix C.

We will focus on proving that conditional on knowing just the trait values of the parents of individual i and the pedigree, the components (Ai+Di) and (RAi+RDi) are both asymptotically normal, but we explain why our proof allows us to extend to the case in which we also know trait values of other individuals. The importance (and surprise) is that given the pedigree relationships between the parents and classical coefficients of quantitative genetics for a base population (assumed to be in linkage equilibrium), knowing the traits of the parents distorts the distribution of their offspring in an entirely predictable way. In particular, this is what we mean when we say that the infinitesimal model continues to hold even with dominance.

The extra challenge compared to the additive case is that, in contrast to the part RAi+RDi, where Mendelian inheritance ensures independence of the summands corresponding to different loci even after conditioning on trait values, when we condition on trait values the terms in Ai+Di will be (weakly) dependent and proving a Central Limit Theorem becomes more involved.

Main results

Recall that the trait values that we observe, and therefore on which we condition, are the sum of a genetic component and an independent environmental component; that is, the observed trait value is

Z~i:=Zi+Ei,

where, for convenience, the {Ei} are independent N(0,σE2)-valued random variables. We suppose that the environmental noise is shared by individuals in a family (so we can think of it as part of the component Ai+Di of the trait value, whose distribution therefore also has a smooth density).

We write Nt for the number of individuals in the population in generation t, (Zt1,…,ZtNt) for the corresponding vector of trait values, and P(t) for the pedigree up to and including generation t. A simple application of the Central Limit Theorem gives that

(Zt1,…,ZtNt)|P(t)

is asymptotically distributed as a multivariate normal random variable as M→∞. More precisely, let (β1,β2,…,βNt)∈RNt, and write Zβ=∑i=1NtβiZti, then using Theorem D2,

|P[Zβ−E[Zβ]Var(Zβ)≤z]−N(z)|≤CMVar(Zβ)(1+C~Var(Zβ)),

for suitable constants C,C~ (which can be made explicit), where N(z) is the cumulative distribution function for a standard normal random variable. The mean and variance of Zβ can be read off from Equations (9), (10), and (11).

Our main results concern the components of the trait values of offspring when we condition on the observed trait values of their parents. The following result follows in essentially the same way as the additive case of Barton et al. (2017).

Theorem 4 The conditioned residuals (RAi+RDi)|P(t),Z~i[1], Z~i[2] are asymptotically normally distributed, with an error of order 1/M. More precisely, for all z∈R,|P[RAi+RDiVar(RAi+RDi)≤z|P(t),Z~i[1],Z~i[2]]−N(z)|≤1MC′Var(RAi+RDi)(1+C′~Var(RAi+RDi))×(1+C(i[1],i[2]))

whereC(i[1],i[2])=C″|Z~i[1]−E[Z~i[1]|P(t−1)]|Var(Z~i[1])+C″|Z~i[2]−E[Z~i[2]|P(t−1)]|Var(Z~i[2])+C‴1Var(Z~i[1])p(Var(Z~i[1]),|Zi[1]−E[Zi[1]|P(t−1)]|)×(1+1Var(Z~i[1]))+C‴1Var(Z~i[2])p(Var(Z~i[2]),|Zi[2]−E[Zi[2]|P(t−1)]|)×(1+1Var(Z~i[2])),

and we have used p(σ2,x) to denote the density at x of a mean zero normal random variable with variance σ2. The constants C′, C′~, C″, C‴ depend only on the bound B on the scaled allelic effects. The variances in the expressions above are all calculated conditional on P(t−1), but not on observed parental trait values.

Put simply, the normal approximation is good to an error of order 1/M; the constant in the error term will be large, meaning that the approximation will be poor, if the within-family variance somewhere in the pedigree is small or if the observed trait values are very different from their expected values. Just as in the additive case, we could prove an entirely analogous result when we condition on any number of observed trait values in the pedigree, except that with dominance this is at the expense of picking up an extra term in the error for each observed trait value on which we condition. The justification required for this is provided by Appendix H.

What is at first sight more surprising is that the shared component of the trait value within a family, i.e. the random variable A+D+E, is also asymptotically normally distributed, even when we condition on observed parental trait values. Note that the randomness of the shared component comes from the fact that the allelic states underlying the parental traits are still random (they are unobserved). In the case of a purely additive trait, it turns out that the shared component can be simply expressed as the average of the two parental traits and therefore conditioning on these traits renders the shared contribution totally deterministic, but such a simplification no longer occurs when we add dominance, due to the nonlinearity of the allelic contributions in D [see Equation (29)]. Our proof of normality uses the fact that we consider the environmental noise to be shared by individuals within the family; in this way we can guarantee that the shared component of the observed trait value also has a smooth density.

We are only going to prove the result for the shared component of a family in generation one that was produced by selfing (i[1]=i[2]). In what follows, for a given function h we write ‖h‖ for the supremum norm of h, and Nμ,σ2(h) for the integral of h with respect to the distribution of an N(μ,σ2) random variable (whenever this quantity makes sense):

Nμ,σ2(h)=12πσ2∫−∞+∞h(z)e−(z−μ)2/(2σ2)dz.

Theorem 5 Let W=A+D+E denote the shared component of the trait value in a family in generation one. Let h be an absolutely continuous function with ‖h′‖<∞, then(30) |E[h(W)|i[1]=i[2],Z~i[1]]−NμW,σW2(h)|≤C‖h′‖M,

where μW is given by Equation (F5), and σW2 is the sum of the variance of the environmental noise and the expression in Equation (F21).

Remark 6 Although we only prove that Ai+Di+Ei is asymptotically normal in this special case of an individual in generation one that is produced by selfing, the same arguments will apply in general. However, the expressions involved become extremely cumbersome. By considering selfing, we capture all the complications that arise in later generations (when distinct parents may nonetheless be related).

We do not record the exact bound on the constant C. It takes the same form as the error function C in Theorem 4, except that the constants C′, C′~, C″, C‴ depend on the inbreeding depression ι, as well as the bound B on the scaled allelic effects. In particular, just as there, the asymptotic normality will break down if the trait value of the parent is too extreme, or if the variance of the trait values among offspring is too small.

Since we are assuming that the environmental noise has a smooth density, convergence in the sense of Equation (30) is sufficient to deduce that the cumulative distribution of Ai+Di+Ei converges.

In Fig. 8, we show the cumulative distribution functions of the additive and dominance parts of the shared and residual components of trait within 10 families after 20 generations of neutral evolution, with M=1,000 loci. All 10 within-family distributions of RA, RD are close to Gaussian; they vary somewhat in slope, since families vary in identity coefficients (see Fig. 5), but this is not apparent in these plots. The normal approximation is better for the residual components than for the shared component. This may be due to the fact that the random variables encoding Mendelian inheritance at different loci are independent and identically distributed, which makes the summands in the expressions for RA and RD more weakly dependent than the summands in A and D, leading to faster convergence to a Gaussian distribution. This also explains why we need a more elaborate approach to show convergence of the shared parts to Gaussians.

Fig. 8. The distributions of the residual (top row: RA , RD) and shared (bottom row: A, D) components of phenotype (M=1,000 loci); for each, the cumulative distribution function is plotted as standard deviations of a Gaussian, z, so that a normal distribution appears as a straight line. These are calculated from families of 1,000 offspring, from multiple pairs of parents, each replicated 10 times, drawn after 20 generations without selection. The residuals are calculated by subtracting values from the family mean, and pooling across the 10 replicates. Thus, for each family there are 10,000 values; the cumulative distribution function is shown for 10 pairs of parents, in 10 colors. The shared component is calculated by taking the mean of each family, and pooling across 100 pairs of parents and across the 10 replicates. Thus, for each plot there are 1,000 points. There is now some deviation from a Gaussian.

Strategy of the derivation

Our first task will be to show that conditional on the pedigree, the distribution of the trait values in generation t is approximately multivariate normal (with an appropriate error bound). Since Mendelian inheritance ensures that (before we condition on knowing any of the previous trait values in the pedigree) the allelic states at different loci are independent, this is a straightforward application of a generalized Central Limit Theorem (generalized because the summands are not required to all have the same distribution). Just as in Barton et al. (2017), we can keep track of the error that we are making in assuming a normal approximation at each generation. In this way we see that, under our assumptions, the infinitesimal model can be expected to be a good approximation for order M generations.

The same Central Limit Theorem guarantees that the joint distribution of (Zi[1],Zi[2],Ai+Di) is asymptotically normally distributed as the number of loci tends to infinity. This certainly suggests that the conditional distribution of Ai+Di given Zi[1], Zi[2] should be (approximately) normal with mean and variance predicted by standard results on conditioning a multivariate normal distribution on some of its marginals (which we recall in Theorem C1). However, this is not immediate. It is possible that the conditioning forces the distribution on to the part of our probability space where the normal approximation breaks down.

To verify that the conditional distribution is asymptotically normal, we shall show that observing the trait value of an individual provides very little information about their allelic state at any particular locus, or any particular pair of loci, and consequently conditioning on parental trait values provides very little information about allelic states in their offspring. This is (essentially) achieved through an application of Bayes’ rule, although some care is needed to control the cumulative error across loci. We use this to calculate the first and second moments of Ai+Di conditional on Z~i[1], Z~i[2]. The fact that they agree with the predictions of Theorem C1 depends crucially on the assumption that dominance is “balanced,” in the sense that the inbreeding depression ι is well defined. This quantity enters not just in the expression for the expected trait value of inbred individuals, but also in our error bounds, c.f. Remark F4.

Of course checking that the first two moments of the conditional distribution of Ai+Di are (approximately) consistent with asymptotic normality is not enough to prove that the conditioned random variable is indeed (approximately) normal. Moreover, we cannot apply our generalized Central Limit Theorem to this term. Instead we use a generalization of Stein’s method of “exchangeable pairs” (outlined in Appendix D), which relies on our ability to control the (weak) dependence between the contributions to Ai+Di from different loci that is induced by the conditioning. We present the details in the case of identical parents (which is the case in which normality is most surprising) in Appendix G.

We only present our results in the case in which we condition on the parental traits of a single individual in generation t. Just as in the additive case, this can be extended to conditioning on any combination of traits in the pedigree up to generation t−1, but the expressions involved become unpleasantly complex. Instead of writing them out, we content ourselves with explaining the only step that requires a new argument. We must show that knowing the traits of all individuals up to generation t−1 does not provide enough information about the allelic states at any particular locus in an individual in generation t to destroy the asymptotic normality of its trait value. This is justified in Appendix H using the fact that, because of Mendelian inheritance, the amount of information gleaned about an allele carried by individual i from looking at the trait value of one its relatives, is proportional to the probability of identity with that individual as dictated by the pedigree.

Asymptotic normality conditional on the pedigree

We first illustrate the application of the generalized Central Limit Theorem by showing that in the ancestral population, the distribution of (Z01,…,Z0N0) is multivariate normal with mean vector (z¯0,…,z¯0) and variance–covariance matrix (σA2+σD2)Id, where Id is the identity matrix and σA2 and σD2 were defined in Table 1.

To prove this, it is enough to show that for any choice of β=(β1,…,βN0)∈RN0,

∑j=1N0βjZj→Zβ,

where Zβ is normally distributed with mean z¯0∑j=1N0βj and variance (σA2+σD2)∑j=1N0βj2. We apply Theorem D2, due to Rinott (1994), which provides control of the rate of convergence as M→∞. It is convenient to write ‖β‖1=∑j=1N0|βj| and ‖β‖22=∑j=1N0βj2. Let us write

Ψl=(ηl(χ^l1)+ηl(χ^l2)+ϕl(χ^l1,χ^l2)),

and we abuse notation by writing Ψlj for this quantity in the jth individual in generation zero. Set El=∑j=1N0βjΨlj. Recalling our assumption that all ηl and ϕl are bounded by some constant B, so that the sum of the scaled effects at each locus is bounded by 3B, we have that |El| is bounded by 3B‖β‖1 for all l. Moreover, since the individuals that found the pedigree are assumed to be unrelated and sampled from an ancestral population in which all loci are in linkage equilibrium, using Equations (22) and (23), we find that

E[∑l=1MEl]=0,Var(∑l=1MEl)=M‖β‖22(σA2+σD2).

Theorem D2 then yields

|P[∑i=1N0βi(Zi−z¯0)‖β‖2σA2+σD2≤z]−N(z)|≤1M‖β‖2σA2+σD2{12π3B‖β‖1+16‖β‖2σA2+σD2(3B)2‖β‖12+10(1‖β‖22(σA2+σD2))(3B‖β‖1)3}.

Here, N is the cumulative distribution function of a standard normal random variable. The right-hand side can be bounded above by

(31) C(‖β‖1)‖β‖2MσA2+σD2(1+1‖β‖22(σA2+σD2)),

for a suitable constant C. In particular, taking βk=0 for k≠j and βj=1, we read off that the rate of convergence to the normal distribution of Z0j as the number of loci tends to infinity is order 1/M. Note that the normal approximation is poor if the variance σA2+σD2 is small.

Exactly the same argument shows that the distribution of (Z1,…,ZNt) of the individuals in generation t converges to that of a multivariate normal, with mean vector (z¯0+ιF11,…,z¯0+ιFNtNt) and variance–covariance matrix determined by Equations (10) and (11).

Our proof of asymptotic normality of Ai+Di conditional on the observed trait values of parents will exploit that the joint distribution of (Ai+Di,Zi[1],Zi[2]) is asymptotically normal, also with an error of order 1/M. This time we show that β1Zi[1]+β2Zi[2]+β3(Ai+Di) is asymptotically normal for every choice of the vector (β1,β2,β3)∈R3. We apply Theorem D2 with

E~l=β1Ψl(i[1])+β2Ψl(i[2])+β3Φli,

where

Ψl(i[1])=ηl(χli[1],1)+ηl(χli[1],2)+ϕl(χli[1],1,χli[1],2),

with a symmetric expression for Ψl(i[2]), and

Φli=12(ηl(χli[1],1)+ηl(χli[1],2)+ηl(χli[2],1)+ηl(χli[2],2))+14(ϕl(χli[1],1,χli[2],1)+ϕl(χli[1],1,χli[2],2)+ϕl(χli[1],2,χli[2],1)+ϕl(χli[1],2,χli[2],2)).

Theorem D2 then shows that the difference between the cumulative distribution function of β1Zi[1]+β2Zi[2]+β3(Ai+Di) and that of a normal random variable with the corresponding mean and variance can be bounded by Equation (31) with ‖β‖22(σA2+σD2) replaced by Var(β1Zi[1]+β2Zi[2]+β3(Ai+Di)), which can be deduced from the expressions for the variance and covariance of Ψli[1], Ψli[2] and Φli that are calculated in Appendix B and recorded in Equations (10), (11), and (16).

Conditioning on trait values of the parents

We suppose that for each i, we know the parents of the individual i and their trait values Zi[1] and Zi[2]. We shall treat the shared components (Ai+Di) and the residuals (RAi+RDi) separately. Both will converge to multivariate normal distributions which are independent of one another.

Mendelian inheritance ensures that the contributions to RAi+RDi from different loci are independent and so normality becomes an easy consequence of Theorem D2 once we have shown that the information gleaned from knowing the trait values only perturbs the distribution by order 1/M. This is checked in Equation (F7) and the proof then closely resembles the proof in the additive setting of Barton et al. (2017) and so we omit the details.

The proof that (Ai+Di) is normal is more involved as once we condition on the trait values in the parents, the contributions Φli for l=1,…,M will all be (weakly) correlated. Our approach uses an extension of Stein’s method of exchangeable pairs which we recall in Appendix D and apply to our setting in Appendix G. This calculation is more delicate, but the key is that our conditioning induces very weak dependence between loci. The deviation from normality is controlled by

1P[Z~i[1]=z1,Z~i[2]=z2,Ai+Di+Ei=w]×∂∂z1P[Z~i[1]=z1,Z~i[2]=z2,Ai+Di+Ei=w],

and the corresponding quantity for the partial derivative with respect to z2 (both to be interpreted as ratios of densities) evaluated at Z~i[1], Z~i[2] respectively. (We recall that Z~ denotes observed trait value.) The normal approximation will break down if the trait values are too extreme or if the pedigree is too inbred.

Discussion

The essence of the infinitesimal model is that the distribution of a polygenic trait across a pedigree is multivariate normal. Necessarily, if some individuals are selected (that is, if we condition on their trait values), there can be an arbitrary distortion away from Gaussian across the population. However, conditional on parental values and on the pedigree, offspring within each family still follow a Gaussian distribution. This was shown in Barton et al. (2017) in the purely additive case, and is extended here to the case with dominance; the only difference being that with dominance, the part of the trait shared by all siblings, A+D, is now still random even when conditioning on the parental traits (observing the parental traits does not give us full information on the contribution of the parental alleles to the average offspring trait as it did in the purely additive case), and the most difficult part of our analysis consists in showing that this shared contribution is also Gaussian. Our results strongly rely on our assumption that inbreeding depression, ι, is finite (it is zero in the purely additive case). Armed with these results, the classic theory for neutral evolution of quantitative traits can be used to predict evolution, even under selection. Theorems 4 and 5 show that this infinitesimal limit holds with dominance, at least over timescales of order square root of the number of loci. Indeed, they show that conditional on the parental traits, the distance between the distributions of the components of the offspring trait and a normal distribution is of the order of 1/M. Hence, the distance between the trait distribution of an individual and the infinitesimal approximation increases in every generation by a factor of order 1/M, and the error bound becomes macroscopic (i.e. order 1) after of the order of M generations.

Our work provides some mathematical justification for the ubiquity of the Gaussian, and the empirical success of quantitative genetics—a success which is remarkable, given the complex interactions that underlie most traits. The limit is not universal: a nonlinear transformation of a Gaussian trait leads to a non-Gaussian distribution, and failure of the infinitesimal model. This is because epistatic and dominance interactions then have a systematic direction, which violates the terms of the Central Limit Theorem. (Recall that in our toy example in the section on modeling Mendelian inheritance, we needed a “balance” in the dominance component, which we see reflected in our main results in the requirement that ι be well defined.) Nevertheless, if the population is restricted to a range that is narrow relative to the extremes that are genetically possible, then the infinitesimal model may be accurate, even if the genotype-phenotype map is not linear. This links to another way to understand our results: if very many genotypes can generate the same phenotype, then knowing the trait value gives us negligible information about individual allele frequencies. To put this another way, the infinitesimal limit implies that selection on individual alleles is weak relative to random drift (Nes∼1), so that neutral evolution at the genetic level is barely perturbed by selection on the trait (Robertson 1960).

If traits truly evolve in this infinitesimal regime, then it will be impossible to find any genomic trace of their response to selection. This extreme view is contradicted by finding an excess of “signatures” of selection in candidate genes, though it might nevertheless be that these signals are generated by alleles with modest Nes, such that the infinitesimal model remains accurate for the trait. Indeed, Boyle et al. (2017) argue that the very large numbers of single nucleotide polymorphisms that are typically implicated in genome-wide association studies for complex traits implies an “omnigenic” view, in which trait variance is largely due to genes with no obvious functional relation to the trait. Frequencies of nonsynonymous and synonymous mutations suggest that selection on deleterious alleles is typically much stronger than drift (Nes≫1; Charlesworth 2015). However, it might still be that selection on the focal trait is comparable with drift, even if the total selection on alleles is much stronger. Whether the infinitesimal model accurately describes trait evolution under such a pleiotropic model is an interesting open question.

In principle, we can simulate the infinitesimal model exactly, by generating offspring from the appropriate Gaussian distributions. For the additive case, this is straightforward, since we only need follow the breeding value of each individual, and the matrix of relationships amongst individuals (e.g. Barton and Etheridge 2011, 2018). However, to simulate the infinitesimal model with dominance, we need to track four-way identities, which is only feasible for small populations (<30, say).

We have not set out the extension of the infinitesimal model to structured populations in detail. In principle, this just requires that we track the identities within and between the various classes of individual. One motivation for the present theoretical work was to extend our infinitesimal model of “evolutionary rescue” (Barton and Etheridge 2018) to include inbreeding depression and partial selfing. This should be feasible, provided that we do not need to track identities between specific individuals, but instead, group individuals according to the time since their most recent outcrossed ancestor—an approach applied successfully by Sachdeva (2019). Already, Lande and Porcher (2015) applied the infinitesimal model to a deterministic model of partial selfing, while Roze (2016) analyzed an explicit multilocus model of partial selfing, allowing for dominance and drift, assuming that all loci are equivalent, and that linkage disequilibria are weak.

One of the most obviously unreasonable assumptions of the classical infinitesimal model, and the extension described here, is that there are an infinite number of unlinked loci. Santiago (1998) showed how loose linkage could be approximated by averaging over pairwise linkage disequilibria. In the additive case, the infinitesimal model can be defined precisely for a linear genome, by assuming that very many genes are spread uniformly over the genome (Sachdeva and Barton 2018). The techniques used in our approach are not robust to (even moderately) high levels of linkage, as groups of genes passed on together will decrease the number of “independent” units of heritable contributions to the trait value, leading to an effective number of loci Meff too low for the Gaussian approximation to be valid (or more precisely, for the bound between the trait distribution and the appropriate Gaussian distribution in Theorems 4 and 5 to be small). In this case, one needs to consider explicit models of recombination that are out of the scope of this work.

The main value of the infinitesimal model may be to show that trait evolution depends on only a few macroscopic parameters; even if we still make explicit multilocus simulations, this focuses attention on those key parameters, and gives confidence in the generality of our results. Quantitative genetics has developed quite separately from population genetics. Although the theoretical synthesis half a century ago (Robertson 1960; Bulmer 1971; Lande 1975) stimulated much subsequent work (empirical as well as theoretical), the failure to find a practicable approximation for the evolution of the genetic variance (e.g. Turelli and Barton 1994) was an obstacle to further progress. The infinitesimal model provides a justification for neglecting the intractable effects of selection on the variance components, and treating them as evolving solely due to drift and migration. This approach may be helpful for understanding evolution in the short and even medium term.

Acknowledgments

We thank the two Reviewers and the Associate Editor for their very useful detailed comments, which helped us to improve the presentation of the results.

Data availability

The code and data produced for this work and used in this article can be found in the public repository (Barton 2023).

Funding

NHB was supported in part by ERC Grants 250152 and 101055327. AV was partly supported by the chaire Modélisation Mathématique et Biodiversité of Veolia Environment—Ecole Polytechnique—Museum National d’Histoire Naturelle—Fondation X.

Appendices

The appendices are organized as follows. Appendix A discusses a simple algorithm to compute identity coefficients. In Appendix B, we derive the mean and covariances of the shared and residual parts of the offspring trait knowing the pedigree (but not the parental traits). In Appendix C, we recall a standard result for conditioning multivariate normal random vectors on their marginal values, while in Appendix D we recall the generalized Central Limit Theorems that will be needed to obtain the normal distribution of the offspring trait components conditional on the parental traits. In Appendix E, we prove some key lemmas on conditional allelic distributions that we use in Appendix F to compute the mean and variance of trait values conditional on the pedigree and on parental traits. The convergence of the shared component of the trait to a Gaussian random variable, as the number of loci tends to infinity, is obtained in Appendix G. Finally, in Appendix H we investigate how information accumulates when we condition on knowing more ancestral traits than those of the parents.

Appendix A: Calculating identity coefficients

Recursions for pairwise identity by descent

Two-way identities are readily expressed as solutions to a recurrence. The recursion for F can be written in terms of a pedigree matrix, Pi,k(t), which gives the probability that a gene in individual i in generation t came from parent k in generation (t−1); each row has two nonzero entries each with value 1/2 (the entries corresponding to the indices of the two parents, since the gene may have been inherited from either parent with the same probability), unless the individual is produced by selfing, in which case there is a single entry with value 1 (that corresponding to the index of the single parent). Observe that the matrices P(t) are totally determined by knowledge of the pedigree. In contrast to Barton et al. (2017), where we focused on haploids, here we necessarily have to deal with diploids. For diploids, the recursion for F is

(A1) Fij(t)=∑k,lPi,k(t)Pj,l(t)Fkl*(t−1),

where

Fkl*=Fklif k≠l,Fkk*=12(1+Fkk).

The quantity Fkl* is the probability of identity of two genes drawn independently from individuals k and l (this independent drawing corresponds to Mendelian inheritance); if k=l, then we may either pick the same gene twice, which happens with probability 1/2 (and since the two genes are identical, they are also identical by descent), or pick the two genes of individual k, again with probability 1/2, and their probability of identity by descent is then Fkk by definition. Restating (A1) in words, the probability that a gene taken in individual i and a gene taken in individual j, both in generation t, are identical by descent is equal to the sum over all potential pairs (k,l) of parents in the previous generation (t−1) of the probability that the gene in i descends from k, the gene in j descends from l and that the “parental” genes in k and l are themselves identical by descent.

Calculating two-, three-, and four-way identities

Several papers have developed algorithms for calculating identity coefficients, given a pedigree (Karigl 1981; Abney 2009; García-Cortés 2015; Kirkpatrick et al. 2019). These assume a single genetic locus, and primarily consider the nine condensed identity coefficients of Fig. B1 that describe the relationship between two diploid individuals. This body of work has developed algorithms that can efficiently calculate identity coefficients involving two individuals, across large pedigrees. Karigl (1982) considers (but does not implement) calculation of identities amongst more than two individuals.

Here, we define and implement a (fairly) simple algorithm that deals with multiple sets of genes across multiple individuals. The corresponding code in Mathematica can be found in supplementary material (Barton 2023). This is unlikely to be as efficient as existing algorithms for identities amongst one set of genes across two individuals; it is limited by the need to calculate and store identities amongst very many sets of ancestral genes, corresponding to the very many routes by which genes may descend through the pedigree.

First, we establish our notation. The two genes in each individual each receive a separate label. Thus, a gene in individual i will have label i={i,1} or i={i,2}. Sets of genes will be generically denoted by S={i1,…,ik}. We define F[S1,S2,…,Sn] to be the probability that the genes contained in each set S1, S2, … , Sn are identical by descent, tracing back to n distinct founders in the ancestral population. For example, F[{i1},{i2,i3},{i4,i5}] is the probability that these three sets of genes, S1={i1}, S2={i2,i3} and S3={i4,i5}, each trace back to three distinct founders: one ancestral to i1, another one ancestral to i2 and i3, and a last one ancestral to i4 and i5. Necessarily, F[{i}]=1 (a single gene traces back to a unique founder), and the probability of identity of genes i1 and i2 satisfies F[{i1,i2}]=1−F[{i1},{i2}]. Identities in generation t are denoted Ft.

Given the pedigree, the identities are defined recursively; Ft is a linear combination of identities Ft−1 in the previous generation. Here, we simply outline the algorithm. A detailed explanation in terms of the Mathematica code is in the supplementary material (Barton 2023).

In generation t=0 all individuals are assumed unrelated and so F0[S1,…,Sn] is set to be 1 if each Sk comprises a single gene and these n genes are all distinct. Otherwise it is set to zero.

The algorithm proceeds in two steps, first identifying the possible parents from which each gene is descended and then the possible genes within that parent. In this way, a list of all possible scenarios is generated, with each scenario having equal probability. A slight twist here is that if a set contains a single gene in a given individual, that gene traces back to one or other parent of the individual, with equal probability; two genes in the same individual must trace back to the two parents, although those may be the same individual if there is selfing. This list contains many permutations that are equivalent, differing only by order; these are tallied to reduce the number of configurations that need to be stored, resulting in a weighted list. This gives a recursion back to the founder generation. The number of generations and size of pedigree is limited by the amount of memory needed to store the intermediate lists.

Appendix B: Conditioning on the pedigree

In this section, we illustrate how to recover the expressions for the mean and variance of the two parts (Ai+Di) and (RAi+RDi) of the trait of individual i from identity coefficients of its parents i[1] and i[2] and the classical coefficients of Table 1. Covariances between families are calculated in the same way. We also calculate the covariance between (Ai+Di) and Zi[1] and Zi[2] (given the pedigree) which will be important for establishing the effect of conditioning on the trait values of the parents. Although these expressions are well known, it seems to be hard to find an explicit derivation such as that presented here. Note that at this stage we are only conditioning on the pedigree, not on the observed trait values and the results in this section do not require us to assume the presence of an environmental noise term.

Fig. B1. All possible four-way identities. The dots represent the four genes across the two parents (each parent corresponding to a row) and lines indicate identity (c.f. Abney et al. 2000).

Notation

Throughout this section, we are going to be calculating quantities conditional on the pedigree. We shall suppress that in our notation.

Mean and variance of Ai+Di

The contribution to the trait Zi from the lth locus is determined by the four alleles χli[1],1, χli[1],2,χli[2],1, and χli[2],2 and the independent Bernoulli random variables Xli and Yli. The mean and variance of (Ai+Di) and (RAi+RDi) will depend on which combinations of these alleles are identical. First, we introduce some notation for the nine possible identity classes. In Fig. B1, the two copies of each gene in each individual are represented by two (horizontally adjacent) dots. Lines between dots represent identity by descent. It is convenient to think of the genes within an individual as being ordered.

Let us define

(B1) Φ(l)=12(ηl(χli[1],1)+ηl(χli[1],2)+ηl(χli[2],1)+ηl(χli[2],2))+14(ϕl(χli[1],1,χli[2],1)+ϕl(χli[1],1,χli[2],2)+ϕl(χli[1],2,χli[2],1)+ϕl(χli[1],2,χli[2],2)),Ψl(i[1])=ηl(χli[1],1)+ηl(χli[1],2)+ϕl(χli[1],1,χli[1],2),

and

Ψl(i[2])=ηl(χli[2],1)+ηl(χli[2],2)+ϕl(χli[2],1,χli[2],2).

For each of the nine possible identity classes between i[1] and i[2], we calculate two quantities from which the mean and variance of (Ai+Di) will readily follow.

id. stateE[1M∑l=1MΦ(l)2|Δ⋅]E[1M∑l=1MΦ(l)|Δ⋅]−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−Δ12σA2+2σADI+σDI2+ι*ιΔ2σA2+σD20Δ354σA2+34σADI+σDI2+ι*4+σD24ι2Δ434σA2+12σD20Δ554σA2+34σADI+σDI2+ι*4+σD24ι2Δ634σA2+12σD20Δ7σA2+σDI2+ι*8+σADI2+σD2414ι*Δ834σA2+σADI4+σDI2+ι*16+316σD214ιΔ9σA22+σD240

To see where these expressions come from, consider for example identity state Δ3, with, say, χl1:=χli[1],1=χli[1],2=χli[2],1≠χli[2],2=:χl2, where “=” here means identical by descent. Then, using Equations (21)–(23),

E[1M∑l=1MΦ(l)2|Δ3]=1M∑l=1ME[(3η(χl1)+η(χl2)2+2ϕ(χl1,χl1)+2ϕ(χl1,χl2)4)2]=54σA2+34σADI+14(σDI2+ι*)+14σD2.

The following quantities can be calculated in the same way. They are important for calculating the covariance between the trait values of parent and offspring (in particular the covariance between (Ai+Di) and Zi[1] and Zi[2]) which will dictate the change in distribution of the trait values within families arising from conditioning on knowing the traits of the parents. We record them here for later reference.

id. stateE[1M∑l=1MΦ(l)Ψl(i[1])|Δ⋅]E[1M∑l=1MΦ(l)Ψl(i[2])|Δ⋅]−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−Δ12σA2+2σADI+σDI2+ι*2σA2+2σADI+σDI2+ι*Δ2σA2+σADI2σA2+σADI2Δ332σA2+σDI2+ι*2+54σADIσA2+σADI4+σD22Δ4σA2+12σADIσA22Δ5σA2+σADI4+σD2232σA2+σDI2+ι*2+54σADIΔ6σA22σA2+12σADIΔ7σA2+σADI4+σD22σA2+σADI4+σD22Δ834σA2+18σADI+14σD234σA2+18σADI+14σD2Δ912σA212σA2

We can express two- and three-way identities between the parents in terms of the four-way identities Δ1,…,Δ9. Recall that we write, for example, F11 for the probability of identity of the two genes in i[1] and F12 for the probability of identity of two genes, one selected at random from i[1] and one from i[2]. In terms of the nine identity states, we have

F11=P[Δ1]+P[Δ2]+P[Δ3]+P[Δ4]F22=P[Δ1]+P[Δ2]+P[Δ5]+P[Δ6]F12=P[Δ1]+12(P[Δ3]+P[Δ5]+P[Δ7])+14P[Δ8]F112=P[Δ1]+12P[Δ3]F122=P[Δ1]+12P[Δ5]F1122=P[Δ1]F~1122=P[Δ2]F~1212=P[Δ7].

Combining the above, we find

E[1M∑l=1MΦ(l)Ψl(i[1])]=σA22(1+F11+2F12)+σADI2(F11+F12+2F122)+σD2(F12−F112)+(σDI2+ι*)F112,

with a symmetric expression for 1M∑l=1ME[Φ(l)Ψl(i[2])]. Similarly,

E[(Ai+Di)]=1M∑l=1ME[Φ(l)]=ιF12

and

1M∑l=1ME[Φ(l)2]=σA22(1+F11+F222+2F12)+σADI(F12+F112+F1222)+σDI2+ι*4(F12+F112+F122+F1122)+σD24(1−F12+F11−F112+F22−F122+F~1122+12F~1212)+14ι*F~1212,

from which, since for l≠m we are assuming E[Φ(l)Φ(m)]=E[Φ(l)]E[Φ(m)],

1ME[∑l=1M∑m=1Ml≠mΦ(l)Φ(m)]=(ιF12)2−ι*F122,

and the expression (6) for the variance of (Ai+Di) follows.

Remark B1 Walsh and Lynch (2018) give an expression for the variance when there is linkage disequilibrium. In their notation, f~ is the probability of identity at two distinct loci. Then for l≠m,E[Φ(l)Φ(m)]=f~E[Φl(χ^l,χ^l)]E[Φm(χ^m,χ^m)],

so that our expression for1ME[∑l=1M∑m=1Ml≠mΦ(l)Φ(m)]

will be multiplied by (f~/F122), resulting (when we subtract E[Ai+Di]2) in an overall expression of (f~−F122)ι2−f~ι* in place of −ι*F122. Correcting for this by adding (f~−F122)(ι2−ι*) to our expression (11) for the variance of Zi (for which we recall that F12 becomes Fii), we recover the expression of Walsh and Lynch (2018).

The covariance between Ai+Di and Aj+Dj

To understand the expression (7) for the covariance between Ai+Di and Aj+Dj for i≠j, consider

E[{ηl(χli[1],1)+ηl(χli[1],2)+ηl(χli[2],1)+ηl(χli[2],2)2+ϕl(χli[1],1,χli[2],1)+ϕl(χli[1],1,χli[2],2)4+ϕl(χli[1],2,χli[2],1)+ϕl(χli[1],2,χli[2],2)4}×{ηl(χlj[1],1)+ηl(χlj[1],2)+ηl(χlj[2],1)+ηl(χlj[2],2)2+ϕl(χlj[1],1,χlj[2],1)+ϕl(χlj[1],1,χlj[2],2)4+ϕl(χlj[1],2,χlj[2],1)+ϕl(χlj[1],2,χlj[2],2)4}].

The 16 terms corresponding to products of additive effects correspond to the 16 different possibilities for the allelic types at locus l if we choose one allele at random from individual i and one from individual j, and the contribution to the expectation will be nonzero precisely if the chosen alleles are identical, in which case they contribute E[ηl(χ^l)2]. Summing over l, the overall contribution of such terms to the covariance will therefore be 2σA2Fij.

Similarly, terms involving one factor of ηl and one ϕl will only be nonzero if all evaluated on the same allelic type, hence the terms multiplied by Fiij and Fijj in Equation (7).

Continuing in this way and using that E[Ai+Di]=ιFii, we recover Equation (7).

The residuals RAi+RDi

The corresponding calculations for the mean and variance of the residuals, RAi+RDi follow exactly the same pattern. It is convenient to consider RAi and RDi separately, and then calculate the covariance. The first of these, corresponding to the additive part is very straightforward since it is only going to depend on pairwise identities.

Recall first that

RAi=1M∑l=1M{(Xi−12)ηl(χli[1],1)+(12−Xli)ηl(χli[1],2)+(Yi−12)ηl(χli[2],1)+(12−Yi)ηl(χli[2],2)}.

Since the Mendelian inheritance is independent of the allelic states, RAi has mean zero; to establish the variance, we must calculate its square. Since inheritance is independent at distinct loci, only the diagonal terms contribute and we find

E[(RAi)2]=1M∑l=1ME[{(Xi−12)ηl(χli[1],1)+(12−Xli)ηl(χli[1],2)+(Yi−12)ηl(χli[2],1)+(12−Yi)ηl(χli[2],2)}2]=14M∑l=1ME[(ηl(χli[1],1))2+(ηl(χli[1],2))2+(ηl(χli[2],1))2+(ηl(χli[2],2))2]−12M∑l=1ME[ηl(χli[1],1)ηl(χli[1],2)+ηl(χli[2],1)ηl(χli[2],2)]=1M∑l=1MVar(ηl(χ^l))−12M∑l=1M(F11+F22)Var(η(χ^l))=(1−F11+F222)σA22.

This is, of course, exactly the expression we would obtain in the purely additive case.

The second residual, RDi, also has mean zero, but its variance will now involve higher order identities. Recall that

RDi=1M∑l=1M{(XliYli−14)ϕl(χli[1],1,χli[2],1)+(Xli(1−Yli)−14)ϕl(χli[1],1,χli[2],2)+((1−Xli)Yli−14)ϕl(χli[1],2,χli[2],1)+((1−Xli)(1−Yli)−14)ϕl(χli[1],2,χli[2],2)}.

Once again, since Mendelian inheritance is independent at different loci, E[(RDi)2] will be entirely determined by the diagonal terms. Note that for independent Bernoulli (parameter 1/2) random variables X and Y,

E[(XY−14)2]=E[12XY+116]=316,

and

E[(XY−14)(X(1−Y)−14)]=−116.

So, taking expectations over the variables Xli and Yli, we find

(B2) E[(RDi)2]=316M∑l=1ME[ϕl(χli[1],1,χli[2],1)2+ϕl(χli[1],1,χli[2],2)2+ϕl(χli[1],2,χli[2],1)2+ϕl(χli[1],2,χli[2],2)2]−216M∑l=1ME[ϕl(χli[1],1,χli[2],1)ϕl(χli[1],1,χli[2],2)+ϕl(χli[1],1,χli[2],1)ϕl(χli[1],2,χli[2],1)+ϕl(χli[1],1,χli[2],1)ϕl(χli[1],2,χli[2],2)+ϕl(χli[1],1,χli[2],2)ϕl(χli[1],2,χli[2],1)+ϕl(χli[1],1,χli[2],2)ϕl(χli[1],2,χli[2],2)+ϕl(χli[1],2,χli[2],1)ϕl(χli[1],2,χli[2],2)].

The first term depends only on pairwise identities and we see immediately that it is

34MF12∑l=1ME[ϕl(χ^l,χ^l)2]+34M(1−F12)∑l=1ME[ϕl(χ^l1,χ^l2)2]=34F12(σDI2+ι*)+34(1−F12)σD2.

The second term in Equation (B2) is most easily calculated conditional on identity class. Let us write Ξ(l) for the summand corresponding to locus l.

identity stateE[1M∑l=1MΞ(l)|Δ⋅]−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−Δ134(σDI2+ι*)Δ234σD2Δ318(σD2+σDI2+ι*)Δ414σD2Δ518(σD2+σDI2+ι*)Δ614σD2Δ718(σD2+ι*)Δ80Δ90

Using our notation for identities, this becomes

−14(F1122+F122+F112)(σDI2+ι*)−14ι*F~1212−14(F11−F112+F22−F122+F~1122+12F~1212)σD2.

Thus,

E[(RDi)2]=14(3F12−F1122−F122−F112)(σDI2+ι*)−14ι*F~1212+14(3(1−F12)−(F22−F122)−(F11−F112)−F~1122−12F~1212)σD2.

The covariance of RAi and RDi

Since RAi has mean zero, it suffices to calculate E[RAiRDi]. We need to establish the mean of

(B3) {(X−12)ηl(χli[1],1)+(12−X)ηl(χli[1],2)+(Y−12)ηl(χi[2],1)+(12−Y)ηl(χli[2],2)}×{XYϕl(χli[1],1,χli[2],1)+X(1−Y)ϕl(χli[1],1,χli[2],2)+(1−X)Yϕl(χli[1],2,χli[2],1)+(1−X)(1−Y)ϕl(χli[1],2,χli[2],2)}.

We have been able to drop the “−1/4” terms in the second bracket since E[RAi]=0.

Now

E[(X−12)XY]=12E[XY]=18,E[(X−12)(1−X)Y]=−12E[(1−X)Y]=−18,

and so the mean of Equation (B3) is that of

18{(ηl(χli[1],1)−ηl(χi[1],2))ϕl(χli[1],1,χli[2],1)+(ηl(χli[2],1)−ηl(χli[2],2))ϕl(χli[1],1,χli[2],1)+(ηl(χli[1],1)−ηl(χli[1],2))ϕl(χli[1],1,χli[2],2)+(ηl(χli[2],2)−ηl(χli[2],1))ϕl(χli[1],1,χli[2],2)+(ηl(χli[1],2)−ηl(χli[1],1)ϕl(χli[1],2,χli[2],1)+(ηl(χli[2],1)−ηl(χli[2],2))ϕl(χli[1],2,χli[2],1)+(ηl(χli[1],2)−ηl(χi[1],1))ϕl(χli[1],2,χli[2],2)+(ηl(χli[2],2)−ηl(χli[2],1))ϕl(χli[1],2,χli[2],2)}.

Taking expectations (conditional on the pedigree) and summing over loci, we find

E[RAiRDi|P(t)]=(F12−F112+F1222)σADI2.

Finally, for two distinct parents, we have found that in generation t, conditional on the pedigree up to time t,

Var(RAi+RDi)=(1−F11+F222)σA22+14(3F12−F112−F122−F1122)(σDI2+ι*)+14(3(1−F12)−(F11−F112)−(F22−F122)−F~1122−12F~1212)σD2+(F12−F112+F1222)σADI−14ι*F~1212.

We can also read off the result for when the two parents are the same from this formula. In that case

F1122=F11=F22=F112=F122,F12=12(1+F11),andF~1212=1−F11.

Thus, Var(RAi+RDi) reduces to

(1−F11)(σA22+38(σDI2+ι*)+14σD2+12σADI)−14ι*.

Appendix C: Conditioning multivariate Gaussian vectors

For ease of reference, we record here a standard result for conditioning multivariate normal random vectors on their marginal values.

Theorem C1 Suppose that[xAxB]∼N([μAμB],[ΣAAΣABΣBAΣBB]).

ThenxA|xB∼N(μA+ΣABΣBB−1(xB−μB),ΣAA−ΣABΣBB−1ΣBA).

The proof can be found, for example, in Brockwell (1996, Proposition 1.3.1 in Appendix A).

Appendix D: Generalized central limit theorems

We shall exploit known techniques for proving both convergence to a normal distribution, and for establishing the rate of convergence, in situations which go beyond the classical setting of independent identically distributed random variables. For convenience we recall the key results that we need here.

We begin with a result of Rinott (1994) on the rate of convergence in a generalized Central Limit Theorem; generalized because the summands are not identically distributed and it allows some dependence between elements in the sum. We do not use this second feature here, but it would be needed to extend our results to include effects that depend on more than one locus, and so for completeness we include it in the statement of the result. It also gives an idea of how quickly the rate of convergence deteriorates if one includes epistasis or higher order dominance effects. This result can be used both to prove asymptotic normality when we condition only on the pedigree (and not on any observed trait values), and to prove asymptotic normality of the residuals (that is the part of the trait distribution within families that is not shared among offspring) conditional on the observed traits of ancestors in the pedigree.

The dependence is captured by a dependency graph.

Definition D1 Let {Xl;l∈V} be a collection of random variables. The graph G=(V,E), where V and E denote the vertex set and edge set respectively, is said to be a dependency graph for the collection if for any pair of disjoint subsets A1 and A2 of V such that no edge in E has one endpoint in A1 and the other in A2, the sets of random variables {Xl;l∈A1} and {Xl;l∈A2} are independent.

The degree of a vertex in the graph is the number of edges connected to it and the maximal degree of the graph is just the maximum of the degrees of the vertices in it.

Theorem D2 Rinott 1994, Theorem 2.2)

Let E1,…,EM be random variables having a dependency graph whose maximal degree is strictly less than D, satisfying |El−E[El]|≤B a.s., l=1,…,M, E[∑l=1MEl]=λ and Var(∑l=1MEl)=σ2>0. Then, for every w∈R,|P[∑l=1MEl−λσ≤w]−N(w)|≤1σ{12πDB+16(Mσ2)1/2D3/2B2+10(Mσ2)D2B3},

where N is the distribution function of a standard normal random variable.

In particular, when D and B are order one and σ2 is of order M, the bound is of order 1/M.

Since we are only allowing for dominance effects that depend on allelic states at a single locus, and we have no epistasis, our dependency graphs will have no edges and so the maximal degree of any vertex will be zero and we may take D=1. Epistasis or higher order dominance effects, will increase the degree. This bound on the accuracy of the normal approximation will decrease rapidly as the number of combinations through which the allelic state at a single locus can influence the trait grows.

Exchangeable pairs

In order to prove the asymptotic normality of the part of the trait value that is shared by all the offspring in a family conditional on parental traits, we require a different approach. Because we are conditioning on the trait values of the parents, there will be weak dependence between all the pairs of loci within the sums defining Ai+Di (and so the dependency graph for the summands would be the complete graph). To check that nonetheless the limit is Gaussian we shall use a variant of Stein’s method of exchangeable pairs, originally introduced in Stein (1986).

Recall that the pair of random variables (W,W′) is called an exchangeable pair if their joint distribution is symmetric. Suppose that E[W]=0, E[W2]=1, (W,W′) is an exchangeable pair and

E[W−W′|W]=λ(W−R),

for some 0<λ<1, where R is a random variable of small order.

Let us write Δ=W−W′ and define

K^(t)=Δ2λ(1{−Δ≤t≤0}−1{0≤t≤−Δ}).

Note that ∫−∞∞K^(t)dt=Δ2/(2λ). In this case, one can show (see Chen et al. 2011, §2.3) that

(D1) E[Wf(W)]=E[∫−∞∞f′(W+t)K^(t)dt]+E[Rf(W)].

Proposition D3 Chen et al. 2011, Proposition 2.4i)

Let h be an absolutely continuous function with ‖h′‖<∞, and F any σ-algebra containing σ(W). If Equation (D1) holds, then|E[h(W)]−N(h)|≤‖h′‖(2πE[|1−K^1|]+2E[K^2]+2E[|R|]),

whereK^1=E[∫−∞∞K^(t)dt|F]=E[Δ22λ|F],andK^2=∫−∞∞|tK^(t)|dt=|Δ|34λ.

Corollary D4 Suppose that (W,W′) is an exchangeable pair with E[W]=μW and Var(W)=σW2 with(D2) E[W′|W]=(1−λ)W+λE[W]−λR

where R is a random variable of small order. Then defining K^1, K^2, h and F as in Proposoition D3,(D3) |E[h(W)]−NμW,σW2(h)|≤‖h′‖(2π1σWE[|σW2−K^1|]+2σW2E[K^2]+2E[|R|]),

where NμW,σW2 denotes the distribution of a normal random variable with mean μW and variance σW2.

Remark D5 Although this result is enough to guarantee that W is asymptotically normal, because we require ‖h′‖<∞, it is not enough to bound even the distance between the cumulative distribution function of W and that of a standard normal random variable with an error of order 1/M. To propagate our argument from one generation to the next requires convergence of the density function of the observed trait value, and once again it is our assumption that there is some environmental noise (with a smooth density) that allows us to guarantee this convergence based on the result proved here.

Appendix E: Key lemmas

Notation E1 Throughout the rest of the appendices, to ease the notation we shall assume that the (Gaussian) environmental noise is subsumed into the trait value Z, so that its distribution can be assumed to have a smooth density. That is, what we call Z below is the observed trait Z~ discussed in the main text. Moreover, when we write P[Z=z], we actually mean the density function of the distribution of Z~ evaluated at the value z (in formula, P[Z=z]:=φZ~(z) with φZ~ the density of Z~). This notation allows us to cover both the case when the allelic distributions are general (potentially concentrated on a finite number of values) and the environmental component is smooth enough that the distribution of their sum is also smooth, and the case when there is no environmental noise but the scaled allelic distributions have a smooth density over [-B,B] (in which case the distribution of the genetic component Z is itself smooth enough for the method below to be employed).

In this section, we prove two key lemmas which will underpin our proof. They will allow us to estimate the effect on the distribution of the allelic types at a particular locus, or particular pair of loci, of knowing the trait value. We shall be using Bayes’ rule. With a slight abuse of notation

P[(χl1,χl2)=(x,x′)|Z=z]=P[Z=z|(χl1,χl2)=(x,x′)]P[Z=z]P[(χl1,χl2)=(x,x′)].

Let us write Ψl(x,x′)=ηl(x)+ηl(x′)+ϕl(x,x′) and Z−l for the trait value of an individual with the effect of locus l removed, then the ratio in this expression becomes

P[Z−l=z−Ψl(x,x′)]P[Z=z].

Of course, this ratio of probabilities should be interpreted as a ratio of density functions. Moreover, bearing in mind our remarks on environmental noise, we are going to suppose that these density functions are sufficiently smooth that we can justify an application of Taylor’s Theorem. Of course, we know that Z−l is approximately normally distributed, using exactly the same argument as for Z, and it is no surprise that the ratio differs from one by something of order 1/M. The importance of the next lemma will become evident when we sum conditional expectations over loci; c.f. Remark E5.

Lemma E2 In the notation above,P[Z−l=z]=P[Z=z]+1ME[Ψl(χl1,χl2)]ddzP[Z=z]+1ME[Ψl(χl1,χl2)]2d2dz2P[Z=z]−12ME[Ψl(χl1,χl2)2]d2dz2P[Z=z]+Cl(z)1M3/2,

where the function Cl(z) in the error term can be bounded independent of l and z.

Remark E3 Conditioning on the pedigree)

Although we have suppressed it in the notation, this lemma holds in any generation, but the expressions E[Ψ(χl1,χl2)]2 and E[Ψ(χl1,χl2)2] should be interpreted as being calculated conditional on the pedigree (which will determine the probability of identity of χl1, χl2).

Proof of Lemma E2 We are going to abuse notation (still further) and imagine that P[χl1=x,χl2=x′,Z=z] has a density with respect to x, x′. Of course, we do not expect that to be true (even with environmental noise), but it makes our expressions easier to parse than using a more mathematically accurate notation. We begin with an application of Taylor’s Theorem (with respect to z):(E1) P[Z−l=z]=∫∫P[χl1=x,χl2=x′,Z=z+1MΨl(x,x′)]dxdx′

(E2) =∫∫P[χl1=x,χl2=x′,Z=z]dxdx′

(E3) +1M∫∫Ψl(x,x′)∂∂zP[χl1=x,χl2=x′,Z=z]dxdx′+12M∫∫Ψl(x,x′)2∂2∂z2P[χl1=x,χl2=x′,Z=z]dxdx′

(E4) +C^l(z)1M3/2.

Provided that P[Z=z] has a uniformly bounded third derivative, our assumption that the terms that make up Ψl are uniformly bounded allows us to deduce that C^l is uniformly bounded in l and z. Notice that the expression in Equation (E2) is just P[Z=z].

Since we are not conditioning on any trait values in the pedigree, and the ancestral population is assumed to be in linkage equilibrium, (χl1,χl2) and Z−l are independent. Combining this observation with Equation (E1), and, once again applying Taylor’s Theorem, we findP[χl1=x,χl2=x′,Z=z]=P[χl1=x,χl2=x′]P[Z−l=z−1MΨl(x,x′)]=P[χl1=x,χl2=x′]∫∫P[χl1=y,χl2=y′,1MZ=z−1MΨl(x,x′)+1MΨl(y,y′)]dydy′=P[χl1=x,χl2=x′]{P[Z=z]1M+1M∫∫(Ψ(y,y′)−Ψ(x,x′))×∂∂zP[χl1=y,χl2=y′,Z=z]dydy′+C~l(x,x′,z)1M},

where the function C~l in the last line is uniformly bounded independent of l and (x,x′,z). (To justify this last statement, recall that we are abusing notation and implicitly subsuming the environmental noise into the distribution of Z. The density function here is actually a convolution of that of the environmental noise, which is smooth, and the true distribution of Z, and is therefore smooth.) Still assuming sufficient regularity, differentiating the previous equation we find(E5) ∂∂zP[χl1=x,χl2=x′,Z=z]=P[χl1=x,χl2=x′]{ddzP[Z=z]+1M∫∫(Ψ(y,y′)−Ψ(x,x′))∂2∂z2P[χl1=y,χl2=y′,Z=z]dydy′+∂∂zC~l(x,x′,z)1M},

and(E6) ∂2∂z2P[χl1=x,χl2=x′,Z=z]=P[χl1=x,χl2=x′]{d2dz2P[Z=z]+1M∂2∂z2C~~l(x,x′,z)},

with ∂∂zC~l and ∂2∂z2C~~l uniformly bounded.

Finally, substituting Equations (E5) and (E6) in Equations (E3) and (E4), we obtainP[Z−l=z]=P[Z=z]+1MddzP[Z=z]∫∫Ψl(x,x′)P[χl1=x,χl2=x′]dxdx′+1Md2dz2P[Z=z]∫∫∫∫(Ψl(y,y′)−Ψl(x,x′))Ψl(x,x′)×P[χl1=y,χl2=y′]P[χl1=x,χl2=x′]dydy′dxdx′+12Md2dz2P[Z=z]∫∫Ψl(x,x′)2P[χl1=x,χl2=x′]dxdx′+C^l(z)1M3/2=P[Z=z]+1ME[Ψl(χl1,χl2)]ddzP[Z=z]+1ME[Ψl(χl1,χl2)]2d2dz2P[Z=z]−12ME[Ψl(χl1,χl2)2]d2dz2P[Z=z]+C^l(z)1M3/2,

as required.  □

We also require an analog of Lemma E2 with which to control the effect of conditioning on the trait value on the distribution of the allelic values at pairs of loci. We write Z−l−m=Z−1M(Ψl(χl1,χl2)+Ψm(χm1,χm2)) for the trait value with the contributions from loci l and m removed. The following lemma follows on iterating the argument that gave us Lemma E2.

Lemma E4 In the notation above,P[Z−l−m=z]=P[Z=z]+1M(E[Ψl(χl1,χl2)]+E[Ψm(χm1,χm2)])ddzP[Z=z]+{1ME[Ψl(χl1,χl2)]2−12ME[Ψl(χl1,χl2)2]+1ME[Ψl(χl1,χl2)]E[Ψm(χm1,χm2)]+1ME[Ψm(χm1,χm2)]2−12ME[Ψm(χm1,χm2)2]}d2dz2P[Z=z]+Cl,m(z)1M3/2,

where the functions Cl,m(z) are uniformly bounded in l, m, z.

Proof of Lemma E4 We iterate the previous result:P[Z−l−m=z]=P[Z−l=z]+1ME[Ψm(χm1,χm2)]ddzP[Z−l=z]+1ME[Ψm(χm1,χm2)]2d2dz2P[Z−l=z]−12ME[Ψm(χm1,χm2)2]d2dz2P[Z−l=z]+Cm(z)1M3/2;

now substitute for P[Z−l=z] and its derivatives.  □

Remark E5 Just as for Lemma E2, the proof of Lemma E4 applies in any generation as long as one interprets the expectations as being taken conditional on the pedigree. We have assumed that our base population is in linkage equilibrium to write E[Ψl(y,y′)Ψm(x,x′)]=E[Ψl(y,y′)]E[Ψm(x,x′)].

We shall only be presenting the detailed proofs for individuals in generation one. To extend to the general case requires an analog of Lemma E2 when we consider the trait values of the two parents of an individual. For completeness, we record that lemma here.

Lemma E6 Let us use P[z1,z2] to denote P[Zi[1]=z1,Zi[2]=z2]. In the following expression, all expectations should be interpreted as taken conditional on the pedigree:P[Z−li[1]=z1,Z−li[2]=z2]−P[z1,z2]=1ME[Ψl(χli[1],1,χli[1],2)]∂∂z1P[z1,z2]+1ME[Ψl(χli[2],1,χli[2],2)]∂∂z2P[z1,z2]+(1ME[Ψl(χli[1],1,χli[1],2)]2−12ME[Ψl(χli[1],1,χli[1],2)2])×∂2∂z12P[z1,z2]

+(1ME[Ψl(χli[2],1,χli[2],2)]2−12ME[Ψl(χli[2],1,χli[2],2)2])×∂2∂z22P[z1,z2]+(2ME[Ψl(χli[1],1,χli[1],2)]E[Ψl(χli[2],1,χli[2],2)]−1ME[Ψl(χli[1],1,χli[1],2)Ψl(χli[2],1,χli[2],2)])×∂2∂z1∂z2P[z1,z2]+O(1M3/2).

Appendix F: Mean and variance of trait values conditional on parental traits

We remind the reader that Notation E1 remains in force.

We now turn to calculating the conditional distribution of the trait values, conditional not just on the pedigree, as we did in Appendix B, but also on the (observed) trait values in the parental generation. We spell out the details in generation one. Here already we can identify the key points, without being overwhelmed by notation. Recall that we are implicitly conditioning not on the exact trait values of the parents, but on the observed trait values when environmental noise is taken into account, so that we can assume that the distribution of parental trait values has a smooth density.

First, we calculate the conditional mean. We distinguish the case of two distinct parents and a family produced by selfing. Recall that we wrote Ai+Di for the component shared by all individuals in the family, with Ai and Di defined in Equations (28) and (29).

 

Generation one: mean trait value, distinct parents

Since the parents are, by assumption, unrelated, we anticipate that the expected value of the dominance component is zero, and so the expected value of the shared component Ai+Di should be the mean value of the parental traits. However, since we are conditioning on knowing the trait values, we do have some information about the allelic types, and we must verify that this does not significantly distort the expectations.

We exploit again the fact that since the parents are unrelated, their trait values (and allelic states at locus l) are independent. Thus,

(F1) P[Zi[1]=z1,Zi[2]=z2|(χli[1],1,χli[1],2χli[2],1χli[2],2)=(x,x′,y,y′)]=P[Z−li[1]=z1−1M(ηl(x)+ηl(x′)+ϕl(x,x′))]×P[Z−li[2]=z2−1M(ηl(y)+ηl(y′)+ϕl(y,y′))].

We now use Lemma E2 and Taylor’s Theorem to deduce that

P[(χli[1],1,χli[1],2)=(x,x′)|Zi[1]=z]=P[(χli[1],1,χli[1],2)=(x,x′)]×{1−1M(Ψl(x,x′)−E[Ψl])1P[Zi[1]=z]ddzP[Zi[1]=z]+O(1M)},

with a symmetric expression for i[2]. Integrating against this expression and using Equations (21), (22), and (23), we find, in an obvious notation,

(F2) E[ηl(χli[1],1)|i[1]≠i[2],Zi[1],Zi[2]]=−1MP′[Zi[1]]P[Zi[1]]Var(ηl(χ^l))+O(1M).

Note that approximating P[Zi[1]] by a normal density and ignoring the environmental component, the order 1/M terms involves 1/(σA2+σD2) and (Zi[1]−z¯0)2/(σA2+σD2), and is controlled through these quantities and our bounds on ηl and ϕl. In particular, the approximation breaks down if the genetic variance is too small or if the trait of the parent is too extreme. Multiplying by 1/M and summing over loci and parents, we arrive at

(F3) E[Ai|i[1]≠i[2],Zi[1],Zi[2]]=−(P′[Zi[1]]P[Zi[1]]+P′[Zi[2]]P[Zi[2]])1M∑l=1MVar(ηl(χ^l))+O(1M).

Remark F1 Since we already checked that the trait Zi[1] is approximately normally distributed, and the same argument evidently gives that Z−li[1] is approximately normally distributed for each l, the derivation above may seem unnecessarily complex. However, in summing the terms in Equation (F2) over loci, we exploited the fact that we could pull the ratio P′[Zi[1]]/P[Zi[1]] outside the sum. Only then did we approximate it by the limiting normal distribution. We could only do this because we expressed everything in terms of the distribution of the whole trait. If we try to approximate the distribution of Z−li[1] directly by a normal distribution, and then sum, we cannot control the error. We shall use this trick repeatedly in what follows.

Similarly,

E[ϕl(χli[1],1,χli[2],1)|i[1]≠i[2],Zi[1]=z1,Zi[2]=z2]=∫…∫ϕl(x,y){1−1M(Ψl(x,x′)−E[Ψl])1P[Zi[1]=z1]×ddz1P[Zi[1]=z1]}×{1−1M(Ψl(y,y′)−E[Ψl])1P[Zi[2]=z2]ddz2P[Zi[2]=z2]}×ν^l(dx)ν^l(dx′)ν^l(dy)ν^l(dy′)+O(1M).

The terms of order one and 1/M vanish as a result of Equations (20), (21), (22), and (23). Multiplying by 1/M and summing over loci, we find that E[Di]=O(1/M).

Recalling that the trait distribution in the ancestral population is (almost) normally distributed with mean z¯0, we see that if we ignore environmental effects, so that the variance of the trait distribution in generation zero is σA2+σD2, then adding z¯0 to the right-hand side of Equation (F3), and substituting

P[Zi[1]]=12π(σA2+σD2)exp(−(Zi[1]−z¯0)22(σA2+σD2)),

we recover that up to an error of order 1/M, the expected trait value among offspring is

z¯0+σA2σA2+σD2(Zi[1]+Zi[2]2−z¯0),

as predicted by Theorem C1.

Remark F2 The breeder’s equation)

Suppose that as a result of environmental noise, the observed trait of each individual in the ancestral population is its genetic trait plus an independent N(0,σE2) random variable. Then assuming normality of the ancestral trait distribution, and using Theorem C1, we find that for unrelated parents the mean trait in generation one is(F4) z¯0+σA2σZ2((Zi[1]+Zi[2])2−z¯0),

where σZ2 is the total variance of the observed trait in the ancestral population; that is σZ2=σA2+σD2+σE2. Equation (F4) is the breeder’s equation.

Mean trait value, same parent

We now turn to the expected trait value in a family in generation one that is produced by selfing. The calculation for the additive term is unchanged, but now we have a nontrivial contribution from the dominance component. We denote the parent Zi[1]. Since Zi[1]=Zi[2], we must calculate E[ϕl(χli[1],1,χli[1],1)|Zi[1]] and E[ϕl(χli[1],1,χli[1],2)|Zi[1]].

Our strategy is as before: we express each of these probabilities in terms of the distribution of the trait value minus the contribution from locus l and we apply Lemma E2. Thus, once again using that in generation zero, before conditioning, the two alleles at locus l in Zi[1] are independent draws from ν^l,

E[ϕl(χli[1],1,χli[1],1)|i[1]=i[2],Zi[1]=z]=∫∫ϕl(x,x)(1−1M(Ψl(x,x′)−E[Ψl])1P[Zi[1]=z]×ddzP[Zi[1]=z])ν^l(dx)ν^l(dx′)+O(1M).

Using Equations (21), (22), (23), we see that on integration the only nonzero contribution comes from the term ηl(x)ϕl(x,x) which can be integrated to yield

E[ϕl(χli[1],1,χli[1],1)|i[1]=i[2],Zi[1]]=E[ϕl(χ^l,χ^l)]−1MP′[Zi[1]]P[Zi[1]]E[ηl(χ^l)ϕl(χ^l,χ^l)]+O(1M).

Similarly,

E[ϕl(χli[1],1,χli[1],2)|i[1]=i[2],Zi[1]]=−1MP′[Zi[1]]P[Zi[1]]E[ϕl(χ^l,χ^2)2]+O(1M).

Multiplying by 1/M and summing over loci, we find that the mean of the term Di in Equation (29), conditional on i[1]=i[2] and on knowing the trait value Zi[1], is

12M∑l=1ME[ϕl(χ^l,χ^l)]−P′[Zi[1]]P[Zi[1]](12M∑l=1ME[ηl(χ^l)ϕl(χ^l,χ^l)]+12M∑l=1ME[ϕl(χ^l,χ^2)2])+O(1M).

Adding on the additive terms that we calculated before and restating everything in terms of the quantities in Table 1, we obtain that for two identical parents

(F5) E[Zi|i[1]=i[2],Zi[1]]=z¯0+(12ι−P′[Zi[1]]P[Zi[1]](σA2+σD22+σADI4))+O(1M).

Notice that the factor of 1/2 in front of ι is the probability of identity F* of the two genes in the offspring.

Of course, there is no surprise here: E[(Ai+Di)|i[1]=i[2]]=ι/2 and

(F6) Cov(Ai+Di,Zi[1]|i[1]=i[2])=1M∑l=1ME[(ηl(χ^l1)+ηl(χ^21)+ϕ(χ^l1,χ^l1)+ϕ(χ^l2,χ^l2)4+ϕ(χ^l1,χ^l2)2)(η(χ^l1)+η(χ^l2)+ϕ(χ^l1,χ^l2))]=σA2+σD22+σADI4.

Thus, up to the error term, Equation (F5) is just

z¯0+E[Ai+Di]+Cov(Ai+Di,Zi[1])(Zi[1]−E[Zi[1]])Var(Zi[1]),

as we expect from the (approximately) bivariate normal distribution of (Ai+Di) and Zi[1].

Variance of the shared parental contribution Ai+Di, generation one

We now turn to the variance of the shared parental contribution. This is where the complications associated with incorporating dominance really start to be felt. In the process of calculating the conditional mean above, we established that conditioning on the parental trait values (and whether or not they are identical) distorts the distribution of the allelic state at a given locus by a factor of order 1/M. This distortion is enough to shift the mean trait (as we see in the breeder’s equation), and, as we shall see, the variance of the sum over loci will have a contribution from linkage disequilibrium.

Conditional variance (Ai+Di), generation one, same parent

First we consider the case in which the parents are the same. We need to calculate the expectation of (Ai+Di)2 conditional upon the parental trait. We begin with the “diagonal” terms, corresponding to a single locus. We take these in three parts. First, proceeding as before,

(F7) E[ηl(χli[1],1)2|i[1]=i[2],Zi[1]=z]=∫∫ηl(x)2(1−1M(Ψl(x,x′)−E[Ψl])×1P[Zi[1]=z]ddzP[Zi[1]=z])ν^l(dx)ν^l(dx′)+O(1M)=E[ηl(χ^l)2]+O(1M).

Notice that the term arising from the Taylor expansion is already of order 1/M, and, since we multiply each of the terms in the sum by 1/M, we have no need to develop the expansion further. Indeed, all terms in the expression for the variance will be multiplied by 1/M and so for the “diagonal” terms in the square of the sum, we only need an expression to leading order.

Remark F3 The error that we are making in discarding the terms arising from the Taylor expansion is 1/M multiplied by a term that depends on P′[Zi[1]]/P[Zi[1]]=−(Zi[1]−E[Zi[1]])/Var(Zi[1]). As usual, the approximation will be poor if the trait value of the parent is too extreme, or the variance is too small.

As a result, for these terms we can calculate with respect to the distribution in the ancestral population and we find

1ME[∑l=1M(ηl(χli[1],1)+ηl(χli[1],2))2|i[1]=i[2],Zi[1]]=2M∑l=1ME[ηl(χ^l)2]+O(1M)=σA2+O(1M).

Similarly, recalling that we are still considering the case of identical parents,

116ME[∑l=1M(ϕl(χli[1],1,χli[1],1)+2ϕl(χli[1],1,χli[1],2)+ϕl(χli[1],2,χli[1],2))2|i[1]=i[2],Zi[1]∑l=1M∑l=1M]=116M∑l=1M(2E[ϕl(χ^l,χ^l)2]+4E[ϕl(χ^l1,χ^l2)2])+O(1M)=18(σDI2+ι*)+14σD2+O(1M),

and

12ME[∑l=1M(ηl(χli[1],1)+ηl(χli[1],2))×(ϕl(χli[1],1,χli[1],1)+2ϕl(χli[1],1,χli[1],2)+ϕl(χli[1],2,χli[1],2))|i[1]=i[2],Zi[1]∑l=1M∑l=1M]=12ME[∑l=1M(ηl(χli[1],1)+ηl(χli[1],2))×(ϕl(χli[1],1,χli[1],1)+2ϕl(χli[1],1,χli[1],2)+ϕl(χli[1],2,χli[1],2))∑l=1M]+O(1M)

=12M∑l=1M2E[ηl(χ^l)ϕl(χ^l,χ^l)]+O(1M)=σADI2+O(1M).

Combining all these terms we find that if the parents are identical, then the contribution to E[(Ai+Di)2|i[1]=i[2],Zi[1]] from the “diagonal” terms is

(F8) σA2+18(σDI2+ι*)+14σD2+12σADI+O(1M).

We must now turn to the contribution from correlations across loci. For this, we must compute

(F9) 1ME[∑l≠m(ηl(χli[1],1)+ηl(χli[1],2))(ηm(χmi[1],1)+ηm(χmi[1],2))|i[1]=i[2],Zi[1]∑l≠m]

(F10) (F10)+12ME[∑l≠m(ηl(χli[1],1)+ηl(χli[1],2))(ϕm(χmi[1],1,χmi[1],1)+ϕm(χmi[1],2,χmi[1],2)+2ϕm(χmi[1],1,χmi[1],2))|i[1]=i[2],Zi[1]∑l≠m](F11)+116ME[∑l≠m(ϕl(χli[1],1,χli[1],1)+ϕl(χli[1],2,χli[1],2)+2ϕl(χli[1],1,χli[1],2))(ϕm(χmi[1],1,χmi[1],1)+ϕm(χmi[1],2,χmi[1],2)+2ϕm(χmi[1],1,χmi[1],2))|i[1]=i[2],Zi[1]∑l≠m].

This time we use Lemma E4.

(F12) 1P[Zi[1]=z]P[Zi[1]=z|(χli[1],1,χli[1],2,χmi[1],1,χmi[1],2)=(x,x′,y,y′)]=1P[Zi[1]=z]P[Z−l−mi[1]=z−1M(ηl(x)+ηl(x′)+ϕl(x,x′)+ηm(y)+ηm(y′)+ϕm(y,y′))1M]=1−1M(Ψl(x,x′)+Ψm(y,y′)−E[Ψl+Ψm])×1P[Zi[1]=z]ddzP[Zi[1]=z]+12M((Ψl(x,x′)+Ψm(y,y′))2−E[(Ψl+Ψm)2])×1P[Zi[1]=z]d2dz2P[Zi[1]=z]−1M(Ψl(x,x′)+Ψm(y,y′)−E[(Ψl+Ψm)])E[Ψl+Ψm]×1P[Zi[1]=z]d2dz2P[Zi[1]=z]+O(1M3/2).

Using that in the ancestral population we are at linkage equilibrium with x,x′ and y,y′ sampled independently from ν^l and ν^m, respectively, multiplying by ηl(x)ηm(y) and integrating against ν^(dx)ν^(dy), the only nonzero term corresponds to the term ηl(x)ηm(y) in (Ψl(x,x′)+Ψm(y,y′))2, so that

(F13) 1ME[∑l≠m(ηl(χli[1],1)+ηl(χli[1],2))×(ηm(χmi[1],1)+ηm(χmi[1],2))|i[1]=i[2],Zi[1]∑l≠m]=4M2P″[Zi[1]]P[Zi[1]]∑l≠mE[ηl(χ^l)2]E[ηm(χ^m)2]+O(1M)=P″[Zi[1]]P[Zi[1]](σA2)2+O(1M).

(The factor of 4 corresponds to the four possible ways of choosing the parents at the two loci.) Similarly, to calculate

E[ηl(χli[1],1)ϕm(χmi[1],1,χmi[1],1)|i[1]=i[2],Zi[1]]

we multiply Equation (F12) by ηl(x)ϕm(y,y) and integrate. Once again, using Equations (21)–(23), we find that most of the terms vanish, leaving only

(F14) −1MP′[Zi[1]]P[Zi[1]]E[ηl(χ^l)2]E[ϕm(χ^m,χ^m)]+12MP″[Zi[1]]P[Zi[1]]{E[ηl(χ^l)3]E[ϕm(χ^m,χ^m)]+E[ηl(χ^l)ϕl(χ^l1,χ^l2)2]E[ϕm(χ^m,χ^m)]+2E[ηl(χ^l)2]E[ηm(χ^m)ϕm(χ^m,χ^m)]}.

Multiplying by 1/(2M) and summing over loci, in the notation of Table 1, the first term yields

−ισA2P′[Zi[1]]P[Zi[1]].

[There are four terms of this form in Equation (F10) and we have taken account of all of them.] The last term gives

σA2σADI2P″[Zi[1]]P[Zi[1]]

[again counting the contribution from all four terms of this form in Equation (F10)].

Now observe that

1M∑m=1ME[ϕm(χ^m,χ^m)]=ιM,

so that summing over loci, the contribution from the first two terms multiplying the second derivative will be O(1/M).

Remark F4 Up to this point, it has been possible to neglect the error terms under the assumption that the within-family variance is not too small and we are not too far out into the tails of the distribution of Zi[1]; the more extreme the trait of the parent, the worse the approximation will be. Now things change. In order for E[Ai+Di] to be finite, we required that the inbreeding depression ι be well defined; here, we see that it also enters into the error terms.

In the same way, we calculate

E[ηl(χli[1],1)ϕm(χmi[1],1,χmi[1],2)|i[1]=i[2],Zi[1]]

by multipling Equation (F12) by ηl(x)ϕm(y,y′) and integrating. The only term to survive integration is

(F15) 12MP″[Zi[1]]P[Zi[1]]E[2ηl(χ^l)2]E[ϕm(χ^m1,χ^m2)2].

There are four terms of this form in Equation (F10), each of which is weighted by 1/(2M) and so, summing over loci, we arrive at an overall contribution of σA2σD2P″[Zi[1]]/P[Zi[1]]. Equations (F14) and (F15) yield that Equation (F10) equals

(F16) 12ME[∑l≠m(ηl(χli[1],1)+ηl(χli[1],2))(ϕm(χmi[1],1,χmi[1],1)+ϕm(χmi[1],2,χmi[1],2)+2ϕm(χmi[1],1,χmi[1],2))|i[1]=i[2],Zi[1]∑l≠m]=−P′[Zi[1]]P[Zi[1]]ισA2+P″[Zi[1]]P[Zi[1]]{σA2σADI2+σA2σD2}+O(1M).

Continuing in this way,

E[ϕl(χli[1],1,χli[1],1)ϕm(χmi[1],1,χmi[1],1)|i[1]=i[2],Zi[1]]

is obtained by multiplying Equation (F12) by ϕl(x,x)ϕm(y,y) and integrating. When we sum the “constant” term over loci we will obtain ι2/M which tends to zero. The remaining nonzero terms are

(F17) −1MP′[Zi[1]]P[Zi[1]]{E[ηl(χ^l)ϕl(χ^l,χ^l)]E[ϕm(χ^m,χ^m)]+E[ϕl(χ^l,χ^l)]E[ηm(χ^m)ϕm(χ^m,χ^m)]}+12MP″[Zi[1]]P[Zi[1]]{E[ϕl(χ^l,χ^l)ϕm(χ^m,χ^m)(ηl(χ^l)2+ηm(χ^m)2+2ηl(χ^l)ηm(χ^m))]+E[ϕl(χ^l1,χ^l1)ϕm(χ^m1,χ^m1)(ηl(χ^l2)2+ηm(χ^m2)2)]}.

The terms in the last line will contribute O(1/M) when we sum, as will the first two terms in the middle line. There are four terms of this form in Equation (F11) and we are multiplying by 1/(16M) and summing over loci, so the top line contributes −ισADIP′[Zi[1]]/(4P[Zi[1]]), similarly the second line will contribute σADI2P″[Zi[1]]/(16P[Zi[1]]).

Now, again using Equation (F12),

E[ϕl(χli[1],1,χli[1],1)ϕm(χmi[1],1,χmi[1],2)|i[1]=i[2],Zi[1]]=∫…∫ϕl(x,x)ϕm(y,y′)[1−1M(Ψl(x,x′)+Ψm(y,y′)−E[Ψl+Ψm])P′[Zi[1]]P[Zi[1]]+12M((Ψl(x,x′)+Ψm(y,y′))2−E[(Ψl+Ψm)2])P″[Zi[1]]P[Zi[1]]−1M(Ψl(x,x′)+Ψm(y,y′)−E[(Ψl+Ψm)])E[Ψl+Ψm]×P″[Zi[1]]P[Zi[1]]]ν^l(dx)ν^l(dx′)ν^m(dy)ν^m(dy′)+O(1M3/2).

There are eight terms of this form in Equation (F11), and we are multiplying by 1/(16M) and summing over loci, so the first term will correspond to a contribution of

−P′[Zi[1]]P[Zi[1]]ισD22.

As usual, terms multiplying the second derivative that involve the locus l only through ϕl(x,x) will contribute O(1/M) to the sum and we find that the nontrivial contributions will be

−1MP′[Zi[1]]P[Zi[1]]E[ϕl(χ^l,χ^l)]E[ϕm(χ^m1,χ^m2)2]+12MP″[Zi[1]]P[Zi[1]]E[2ηl(χ^l)ϕl(χ^l,χ^l)]E[ϕm(χ^m1,χ^m2)2].

There are eight terms of this form in Equation (F11), so multiplying by 1/(16M) and summing over loci gives

(F18) −ισD22P′[Zi[1]]P[Zi[1]]+σADIσD24P″[Zi[1]]P[Zi[1]].

Finally, when we scale and sum over loci, the only nontrivial term in our expression for

E[ϕl(χli[1],1,χli[1],2)ϕm(χmi[1],1,χmi[1],2)|i[1]=i[2],Zi[1]]

is

+12MP″[Zi[1]]P[Zi[1]]E[2ϕl(χ^l1,χ^l2)2]E[ϕm(χ^m1,χ^m2)2].

There are four terms of this form, and so multiplying by 1/(16M) and summing gives

(F19) (σD2)24P″[Zi[1]]P[Zi[1]].

Combining Equations (F17), (F18), and (F19), we find that Equation (F11) is

(F20) 116ME[∑l≠m(ϕl(χli[1],1,χli[1],1)+(ϕl(χli[1],2,χli[1],2)+2ϕl(χli[1],1,χli[1],2))(ϕm(χmi[1],1,χmi[1],1)+(ϕm(χmi[1],2,χmi[1],2)+2ϕm(χmi[1],1,χmi[1],2))|i[1]=i[2],Zi[1]∑l≠m]=−P′[Zi[1],1]P[Zi[1],1](ισADI4+ισD22)+P″[Zi[1],1]P[Zi[1],1](σADI216+σADIσD24+(σD2)24)+O(1M).

Adding Equations (F8), (F13), (F16), and (F20) yields E[(Ai+Di)2], and subtracting the square of Equation (F5), we obtain

Var(Ai+Di|i[1]=i[2],Zi[1])=σA2+18(σDI2+ι*)+14σD2+12σADI−P′[Zi[1]]P[Zi[1]]{ισA2+ισADI4+ισD22}+P″[Zi[1]]P[Zi[1]]{(σA2)2+σA2σADI2+σA2σD2+(σD2)24+σADI216+σD2σADI4}−(ι2−P′[Zi[1]]P[Zi[1]](σA2+σD22+σADI4))2+O(1M).

Now if we substitute the Gaussian density for Zi[1], observing that

P″[Zi[1]]P[Zi[1]]−(P′[Zi[1]]P[Zi[1]])2=−1σA2+σD2,

we see that the variance reduces to

(F21) −ι24+1σA2+σD2(σA2+σD22+σADI4)2+σA2+18(σDI2+ι*)+14σD2+12σADI+O(1M).

Again, that was a lot of work to recover exactly the expression that we expected from conditioning the multivariate normal random variable ((Ai+Di),Zi[1]) on its second argument. However, in the process, we have identified where the normal approximation to the conditioned process will break down. The bounds that we have obtained will be poor if the trait value of either parent is too extreme, or if the pedigree is too inbred (as a result of which the variance of trait values will be small and inbreeding depression may be high).

Of course, we have not proved that the conditional distribution of (Ai+Di) converges to a normal, we have just checked that the first two moments are asymptotically what we would expect. We defer the proof of normality until we have calculated the conditional variance of (Ai+Di) in the (much simpler) case of two distinct parents.

Conditional variance (Ai+Di), generation one, distinct parents

If the parents are distinct, then the expressions are much simpler. First

14ME[∑l=1M(ηl(χli[1],1)+ηl(χli[1],2)+ηl(χli[2],1)+ηl(χli[2],2))2|i[1]≠i[2],Zi[1],Zi[2]∑l=1M]=1M∑l=1ME[ηl(χ^l)2]+O(1M)=σA22+O(1M).

Next

116ME[∑l=1M(ϕl(χli[1],1,χli[2],1)+ϕl(χli[1],1,χli[2],2)+ϕl(χli[1],2,χli[2],1)+ϕl(χli[1],2,χli[2],2))2|i[1]≠i[2],Zi[1],Zi[2]∑l=1M]=14M∑l=1ME[ϕl(χ^l1,χ^l2)2]+O(1M)=14σD2+O(1M).

Finally,

14ME[∑l=1M(ηl(χli[1],1)+ηl(χli1,2)+ηl(χli[2],1)+ηl(χli[2],2))×(ϕl(χli[1],1,χli[2],1)+ϕl(χli[1],1,χli[2],2)+ϕl(χli[1],2,χli[2],1)ϕl(χli[1],2,χli[2],2))|i[1]≠i[2],Zi[1],Zi[2]∑l=1M]=O(1M).

We now turn to the off-diagonal terms. We need to be able to calculate the conditional expectation of

(F22) [12(ηl(χli[1],1)+ηl(χli[1],2)+ηl(χli[2],1)+ηl(χli[2],2))+14(ϕl(χli[1],1,χli[2],1)+ϕl(χli[1],1,χli[2],2)+ϕl(χli[1],2,χli[2],1)+ϕl(χli[1],2,χli[2],2))12]×[12(ηm(χmi[1],1)+ηm(χmi[1],2)+ηm(χmi[2],1)+ηm(χmi[2],2))+14(ϕm(χmi[1],1,χmi[2],1)+ϕm(χmi[1],1,χmi[2],2)+ϕm(χmi[1],2,χmi[2],1)+ϕm(χmi[1],2,χmi[2],2))12]

given the trait values in the (unrelated) parents i[1] and i[2]. Because the parents are distinct, and they are in generation zero, as in Equation (F1) in our calculation of the conditional mean, we can exploit the fact that the trait values Zi[1] and Zi[2] are independent so that the joint probability that

(χli[1],1,χli[1],2,χmi[1],1,χmi[1],2)=(x,x′,y,y′),

conditional on Zi[1],Zi[2], is just the same as if we only condition on Zi[1]. Recalling Equation (F1), we can calculate the conditional expectation of Equation (F22) using Equation (F12). None of the genes at either locus are identical by descent, and so integrating against the term of order 1/M in the Taylor expansion in Equation (F12) gives zero, but since we are calculating the conditional expectation of O(M2) terms, each of which is of order 1/M, we can expect to see a contribution from the term of order 1/M. All the terms involving the dominance components vanish, as do those terms involving only one copy of the additive component at one of the loci. In total we find that the conditional expectation of Equation (F22) is

1ME[ηl(χ^l)2]E[ηm(χ^m)2]×{P″[Zi[1]]P[Zi[1]]+P″[Zi[2]]P[Zi[2]]+2P′[Zi[1]]P[Zi[1]]P′[Zi[2]]P[Zi[2]]}.

Summing over loci (and noting that we may include the diagonal terms and only incur an error of order 1/M), we find that, in the case of different parents, the variance of the shared terms Ai+Di, conditional on the trait values of the parent is

12σA2+14σD2−((P′(Zi[1])P(Zi[1])+P′(Zi[2])P(Zi[2]))σA22)2+{P″[Zi[1]]P[Zi[1]]+P″[Zi[2]]P[Zi[2]]+2P′[Zi[1]]P[Zi[1]]P′[Zi[2]]P[Zi[2]]}(σA22)2+O(1M).

Once again we see that if we approximate the distribution of Zi[1] and Zi[2] by that of independent normal random variables with mean z¯0 and variance σA2+σD2, most of these terms cancel and we are left with

σA22+σD24−σA42(σA2+σD2),

exactly as predicted by Theorem C1.

The general case

So far we have only dealt with generation one, where expressions are simplified by the fact that Zi[1], Zi[2] are either identical or independent. More generally, we can perform entirely analogous calculations using Lemma E6 in place of Lemma E2. In the interest of sanity, we omit the details.

Appendix G: Convergence to normal of (Ai+Di) conditional on parental traits

Notation G1 We remind the reader that Notation E1 remains in force. Moreover, since the environmental noise Ei is assumed to be shared by all offspring of the couple i[1], i[2], with this convention we can also assume that the distribution of Ai+Di has a smooth density.

We have verified that the first two moments of the conditional distribution converge to the limits that we would expect if the limit of (Ai+Di) were multivariate normal, but this is not sufficient. To prove that the conditional distribution is indeed asymptotically normal, we appeal to Proposition D3, or rather Corollary D4. We perform the calculation in the case of identical parents, the case of distinct parents being analogous (and less surprising). For definiteness, we consider only generation one. The same argument will work in any generation, but the calculations become considerably more involved, c.f. Lemma E6.

Recall that Ai+Di=∑l=1MΦ(l)/M with Φ defined in Equation (B1). Since we are considering the case of a single parent, Zi[1]=Zi[2]. We shall write Φl(χl1,χl2) when we need to specify the alleles at locus l in Zi[1] on which this is evaluated.

Writing W=∑l=1MΦ(l)/M (as a shorthand for Ai+Di), we write

W^=1M∑l=1MΦ(l)|i[1]=i[2],Zi[1];

that is W^ is the random variable W in the ith individual, conditional on it being produced by selfing and on the parental trait value. This is the quantity that we should like to prove is normally distributed. The first step is to find a suitable exchangeable pair. We write Φ^(l) for the conditioned version of Φ(l).

For each l∈{1,…,M}, let Φ^*(l) be an independent draw from the conditional distribution of Φ^(l) given the sum of Φ^(m) over all m≠l; that is, in an obvious notation, Φ^*(l) has the same distribution as

Φ(l)|∑m≠lΦ^(m),i[1]=i[2],Zi[1].

Now let L be a uniform random variable on {1,…,M} and define

W^′=W^−(Φ^(L)−Φ^*(L))M.

Then (W^,W^′) is an exchangeable pair.

Observe that

(G1) E[W^−W^′|W^]=E[1M1M∑l=1M(Φ^(l)−Φ^*(l))|W^]=1MW^−1M1M∑l=1ME[Φ^*(l)|W^]:=1MW^−T(W^).

Remark G2 We wish to apply Corollary D4. Our first instinct is to write E[W^′|W^]=W^(1−1/M)+T(W^) and take λ=1/M in Equation (D2). This will not suffice, as, with this choice, the first term on the right of Equation (D3) will be too big. As we shall see, the resolution is to take a larger value of λ which captures the dependence of W^′ on W^.

Before we can apply Corollary D4, we need to investigate T. The first step is to establish the distribution of χl1,χl2 conditional on i[1]=i[2], Zi[1] (which we shall for the rest of this section abbreviate to Z) and W−l.

Keeping in mind Notation G1, and recalling that (Z,W) is shorthand for (Zi[1],Ai+Di), we write P[z,w] for the density function of (Z,W) evaluated at (z,w) and Pz[z,w], Pw[z,w], and so on, for the corresponding partial derivatives. The proof of the following lemma mirrors those of Appendix E.

Lemma G3 The (unconditional) distribution of (Z−l,W−l) can be written asP[Z−l=z,W−l=w]=P[z,w]+1ME[Ψl]Pz[z,w]+1ME[Φl]Pw[z,w]+1M(E[Ψl]2−E[Ψl2])Pzz[z,w]+2M(E[Φl]E[Ψl]−E[ΦlΨl])Pzw[z,w]+1M(E[Φl]2−E[Φl2])Pww[z,w]+12ME[Ψl2]Pzz[z,w]+1ME[ΦlΨl]Pzw[z,w]+12ME[Φl2]Pww[z,w]+O(1M3/2).

Proof of Lemma G3 The key, as usual, is Taylor’s Theorem.P[Z−l=z,W−l=w]=∫∫P[χl1=x,χl2=x′,Z=z+1MΨl(x,x′),W=w+1MΦl(x,x′)]dxdx′=∫∫P[χl1=x,χl2=x′,Z=z,W=w]dxdx′+1M∫∫Ψl(x,x′)∂∂zP[χl1=x,χl2=x′,Z=z,W=w]dxdx′+1M∫∫Φl(x,x′)∂∂wP[χl1=x,χl2=x′,Z=z,W=w]dxdx′+12M∫∫Ψl2(x,x′)∂2∂z2P[χl1=x,χl2=x′,Z=z,W=w]dxdx′+1M∫∫Ψl(x,x′)Φl(x,x′)∂2∂z∂wP[χl1=x,χl2=x′,Z=z,W=w]dxdx′+12M∫∫Φl2(x,x′)∂2∂w2P[χl1=x,χl2=x′,Z=z,W=w]dxdx′+O(1M3/2).

Now writeP[χl1=x,χl2=x′,Z=z,W=w]=P[χl1=x,χl2=x′]×P[Z−l=z−1MΨl(x,x′),W−l=w−1MΦl(x,x′)].

Using the notation P[x,x′,z,w]:=P[χl1=x,χl2=x′,Z=z,W=w], we substitute from above and apply Taylor’s Theorem to obtain,P[x,x′,z,w]=P[χl1=x,χl2=x′]×∫∫P[y,y′,z−1MΨl(x,x′)+1MΨl(y,y′),w−1MΦl(x,x′)+1MΦl(y,y′)]dydy′=P[χl1=x,χl2=x′]{P[Z=z,W=w]1M∫∫+1M∫∫(Ψl(y,y′)−Ψl(x,x′))∂∂zP[y,y′,z,w]dydy′+1M∫∫(Φl(y,y′)−Φl(x,x′))∂∂wP[y,y′,z,w]dydy′+O(1M)}.

Differentiating with respect to z (and assuming sufficient regularity),∂∂zP[χl1=x,χl2=x′,Z=z,W=w]=P[χl1=x,χl2=x′]×{∂∂zP[Z=z,W=w]+O(1M)+1M∫∫(Ψl(y,y′)−Ψl(x,x′))∂2∂z2P[y,y′,z,w]dydy′+1M∫∫(Φl(y,y′)−Φl(x,x′))∂2∂z∂wP[y,y′,z,w]dydy′},

and similarly∂∂wP[χl1=x,χl2=x′,Z=z,W=w]=P[χl1=x,χl2=x′]×{∂∂wP[Z=z,W=w]+O(1M)+1M∫∫(Ψl(y,y′)−Ψl(x,x′))∂2∂z∂wP[y,y′,z,w]dydy′+1M∫∫(Φl(y,y′)−Φl(x,x′))∂2∂w2P[y,y′,z,w]dydy′}.

As in the proof of Lemma E2, we only require the second derivatives to leading order∂2∂z2P[χl1=x,χl2=x′,Z=z,W=w]=P[χl1=x,χl2=x′]∂2∂z2P[Z=z,W=w]+O(1M),

with similar expressions for the other second partial derivatives. Substituting back into the first display yields the result. □

Lemma G4 The conditional distribution of χl1,χl2 given Z and W−l is given by(G2) P[χl1=x,χl2=x′|Z=z,W−l=w−l]=P[χl1=x,χl2=x′]{1−Ψl(x,x′)MP[z,w−l]Pz[z,w−l]+E[Ψl]MP[z,w−l]Pz[z,w−l]}+O(1M).

Proof of Lemma G4 This is just an application of Bayes’ rule:(G3) P[χl1=x,χl2=x′|Z=z,W−l=w−l]=P[Z=z,W−l=w−l|χl1=x,χl2=x′]P[Z=z,W−l=w−l]P[χl1=x,χl2=x′]=P[Z−l=z−Ψl(x,x′)M,W−l=w−l]P[Z=z,W−l=w−l]P[χl1=x,χl2=x′].

Using Lemma G3 and Taylor’s Theorem,P[Z−l=z−Ψl(x,x′)M,W−l=w−l]=P[Z=z−Ψl(x,x′)M,W−l=w−l]+1ME[Ψl]Pz[z−Ψl(x,x′)M,w−l]+1ME[Φl]Pw[z−Ψl(x,x′)M,w−l]=P[Z=z,W=w−l]−Ψl(x,x′)MPz[z,w−l]+E[Ψl]MPz[z,w−l]+E[Φl]MPw[z,w−l]+O(1M).

When we integrate this expression with respect to x and x′, to calculate the denominator in Equation (G3), we recover Pz[z,w−l]+E[Φl]MPw[z,w−l]+O(1M) (since the expectation of Ψl cancels). Expanding the ratio in Equation (G3) in powers of 1/M, the terms involving E[Φl] cancel, and the result follows.  □

Finally, we are in a position to calculate the quantity T(W^) that was defined in Equation (G1). Recall that Φ^l* is an independent draw from the conditional distribution of Φ^l given W−l and Z, so using Equation (G2),

(G4) E[Φ^l*|W^]=E[Φl]−(E[ΦlΨl]−E[Φl]E[Ψl])×1ME[1P[Z,W−l]∂∂zP[Z,W−l]|W=W^]+O(1M).

Conditioning only on i[1]=i[2], using the calculations in Appendix B and Equation (F6), by an application of Theorem D2 (up to an error of order 1/M) the joint distribution of (Ai+Di,Zi[1]) is approximately that of a bivariate normal.

We will need that for a bivariate normal distribution with mean vector (μZ,μW) and covariance matrix

(σZ2Cov(Z,W)Cov(Z,W)σW2)

the density function takes the form

p(z,w)=12πσZσW1−ρ2×exp{−12(1−ρ2)((z−μZ)2σZ2−2ρ(z−μZ)(w−μW)σZσW+(w−μW)2σW2)},

where ρ=Cov(Z,W)/(σZσW). Differentiating, we find

(G5) 1p(z,w)∂∂zp(z,w)=1(1−ρ2){ρ(w−μW)σZσW−(z−μZ)σZ2}.

Recall the definition of T from Equation (G1). Multiplying Equation (G4) by 1/M, observing that E[W−l|W,Z]=W+O(1/M) (and since Φl is uniformly bounded independent of l, the error is bounded independent of l), and then averaging out over l as in the definition of T(W^), on substituting Equation (G5) and Cov(Z,W)=ρσZσW, we find

(G6) T(W^)=1ME[W]+Cov(Z,W)M{Z−E[Z]σZ2(1−ρ2)−ρ1−ρ2W^−E[W]σZσW}+O(1M3/2)=1ME[W]+ρσZσWM{Z−E[Z]σZ2(1−ρ2)−ρ1−ρ2W^−E[W]σZσW}+O(1M3/2)=1ME[W]+1MρσWσZZ−E[Z]1−ρ2−1Mρ21−ρ2(W^−E[W])+O(1M3/2).

Using the approximation for the conditional distribution of (χl1,χl2), given Z obtained in Appendix E,

E[W^]=E[W|Z]=E[W]+ρσwσz(Z−E[Z])+O(1M),

so we can rewrite Equation (G6) as

T(W^)=1M1(1−ρ2)E[W^]−1Mρ21−ρ2W^+O(1M3/2).

Substituting in Equation (G1),

(G7) E[W^−W^′|W^]=1M11−ρ2(W^−E[W^])+O(1M3/2).

We are going to apply Corollary D4 to (W^,W^′) with F=σ(W^). We set λ=1/(M(1−ρ2)) and observe from Equation (G7) that we may take a remainder term R with E[|R|] of order 1/M1/2 in Equation (D3). Moreover,

K^2=12λ|Δ|32,

and so, since by construction |Δ|<C/M, E[K^2] is also order at most 1/M1/2.

Since with these definitions

K^1=M(1−ρ2)2E[(W^−W^′)2|W^],

it remains to control

(G8) E[|σW^2−M(1−ρ2)2E[(W^−W^′)2|W^]|].

Again using the results of Appendix E,

σW^2=(1−ρ2)σW2+O(1M),

(the first term being the conditional variance if the random variables were distributed exactly as a bivariate normal), whereas

E[E[(W^−W^′)2|W^]]=E[(W^−W^′)2]=E[1M2∑l=1M(Φ^l−Φ^l*)2]=21MσW2+O(1M3/2).

(Note that we see the unconditioned σW2 in this second expression since it involves only diagonal terms.)

To control Equation (G8), observing that, by the Cauchy–Schwarz inequality,

E[|E[M(W^−W^′)2]−E[M(W^−W^′)2|W^]|]≤Var(E[M(W^−W^]′)2|W^])1/2,

it suffices to control

Var(E[M(W^−W^′)2|W^]).

In particular, we should like to show that this expression is of order O(1/M).

Now we use the standard decomposition of conditional expectations: for two random variables X and F,

Var(X)=E[E[X2|F]−(E[X|F])2+(E[X|F])2]−E[E[X|F]]2=E[Var(X|F)]+Var(E[X|F]).

So

Var(E[X|F])=Var(X)−E[Var(X|F)].

For us, X=M(W^−W^′)2=(ΦL−ΦL*)2, and F=W^, so

Var(X)=1M∑l=1ME[(Φl−Φl*)4]−(1M∑l=1ME[(Φl−Φl*)2])2,

and we seek

1M∑l=1ME[(Φl−Φl*)4]−(1M∑l=1ME[(Φl−Φl*)2])2−1M∑l=1ME[E[(Φl−Φl*)4]|W^]+E[(1M∑l=1ME[(Φl−Φl*)2|W^])2],

where the expectation is with respect to the distribution of W^. By the tower property, the terms involving (Φl−Φl*)4 cancel, leaving

(G9) 1M2∑l=1M∑m=1M{E[E[(Φl−Φl*)2|W^]E[(Φm−Φm*)2|W^]]−E[(Φl−Φl*)2]E[(Φm−Φm*)2]}.

Expanding E[(Φl−Φl*)2|W^]E[(Φm−Φm*)2|W^] in an entirely analogous way to Equation (G4), when we take expectations, using the tower property of conditional expectations, the part of the product that is an affine function of W^ will cancel in Equation (G9), leaving quadratic (and higher order) terms, each of which is of order O(1/M) in the summand. Overall then Equation (G9) is O(1/M), and applying Corollary D4, the proof that Ai+Di is normal with an error of order 1/M is complete.

The residuals, generation one

Proving that RAi+RDi is normal is much simpler. Since Mendelian inheritance is independent across loci, we are able to use Theorem D2 in much the same way as in generation zero. A combination of Lemma E2 and Bayes’ rule suffices to show that the variance is not affected by conditioning on parental trait values, after which the proof proceeds essentially as in the additive case and so is omitted.

Appendix H: Generation t, accumulation of information

If we wanted to prove a strict analog of the results of Barton et al. (2017) in the additive case, then we would want to condition not just on the trait values of the parents, but on the trait values of an arbitrary collection of individuals in the pedigree. Such a proof can follow essentially the same lines as above, although the calculations are considerably longer to write out. The only thing that must be checked is that we do not accumulate too much information from knowing those trait values; it is this that controls for how long the infinitesimal approximation will remain accurate. This requires more care than the additive case of Barton et al. (2017), so we present the argument here.

Recall that we write P(t) for the pedigree up to and including generation t and Z(t) for the corresponding vector of trait values of all individuals in P(t). We would like to understand the distribution of the allelic types χl1(j*),χl2(j*) at locus l of an individual j* in generation t, conditional on knowing the trait values of all individuals in the pedigree up to generation t−1. That is, we would like to estimate

(H1) P[(χl1(j*),χl2(j*))=(x,x′)|P(t),Z(t−1)=(zj)j∈P(t−1)]=P[Z(t−1)=(zj)j∈P(t−1)|(χl1(j*),χl2(j*))=(x,x′),P(t)]P[Z(t−1)=(zj)j∈P(t−1)|P(t)]×P[(χl1(j*),χl2(j*))=(x,x′)|P(t)].

To estimate the numerator in the fraction, we partition over the possible patterns of identity at locus l in the pedigree, conditional on that pedigree; that is we condition on the values of the Bernoulli random variables that determine Mendelian inheritance at locus l across the pedigree. We denote this Ml(t) and abuse notation by writing (Ml1(j),Ml2(j)) for the allelic states at locus l in individual j∈P(t−1) conditional on Ml(t). More precisely, if χl1(j*)=x and χl2(j*)=x′, (Ml1(j),Ml2(j))=(y,y′),(y,x′),(x,y′),(x,x′) according to whether j is identical by descent with the chosen individual j* on neither chromosome, one chromosome or both chromosomes. We use EMl when we wish to emphasize that we are taking the expectation with respect to this quantity. We proceed as in Lemma E6:

P[Z(t−1)=(zj)j∈P(t−1)|(χl1(j*),χl2(j*))=(x,x′),P(t),Ml(t)]=P[Z−lj=zj−1MΨl(Ml1(j),Ml2(j)),∀j∈P(t−1)|P(t)]=E[P[Zj=zj−1MΨl(Ml1(j),Ml2(j))+1MΨl(χl1(j),χl2(j)),∀j∈P(t−1)|P(t)]],

where in the last line the expectation is taken with respect to the unconditional law of the random family {(χl1(j),χl2(j)),j∈P(t−1)}.

Substituting in Equation (H1), in an obvious notation,

P[(χl1(j*),χl2(j*))=(x,x′)|P(t),Z(t−1)=z]=P[(χl1(j*),χl2(j*))=(x,x′)|P(t)]×(1−∑j∈P(t−1)1M{E[Ψl(Ml1(j),Ml2(j))|P(t)]−E[Ψl(χl1(j),χl2(j))|P(t)]}PZj[z]P[z]∑j∈P(t−1))+O(1M).

In particular, the summand will vanish if j and j* are not identical by descent in at least one copy at locus l, since then the allelic states at locus l in individuals j and j* are independent. Furthermore, the more distant the relationship between j and j* (that is, the smaller the probability of their being identical by descent), the less information we glean about the allelic states in j* from observing the trait value of individual j, resulting in a small contribution of the jth term to the difference between the conditional and unconditional laws of (χl1(j*),χl2(j*)). The infinitesimal model can be expected to break down for an individual if we know that one of its close relatives had a particularly extreme trait value, or if the pedigree is particularly inbred (so that there is little variation between offspring).

Appendix I: Supplementary material and codes

The following supplementary material can be found in the public repository (Barton 2023):

The Mathematica notebook Algorithm for calculating identities.nb, comprising a set of codes to compute the identity coefficients of Appendix A.

The Mathematica notebook Infinitesimal with dominance.nb, accompanying and complementing the simulations and figures presented in the paper.

The different datasets from which the numerical examples in this paper can be reproduced.
==== Refs
Literature cited

Abney  M . A graphical algorithm for fast computation of identity coefficients and generalized kinship coefficients. Bioinformatics. 2009;25 :1561–1563. doi:10.1093/bioinformatics/btp185 19359355
Abney  M, McPeek  MS, Ober  C. Estimation of variance components of quantitative traits in inbred populations. Am J Hum Genet. 2000;66 :629–650. doi:10.1086/302759 10677322
Barton  NH . The Infinitesimal Model with Dominance—Codes and Data. Vienna: ISTA; 2023.
Barton  NH, Etheridge  AM. The relation between reproductive value and genetic contribution. Genetics. 2011;188 :953–973. doi:10.1534/genetics.111.127555 21624999
Barton  NH, Etheridge  AM. Establishment in a new habitat by polygenic adaptation. Theor Pop Biol. 2018;122 :110–127. doi:10.1016/j.tpb.2017.11.007 29246460
Barton  NH, Etheridge  AM, Véber  A. The infinitesimal model: definition, derivation, and implications. Theor Pop Biol. 2017;118 :50–73. doi:10.10.16/j.tpb.2017.06.00128709925
Boyle  EA, Li  YI, Pritchard  JK. An expanded view of complex traits: from polygenic to omnigenic. Cell. 2017;169 :1177–1186. doi:10.1016/j.cell.2017.05.038 28622505
Brockwell  PJ, Davis  RA. Introduction to Time Series, Forecasting. Springer Texts in Statistics. New York: Springer-Verlag; 1996.
Bulmer  MG . The effect of selection on genetic variability. Am Nat. 1971;105 :201–211.
Charlesworth  B . Causes of natural variation in fitness: evidence from studies of Drosophila populations. Proc Natl Acad Sci USA. 2015;112 (6 ):1662–1669. doi:10.1073/pnas.1423275112 25572964
Chen  L, Goldstein  L, Shao  Q-M. Normal Approximation by Stein’s Method. Berlin, Heidelberg: Springer; 2011.
Fisher  RA . The correlation between relatives on the supposition of Mendelian inheritance. Proc R Soc Edinb. 1918;52 :399–433.
García-Cortés  LA . A novel recursive algorithm for the calculation of the detailed identity coefficients. Genet Sel Evol. 2015;47 :33. doi:10.1186/s12711-015-0108-6 25926309
Hill  WG, Barton  NH, Turelli  M. Prediction of effects of genetic drift on variance components under a general model of epistasis. Theor Pop Biol. 2006;70 :56–62. doi:10.1016/j.tpb.2005.10.001 16360188
Karigl  G . A recursive algorithm for the calculation of identity coefficients. Ann Human Genet. 1981;45 :299–305. doi:10.1111/j.1469-1809.1981.tb00341.x 7305283
Karigl  G . A mathematical approach to multiple genetic relationships. Theor Popul Biol. 1982;21 :379–393. doi:10.1016/0040-5809(82)90025-9
Kirkpatrick  B, Ge  S, Wang  L. Efficient computation of the kinship coefficients. Bioinformatics. 2019;35 :1002–1008. doi:10.1093/bioinformatics/bty725 30165566
Lande  R . The maintenance of genetic variability by mutation in a polygenic character with linked loci. Genet Res. 1975;26 :221–235. doi:10.1017/S0016672300016037 1225762
Lande  R, Porcher  E. Maintenance of quantitative genetic variance under partial self-fertilization, with implications for evolution of selfing. Genetics. 2015;200 (3 ):891–906. doi:10.1534/genetics.115.176693 25969460
Lange  K . Central limit theorems of pedigrees. J Math Biol. 1978;6 :59–66. doi:10.1007/BF02478517
Rinott  Y . On normal approximation rates for certain sums of dependent random variables. J Comput Appl Math. 1994;55 :135–143. doi:10.1016/0377-0427(94)90016-7
Robertson  A . A theory of limits in artificial selection. Proc R Soc Lond B. 1960;153 :234–249. doi:10.1098/rspb.1960.0099
Roze  D . Background selection in partially selfing populations. Genetics. 2016;203 :937–957. doi:10.1534/genetics.116.187955 27075726
Sachdeva  H . Effect of partial selfing and polygenic selection on establishment in a new habitat. Evolution. 2019;73 :1729–1745. doi:10.1111/evo.13812 31339550
Sachdeva  H, Barton  NH. Introgression of a block of genome under infinitesimal selection. Genetics. 2018;209 :1279–1303. doi:10.1534/genetics.118.301018 29895560
Santiago  E . Linkage and the maintenance of variation for quantitative traits by mutation-selection balance: an infinitesimal model. Genet Res. 1998;71 (2 ):161–170. doi:10.1017/S0016672398003231
Stein  C . Approximate Computation of Expectations. Lecture Notes—Monograph Series. Institute of Mathematical Statistics; 1986.
Turelli  M, Barton  NH. Genetic and statistical analyses of strong selection on polygenic traits: what, me normal?  Genetics. 1994;138 (3 ):913–941. doi:10.1093/genetics/138.3.913 7851785
Walsh  JB, Lynch  M. Evolution and Selection of Quantitative Traits. Oxford: Oxford University Press; 2018.
