Model:

We consider a population in which interactions occur within “groups” of n individuals each and affect phenotypic values of individuals. In classical quantitative genetic theory, the phenotypic trait value of an individual, P, is the sum of a heritable component, A, referred to as “breeding value” and a residual component, E, referred to as “environment”: P = A + E. With interactions among n individuals, however, the observed phenotypic trait value of an individual is due to two unobserved phenotypic effects: a direct effect originating from the focal individual and the sum of the associative effects originating from of each of its n − 1 group members (Griffing 1967). Associative effects may be interpreted as heritable environmental effects provided by associates to the focal individual (Wolf 2003). Both phenotypic direct and associative effects are the sum of a heritable component (A) and a nonheritable component (E), so that the observed phenotypic trait value of an individual is given by(1)(Griffing 1967), in which i denotes the focal individual and j one of its associates (i.e., group members). The first two terms in Equation 1 represent effects due to the focal individual; the last two terms represent the effects due to its n − 1 associates. The effect is the direct breeding value (DBV) and E_D,i is the direct nonheritable value of individual i, whereas is the associative breeding value (SBV) and E_S,j the associative nonheritable value of associate j. The DBV and the SBV represent the heritable components underlying the observed phenotypes.

Note that nonheritable associative effects, indicated by in Equation 1, may contribute to the observed phenotype. This is because the full associative effect of an individual is a phenotypic effect, which may contain both heritable and nonheritable components. The genes of an individual do not directly affect the associates; the associates experience a phenotype. We emphasize those nonheritable associative effects, because they modify the statistical model required to obtain unbiased estimates of genetic parameters (see below).

Heritable variation:

Because each individual interacts with n − 1 associates, the total heritable impact of an individual on the population mean, referred to as its total breeding value (TBV), equals the sum of its DBV and n − 1 times its SBV: TBV_i = DBV_i + (n − 1)SBV_i. Response to selection equals the change per generation of the TBV. Hence, the TBV is a generalization of the classical breeding value, to account for interactions. The total heritable variation available for response to selection, therefore, equals the variance of the TBV,(2)(Bijma et al. 2007), in which is the variance of DBV, is the variance of SBV, and is the covariance between DBV and SBV. The sign of the covariance between DBV and SBV is a measure of competition vs. cooperation. Negative may be interpreted as “heritable competition,” in the sense that individuals with positive breeding values for their own phenotype (DBV) have on average a negative heritable impact on the phenotypes of their associates (SBV). Conversely, positive may be interpreted as “heritable cooperation.”

Classical quantitative genetic methodology to estimate additive genetic variances of traits rests on similarity of phenotypes' relatives. Associative breeding values of individuals, however, are not expressed in their own phenotypes, but in the phenotypes of their associates. Thus only the variance of DBV surfaces in classical analyses; the remaining part of the full heritable variance, , is hidden. Depending on group size and the heritability of interactions, this hidden heritable variance may be substantially larger than the variance of DBVs. In those cases, classical methods will substantially underestimate the total heritable variation. With specific family structures, the SBV may also be confounded with the DBV (see discussion and Wolf 2003), which would make the result of classical analyses an inaccurate estimate of even the direct genetic variance. The total heritable variation (Equation 2) does not depend on relatedness among associates. This is because all components of the TBV of an individual refer to its own genes. (See Bijma et al. 2007 for the effect of population structure.)

Response to selection:

In classical quantitative genetic theory, response to selection is in the same direction as selection pressure; heritability is merely a scaling factor affecting the magnitude of response, not its direction. With interactions among individuals, however, not only the magnitude, but also the direction of response depend on the genetic parameters (Griffing 1967; Kirkpatrick and Lande 1989; Moore et al. 1997). To illustrate this phenomenon, we consider response to selection among unrelated individuals, which follows from substituting r = g = 0 into Equation 5 of Bijma et al. (2007), giving , in which is the selection gradient. (See also Griffing 1976a.) This result shows that the direction of response depends on the covariance between DBV and SBV (). Strong competition, i.e., large and negative , will yield response in the direction opposite to the direction of selection. With interactions, therefore, the genetic parameters play a key role in the outcome of selection; they are no longer merely a scaling factor.

Derivations:

This section describes the derivations of Equations 4 and 5. Consider two individuals, i with group members denoted by j and k with group members denoted by l. The covariance between phenotypes of i and k follows from substituting Equation 1 for both P_i and P_k, giving , which consists of a genetic and an environmental component. Analogous to the basic animal model, we account for the genetic components by including direct and associative breeding values in the model, together with their covariance structure. This yields Equation 4. For example, in Equation 4, the ith element of y corresponds to P_i, the ith element of Z_Da_D corresponds to A_D,i, the ith element of Z_Sa_S corresponds to , and the ith element of e corresponds to . As usual in an animal model, we account for covariances between breeding values in a_D and a_S by incorporating the relatedness matrix A in the statistical analyses. (See below and the mixed-model equations in Muir 2005.)

What remains to be done is the covariance structure of the vector of residuals e. The covariance between residuals in e originates from the environmental component of the covariance between observations and can be partitioned into . As usual in quantitative genetics, we assume that E-terms of distinct individuals are independent. (Environmental effects common to multiple individuals are accounted for by the fixed effects Xb.) Thus a nonzero residual covariance occurs only when it involves the same individual, which happens either when i and k are the same individual or when i and k are group members. When i and k are the same individual, , which is the residual variance. The first term is due to the direct effect of the individual itself, and the second term is due to its n − 1 group members. When i and k are group members, their residuals are correlated for two reasons. First, for groups larger than two, i and k will have n − 2 group members in common, and thus they have n − 2 identical associative effects in their phenotype, giving . Second, the direct effect of i is expressed in its own phenotype, whereas the associative effect of i is expressed in the phenotype of k, and vice versa for k. As with any two traits of a single individual, the direct and associative effect originating from a single individual may be correlated Walsh (Lynch and Walsh 1998), giving . Thus the full covariance between residuals of group members equals . The represents the nonheritable relationship between the direct and associative effect of an individual and is a nonheritable measure of competition vs. cooperation (analogous to , see above). A negative may, for example, arise from variation in the juvenile environment experienced by individuals. Individuals having experienced good (bad) environments may be more (less) competitive, in which case they will have higher (lower) phenotypic values at the expense (to the advantage) of their associates. In Equation 5, we accounted for the nongenetic covariance between group members by fitting the correlation ρ between records of group members.

Experimental design for application:

The use of Equations 4 and 5 does not require a specific experimental design; groups may be composed of any combination of related and/or unrelated individuals and may vary in size. To calculate the relationship matrix A, information on at least one generation of family relationships will be required, so that one can distinguish full sibs, half sibs, and unrelated individuals. Information on more generations of pedigree, which provides information on more distant relationships, will add to the accuracy of results, but is not required for application. With groups of equal size, there will be no information in the data to separate and , because different combinations of those two parameters may yield an identical residual correlation. The inability to estimate and with equal group size is no problem, because the residual correlation itself, which included these components as part of a ratio, is estimable and sufficient in the analysis to avoid bias in the estimates of the genetic parameters. There are population structures from which the parameters of interest cannot be estimated. (See discussion on the work of Wolf 2003).

We used Monte Carlo simulations to validate our statistical methodology (appendix). Table 1 shows that estimated genetic parameters obtained from simulated data are unbiased; true values of genetic parameters are within ±3 standard errors of mean estimated values in all cases. These results confirm that Equations 4 and 5 correspond to the covariance structure that arises from Equation 1 and that there is no need for a specific structure of the family relationships within and between groups.

Material:

The methodology presented above was applied to analyze mortality in a commercial population of layer chickens. Field data on this population were provided by Hendrix Genetics and consisted of observations on survival days of a single generation of 3800 hens bred from 36 sires and 287 dams, which had been mated at random. Each sire had been mated to ∼8 dams, and each dam contributed on average 13.2 female offspring. Five generations of pedigree were included in the calculation of the matrix of relatedness (A).

Mortality in commercial chicken populations is largely due to cannibalistic pecking behavior (e.g., Craig and Muir 1996). The direct effect of an individual is, therefore, related to its ability to survive by avoiding being pecked, whereas its associative effect measures its effect on survival of its cage members, which is closely related to its own pecking behavior. In commercial egg production, individuals are frequently beak trimmed, to reduce mortality due to pecking behavior. In this data set, however, individuals had not been beak trimmed, which will be obligatory in the European Union from 2011 onward.

Data consisted of four age groups of approximately equal numbers of hens, differing by 2 weeks of age each. At an age of 20 weeks, individuals were allocated randomly to 950 standard commercial battery cages, four individuals per cage. Due to chance, some cages contained full or half sibs, but most cages contained unrelated individuals only. Facilities consisted of eight rows of cages. Each row was 40 cages long and consisted of 3 cages on top of each other, giving 120 cages per row. The eight rows were grouped into four double rows; both rows of a double row were put with the backs against each other. Each age group was allocated to a single double row. Back fences of cages consisted of latticework allowing limited contact between “backdoor neighbors,” i.e., between the individuals of two corresponding cages, each in one of the two rows of a double row. Adjacent cages in the same row were separated by a closed fence fully preventing contact. Individuals were kept 60 weeks and culled at an age of 80 weeks. Survival (0/1) was recorded daily on all individuals. For each individual, survival days were defined as age at death in days. Individuals surviving until 80 weeks of age had 80 × 7 = 560 survival days. Fifty-four percent of individuals survived until 80 weeks of age. Mean survival time was 454 days, with a standard deviation of 122 days. Inspection of dead hens showed that the vast majority of chickens had died due to being pecked. This observation was supported by previous data on beak-trimmed individuals of the same line kept in the same facilities, which showed substantially lower mortality (<10%). Thus the absence of beak trimming substantially amplified the consequences of social interactions.

Data were analyzed with Equations 4 and 5 as a model, using restricted maximum likelihood (REML) (Patterson and Thompson 1971) as implemented in the ASReml software package (Gilmour et al. 2002). The model included a fixed effect for each combination of row (1 through 8) and level (top, middle, and bottom) of the cage and for the number of surviving birds in the back cage, to account for a possible associative effect due to the backdoor neighbors. As usual with an animal model, this analysis directly yields estimates of the parameters of interest and utilizes the full data.

Results:

When using a conventional model without associative effects, estimated heritability for survival days in the chicken line was 6.7%, which is in the range common for survival traits in livestock (Table 2; Beaumont et al. 1997; Carlen et al. 2005; Gitterle et al. 2005). The total heritable variance estimated from the full model (Equations 2, 4, and 5), however, was 2742 days², which corresponds to 20% of the phenotypic variance. This result indicates a “heritability” threefold greater than that obtained using classical methods, i.e., 20% vs. 6.7%. Two-thirds of the heritable variation, therefore, is due to social interactions among individuals and is hidden from classical analyses.

At first glance, one might expect direct and associative effects for mortality due to pecking behavior to be negatively correlated, because dominant individuals may kill others and survive themselves. However, we found a small but nonsignificant positive correlation of +0.28 (P = 0.016), suggesting that individuals benefit from not harming others. In this population, there seems to be mutual benefit of behaving peacefully and a mutual cost of behaving aggressively. A positive covariance between DBV and SBV increases the total heritable variation (Equation 2) and therefore increases the potential of the population to respond to selection.

Figure 1 shows predicted responses to multilevel selection applied to the chicken population, using equations presented in Bijma et al. (2007). (Details and an example of calculation are in the appendix.) In Figure 1, the strength of multilevel selection is indicated by a factor g, where g = 0 represents selection among individuals, g = 1 represents selection among groups, and values of g between 0 and 1 represent a mix of selection among individuals and groups (Bijma et al. 2007). Classical theory suggests a response of only 7.8 days of survival. Predictions accounting for heritable interactions, however, yield substantially higher responses. Selection among individuals applied to a population composed of groups of unrelated individuals yields an expected response of ∼11 days of survival (g = 0, r = 0). Mild multilevel selection applied to a population composed of groups of full sibs more than doubles the predicted response (g ≥ 0.3, ). The maximum response that can be obtained equals 23 days (g = 0.6–0.8 and ), which is nearly threefold greater than suggested by classical theory. (Clones do not occur in reality.) The lines in Figure 1 are curved, which indicates that a degree of between-group selection (g) that maximizes response exists, which was also observed by Griffing. This value of g may be applied in artificial breeding programs, where the aim is to maximize response to selection. Note that results in Figure 1 are specific to the chicken population and to the environment in which it was kept, such as the density of housing; other populations may have other genetic parameters and other relationships among relatedness, multilevel selection, and response. Results in Figure 1 should, therefore, not be taken to conclude that relatedness has a bigger impact than multilevel selection in general.

• Download figure
• Open in new tab
• Download powerpoint

Figure 1.—

Predicted responses to selection for the trait survival days in the chicken population, with groups consisting of one of four types of relatives. Predictions are obtained using results in Table 2 and a selection intensity of 1 unit . (Methods are in the appendix). The g denotes the strength of multilevel selection. The prediction from classical theory is . Results for clones, which do not occur in chicken, are included to enable generalization to other species.

Monte Carlo simulation used for Table 1:

A population of five discrete generations was simulated. First, a base generation of 250 sires and 500 dams was generated. DBVs of base individuals were sampled from and nongenetic direct effects from . Associative effects of individuals were simulated conditional on their direct effects; the SBV of the ith individual of the base generation was sampled as , in which r_A is the correlation between DBV and SBV, and z is a random variable taken from a standard normal distribution, . Nongenetic indirect effects were sampled as , in which r_E is the nongenetic correlation between direct and associative effects. Phenotypes were created according to Equation 1 and parents for the next generation were taken at random from the base generation.

Second, four additional generations were generated. DBVs of individuals in those generations were sampled as and SBVs as , in which MS denotes Mendelian sampling components of breeding values, which were sampled as and . Nongenetic effects were sampled in the same way as in the base generation. Phenotypic values from the fifth generation were used to estimate variance components, and all five generations were included in the calculation of the relatedness matrix (A). Fifty replicates were simulated for each set of genetic parameters and estimates were averaged over replicates.

Example of calculation of response to selection for the chicken line:

Estimated genetic parameters for the chicken line were , , , , and ρ = 0.09. In the data, groups consisted of four unrelated individuals, so that the phenotypic variance equals (Table 2).

In the following, we give an example of calculation of response to selection for a population consisting of groups composed of four full sibs (n = 4, r = ), with substantial selection pressure given to group performance (g = 0.8), and for an intensity of selection of 1 unit (ι = 1). The general expression for response to selection is , in which C denotes the selection criterion, , and g the relative emphasis given to group performance in the selection criterion (Bijma et al. 2007). Thus, calculation of response to selection requires calculating the variance of the total breeding values, , the covariance between phenotypic values and TBVs of individuals, , and the selection gradient, . We start with the selection gradient and subsequently calculate and .

The denominator of the selection gradient, the standard deviation of the selection criterion, , depends on the phenotypic variance, , and on the covariance between phenotypes of group members, . The general expression for the phenotypic variance follows from Equation 1, giving . The covariance between phenotypes of associates also follows from Equation 1, . Next, it follows from that . Thus the selection gradient equals .

The variance of the total breeding value equals (Equation 2). The covariance between phenotypic values and TBVs of individuals equals (Bijma et al. 2007). Finally, response to selection follows from .

Expected sire and dam covariances in the design of Wolf (2003):

In the experimental setup of Wolf (2003), flies were kept in pairs of two individuals that were full sibs, half sibs, or unrelated. With unrelated individuals, the same two half-sib families were always combined in a pair. Dams were nested within sires; each sire was mated to exactly two dams and each dam to a single sire only.

Investigation of the expectations of sire and dam variances presented in Wolf (2003) shows that they are incomplete. Consider, for example, the case where tube members are half sibs. In this case, the sire covariance between two individuals in different groups (tubes) equals , instead of as listed in Table 1 of Wolf (2003), i.e., a difference of . This extra originates from the full-sib relationship that exists between a focal individual and the tube member of a half sib of the focal individual, which is not accounted for in Wolf (2003). When each sire is mated to exactly two dams, the group member of a half sib of the focal individual must be a full sib of the focal individual (in the case where tube members are half sibs). Note that is not a variance, but a covariance. For example, if estimates presented in Wolf (2003) would have been the true values, so that , , and , then the expected sire variance would equal −2666, i.e., a negative value. Hence, a sire plus dam-within-sire model is inappropriate, because the expected sire variance may be negative for combinations of direct and indirect genetic variances that are biologically possible.

To investigate whether direct and associative genetic variances can be estimated from the data structure used in Wolf (2003), we derive general expressions for expectations of sire and dam covariances. We distinguish between (i) covariances between individuals in the same group and (ii) covariances between individuals in different groups. The expected covariance between two individuals, i and j, in the same group is(A1)and in different groups is(A2)in which is relatedness between i and j, is relatedness between i and the group member of j, is relatedness between j and the group member of i, and is relatedness between the group members of i and j. Table A1 shows expected covariances between full sibs, half sibs, and unrelated individuals for each of the three tube compositions used in Wolf (2003). The expectations of sire and dam covariances given in Table 3 differ from those presented in Wolf (2003), primarily because covariances due to group members of full and half sibs are incomplete in Wolf (2003).

It follows from Equations A1 and A2 that and can be separated only when . Thus the information in the data to separate and depends on whether and how frequently occurs. It follows from Table A1 that occurs only in the situation where groups consist of two unrelated individuals and only when relatedness between the two group members of both individuals differs from that between both individuals considered. (see Table A1 footnote). Whether or not that situation occurs in the data depends on details of the experimental implementation that are not described in Wolf (2003). Thus, the data structure used in Wolf (2003) provides either very little or no information at all to separate the direct and the associative heritable variance.

To validate this theoretical result we performed Monte Carlo simulation of the data structure used in Wolf (2003), where, in the case of groups consisting of unrelated individuals, group members of full sibs were full sibs as well, and group members of half sibs were half sibs as well. Hence, theoretically this data structure should not contain any information to separate and . Statistical analyses of the simulated data using Equations 4 and 5 as a model showed that the matrix of second derivatives of the likelihood (Fisher's information matrix) was singular, indicating that there are an infinite number of combinations of parameters that yield the identical maximum likelihood. Indeed identical maximum likelihoods were observed for multiple combinations of and , which confirms the theoretical result derived above.