Estimation of additive genetic and environmental sources of quantitative trait variation using data on married couples and their siblings

Summary

Twin studies have been used to understand the sources of genetic and environmental variation in body height, body weight and other common human quantitative traits. However, it is rather unclear whether these two sources of variation could be really separated in practice. Here, we consider a special study design where phenotype data from married couples and their siblings have been collected. The marital status gives information about the shared environment, while siblings give information about both genetic and environmental variation. To dissect sources of variation and to allow some deviations and pedigree errors in the data, we model such data using a robust polygenic model with finite genome length assumption. As a summary, we provide the estimates for age-dependent proportions of total variation which are due to polygenic and environmental effects. Here, these estimates are provided for body height, weight, systolic blood pressure and total serum cholesterol measured from subjects of the Indian Migration Study.

1. Introduction

Both genetic and common environmental effects are important in understanding sources of variation in quantitative phenotypes (e.g., risk or susceptibility factors for diseases). However, due to lack of evidence for the effects of shared environment, family studies or twin studies may result in over-estimation of genetic effects (Hopper, 2000). Because of this, an extended twin design is often preferred (e.g., Stoel et al., 2006) including phenotypes of the spouses (Eaves et al., 1999). To account for the shared environmental effects, a large number of different types of individuals, their siblings and spouses need to be studied and information on the number of years spent together in a specific/shared environment, age up to which the environment (household) was shared and phenotypes need to be collected. To summarize, the estimated division between the genetics and environment in a population-specific manner, heritability is usually used as a summary statistics to describe the proportion of phenotypic variation attributable to genetic factors.

For quantitative traits, heritability estimation is often performed using a polygenic model (Henderson, 1984; Lynch and Walsh, 1998; Abney et al., 2000). A polygenic model (or infinite locus or infinitesimal model) has been defined (Fisher, 1918) and used by assuming infinite number of additive unlinked loci (and infinite genome length), which implies 0·5 as a constant degree of genetic relationship between full siblings (Fisher, 1918). However, as is well known, the length of the genome is finite and specific to different species. The size of the human genome corresponds to an equivalent of 85 independent loci (Visscher et al., 2006). Because of this, Visscher et al. (2006) first estimated the actual values of relationship coefficients (entries of the covariance matrix) using a set of neutral markers and then incorporated these values into the relationship matrix instead of a fixed 0·5 for full siblings. This practice resulted in much higher heritability estimates, than otherwise expected, based on collected phenotypes for human height, which is a typical quantitative trait (see Visscher et al., 2006; Xu, 2006). The obtained heritability of h ²≈0·8 was close to estimates from twin studies (Silventoinen et al., 2003).

For computational easiness, a finite polygenic model has been proposed as a finite locus approximation for polygenic model (Thompson and Skolnick, 1977; Du et al., 1999; Du and Hoeschele, 2000). However, the drawback in such models is that the heritability estimates are clearly dependent on the number of loci assumed in the model. Here, we propose another kind of finite approximation for a polygenic model, where instead of assuming infinite genome length, elements of the relationship matrix corresponding to full siblings (degree of genetic relationship between siblings) are treated as random variables and estimated simultaneously with the variance components. Unlike Visscher et al. (2006), we do not utilize molecular marker information in estimation, but we assume the mean 0·5 and known standard deviation 0·0384 (cf. Guo, 1996; Visscher et al. 2006) for the unknown coefficients of genetic relationship between full siblings. This formulation provides a more flexible model, which should be robust for deviations and pedigree errors in the data.

In the Indian Migration Study (IMS) (Lyngdoh et al., 2006), factory workers from four different parts of India are examined for various anthropometric and biochemical measurements apart from other characteristics. The spouses of the factory workers and siblings or relatives (in a very few cases when it was not possible to involve siblings) of both the factory workers and their spouses are examined. The main aim of IMS is to study the migration pattern in India; data on the place of origin, age when the place of origin was abandoned, age at the time of marriage and the number of years spent with the spouse in a common household are collected. It is to be noted that the siblings or relatives are taken from the place of origin of the respective factory worker or the spouse. Hence, IMS provided a unique opportunity to study the genetic as well as the shared environmental effects by modelling the dependency between the factory worker and his sibling/relative, factory worker and his spouse and the spouse and her sibling/relative (see Fig. 1).

Fig. 1. IMS study subjects and their dependency due to genetic and environmental sharing.

The present work is motivated by IMS. In the setting of a simple linear mixed model, the effects of shared environment is modelled as a function of the number of years the environment was shared. In the present context, an environment refers to a common household in either a rural or an urban region.

In section 2, we construct a model describing the genetic and environmental effects on the phenotypes. We also define age-dependent proportions of total variations that are due to polygenic and environmental effects as meaningful summary statistics in the present context. In section 3, a Bayesian approach to the estimation of parameters is discussed and, in section 4, the proposed approach is applied to the IMS data. Finally, we conclude with the discussion section.

2. Model

In this section, we construct a model under the setting of the IMS, where the observation unit of the analysis is a group of four individuals defined by a male factory worker, his spouse, his sibling and his spouse's sibling. However, the idea presented here is very general and can be extended easily to other data types by appropriately specifying the relationships between the individuals within the unit of the analysis.

In the sequel, subscripts (1, 2, 3, 4) are used for the factory worker, his spouse, his sibling and his spouse's relative, respectively. Let a 4×1 vector of quantitative phenotypes corresponding to the factory worker i, i=1, 2,…, n and related observations be y _i=(y _i1, y _i2, y _i3, y _i4)^t, where t denotes the transpose. That is, y _i1 corresponds to the measurement of the factory worker i, y _i2 corresponds to the measurement of the spouse of the factory worker i and so forth. Let x ₂=(x _i1, x _i2, x _i3, x _i4)^t be the dummy variables or continuous measurements of the covariates of interest.

Consider a simple linear mixed model

(1)

where μ is the overall mean and β is the effect of covariates x _i. Further, the additive genetic value or the polygenic effect vector, G _i=(G _i1, G _i2, G _i3, G _i4)^t, is assumed to be jointly (multivariate) normally distributed with the mean vector of zeros and variance–covariance matrix

(2)

Here, A _i is the additive genetic relationship matrix (Henderson, 1984), which characterizes the degree to which related individuals jointly share their genomes identical-by-descent (IBD) and is given by

and sibling covariances c _i(j, j′)∊[0, 1], (j, j′)=(1, 3), (2, 4) and i=1, 2,…, n are either taken from the common literature (when there is no inbreeding) as fixed constants equal to 0·5 or here are assumed to be independent and identically distributed random variables with expectation 0·5 and known variance (0·0384)². Note that the fixed value 0·5 was originally derived as an asymptotic limiting value by assuming infinitely many loci (infinite genome length). Because all the genomes are of finite length in practice, the amount of IBD-sharing varies from one pair of siblings to another. Substantial deviations may also be expected in the presence of pedigree errors in the data. The detailed motivation for such a covariance structure is provided in Appendix A.

Further, the environment effect is modelled through the cumulative effects of the locations E _i=(E _i1, E _i2, E _i3, E _i4)^t. The effects E _i are assumed to be drawn randomly from a multivariate normal distribution with mean vector zero and variance–covariance matrix Σ_ei, which has elements

(3)

where the sum is over all the locations or households shared by j and j′, D _i^l(j, j′) is the time (in years) (j, j′) spent together at location l for observation i, and σ_e² is the constant multiplier of the duration matrix D _i. Note that if (j, j′) did not share any location, then D _i^l(j, j′)=0 and the corresponding observations on E _i would be independent. This is the case for (j, j′) equal to {(1, 4), (4, 1), (3, 4), (4, 3)}. Similar identity is expected between observation (y _i2, y _i3), but in many Indian families siblings share the household even after they are married and hence the dependency between (y _i2, y _i3) and (y _i2, y _i3). This property is lacking from the settings of studies from developed countries. The rationale behind such a covariance structure is explained in Appendix A.

Excluding the interaction between G and E and under the assumption of independent and identical errors ε_i, the overall variance–covariance of y _i, viz. Σ_i is then the sum of appropriate covariance terms of the polygenic effect, shared environmental effects and the model (1) error variance (σ_m²).

(4)

where I is a 4×4 identity matrix. The mean vector of y _i is μ_i=μ+βx _i.

In twin studies the heritability is traditionally defined directly from the observed differences between the measured correlations among monozygotic and dizygotic twins (Falconer, 1989). Here, we assume no dominance variance and thus focus on ‘narrow sense’ heritability (σ_g²/σ_y²). It is clear from (4) that the phenotypic variability (σ_y²) is modelled as an age-dependent quantity (σ_g²+age σ_e²+σ_m²) and this corresponds to the diagonal elements of the Σ_i. Note that each individual has lived with oneself for the number of years which are identical to his or her current age. In the present situation, the phenotypic variation turns out to be a function of the age when the phenotype was measured.

For a given age, we define the proportion of the total variation due to polygenic effect as

(5)

and the proportion of total variation due to environmental effect as

(6)

It is clear that heritability v _g(age) is a decreasing function of age, but the rate at which it decreases with age depends on the environment variability σ_e². In general, the heritability estimates of the same trait assessed in an environment with larger variability would be lower. Similarly, v _e(age) is an increasing function of age and the rate is dependent on (σ_g²+σ_m²).

Omitting the age-dependency, a rather crude estimator of the heritability can also be given as

(7)

(8)

where is the empirical phenotypic variability, which can be estimated using sample variance of the observed phenotype data. Note that in expressions (5) and (6), σ_g² is estimated under the proposed model, while in expressions (7) and (8), is estimated externally from the observed data. The other two variances (σ_e², σ_m²) are as defined earlier under the proposed model.

3. Bayesian computation

The likelihood function is simply the product of n-independent 4-variate normal density

(9)

where θ=(μ, β, σ_g², σ_e², σ_m²) and |Σ_i| is the determinant of the matrix Σ_i. A priori independent normal priors are assigned to μ and βs, while mutually independent Gamma priors are assigned to the inverses of the three variance components. For discussions of informative and uninformative priors, see Gelman (2006) and Dongen (2006). The posterior distribution is proportional to the product of the likelihood and the prior densities

(10)

where k is the number of covariates, N(z; 0,10000) is univariate normal density with mean 0 and variance 10 000 evaluated at z and Γ(z; 0·01, 0·01) is the Gamma density with shape and scale parameters equal to 0·01 evaluated at z. The Gibbs sampler, applied to the posterior distribution of θ, involves specification of the following full conditional distributions:

where all the variables are being conditioned on their most recently sampled values. If sampling from any of these full conditional distributions is difficult, then a Metropolis–Hastings step can be used. The samples from the posterior distribution of θ (10) are generated by using the OpenBUGS 2.2.0 software (Spiegelhalter et al., 2005; Thomas et al., 2006).

Let (σ_g^2(l), σ_e^2(l), σ_m^2(l)), l=1,…, M be the Markov Chain Monte Carlo (MCMC) samples from the posterior distribution. The sample for the proportions at MCMC round l, l=l,…, M is

The posterior median and credible intervals are then estimated using these samples.

When c _i(1, 3) and c _i(2,4) are allowed to be randomly distributed according to N(0·5, 0·0384²), the above posterior distribution (10) is multiplied by

The full conditional distributions of c _i(1, 3) and c _i(2, 4) are

where

These unknown factors then get updated at each iteration. The missing data in y (under the missing at random assumption) are handled naturally in the Bayesian computation and they get updated or predicted at each iteration. The missing data in the covariates can be handled by introducing a prior distribution for covariates.

4. Data analyses

The IMS data are analysed here using three different models: (1) assuming a model y _i=μ+ε_i and a general variance–covariance matrix Σ for each i and using the Wishart prior for it, (2) the proposed model y _i=μ+βx _i+G _i+E _i+ε_i with 0·5 as the constant degree of genetic relationship, that is, C _i(1, 3) and C _i(2, 4) are equal to 0·5 for all i, and (3) the proposed robust model y _i=μ+βx _i+G _i+E _i+ε_i when the degree of genetic relationship is random and distributed according to the normal distribution with mean 0·5 and variance (0·0384)². We also estimate the proportions v _g(age) and v _e(age).

We analyse 372 (=n) household data each of which consists of four observations as shown in Figure 1. For the purpose of illustration, four phenotypes are analysed, viz. body height, weight, systolic blood pressure and total serum cholesterol. The covariates used are age when the phenotypes were measured, sex (1=man, 0=woman) and location (categorized as urban and rural). The total number of observations on the phenotypes is 1488. For cholesterol, 167 observations were missing, while for systolic blood pressure, four observations were missing. All the observations were available for height and weight. A natural question that arises is whether the correlation between phenotypes of spouses is larger for spouses that have been together for longer? To check this, a regression analysis of the phenotype of the factory worker y _i1 on the age when the phenotype was collected, the phenotype of the spouse y _i2, the duration they have lived together D _i(1, 2) and the interaction y _i2∗D _i(1, 2) gave positive regression coefficients for the phenotype of the spouse, and the duration they have lived together. Hence, it may be concluded that the spouse correlation is larger for spouses that have been together longer. We also looked at the phenotypic variances among those below and above 40 years of age to see whether phenotypic variances increase with age and by grouping the pairs of observations, ((y _i1, y _i2), (y _i1, y _i3), (y _i2, y _i4)), according to the number of years spent in 5-year groups. We noticed that in our data, the phenotypic variances among those who are below 40 years of age are smaller than those above 40 years of age. There was a clear increase in the phenotypic covariations in the initial years of living together (up to 10 years) and for the later years no clear pattern was observed. This may indicate a more general model, where the shared environmental effect decreases with age (see the Discussion section).

In Appendix B, OpenBUGS code for analysing the proposed model is given. The model (1) gave an estimate of Σ with positive covariance between individuals 1 and 4. This is not desirable since these individuals neither are genetically related nor have a shared environment. Hence, this may indicate a possible lack of fit for models (2) and (3), which assume a specific structure of dependency under random mating (see the Discussion section).

Table 1 gives the posterior median and 95% credible intervals for three variance components and for four phenotypes using models (2) and (3) and Table 2 gives corresponding estimates of the proportions of variability due to polygenic effect and environmental effect at the age of 30 years along with their crude estimates. The models show similar results except for slight differences for cholesterol. This is due to the missing data for cholesterol. It is our observation that the convergence was faster with model (3) when there were missing data for phenotypes. For all the phenotypes considered, the variation due to the shared environment is less compared to the polygenic and residual variability. This is because the environmental variation is not simply σ_e² but with a multiplier D _i. The latter are rather close for height and weight, while polygenic variability is higher than the residual variability for systolic blood pressure and cholesterol.

Table 1. Estimated variance components. Posterior median (95% credible intervals) of three variances (σ_e², environment; σ_g², genetic; σ_m², error) for various quantitative phenotypes using models with 0·5 as the constant degree of genetic relationship (2) and when the degree of genetic relationship is random (3). The abbreviation ‘SystBP’ is used for systolic blood pressure

Table 2. Estimated proportion of the total variation due to polygenic effect (v_g(age)) and due to environment effect (v_e(age)) at the age of 30 and crude estimates of these proportions ṽ_g and ṽ_e. Posterior median (95% credible intervals) of the two proportions for various quantitative phenotypes using models with 0·5 as the constant degree of genetic relationship (2) and when the degree of genetic relationship is random (3). The abbreviation ‘SystBP’ is used for systolic blood pressure

Figures 2(a)–(h) show the posterior median and 95% credible intervals for the proportions v _g(age) and v _e(age) for different phenotypes and for different values of age. The curves due to the two models are overlapping in most cases. It is clear from Figure 3(a), which shows the posterior medians of v _g(age), that the trend in the proportion is similar between weight, systolic blood pressure and cholesterol, while for height there is sudden decrease with age. This is because participants reach their maximal height in early adult life and from then onwards there is no further growth. Till age 20, the proportion of total variability due to polygenic effect is higher in systolic blood pressure compared to weight, but after the age of 30 the trend reverses. Similarly, v _e(age) are shown in Figure 3(b) where height behaves differently compared to other three phenotypes. We also carried out the analysis using log-transformation for height and the proportions of total variations due to polygenic and environmental effects were close to the results obtained without log-transformation.

Fig. 2. Posterior median (middle curves) and 95% credible intervals (upper and lower curves) of the proportion of total variation due to polygenic and shared environmental sources for various phenotypes using models with 0·5 as the constant degree of genetic relationship (2) and when the degree of genetic relationship is random (3). (a) Height – polygenic, (b) height – environmental, (c) weight – polygenic, (d) weight – environmental, (e) systolic blood pressure – polygenic, (f) systolic blood pressure – environmental, (g) total serum cholesterol – polygenic and (h) total serum cholesterol – environmental. The 2·5% limit, median and the 97·5% limit are denoted by square, polygon and triangle symbols, respectively, for model (2) and by star, cross and polygon symbols, respectively, for model (3). The vertical line corresponds to the age limit of 65 in the data.

Fig. 3. Posterior median of the proportion of total variation due to (a) polygenic and (b) shared environment sources for various phenotypes when the degree of genetic relationship is random (cross, height; polygon, weight; triangle, total serum cholesterol; square, systolic blood pressure). The vertical line corresponds to the age limit of 65 in the data.

The posterior median and 95% credible intervals are given for some of the regression parameters in Table 3. Only the sex effect is seen on height with men tending to be taller than women. Body weight is affected by age and sex. With age, body weight increases and men have a higher weight compared to women. A similar trend is seen for systolic blood pressure with respect to age and sex. The regression coefficient of sex is negative for cholesterol, indicating that women have higher cholesterol compared to men.

Table 3. Estimates of regression coefficients. Posterior median (95% credible intervals) estimates for the regression parameters (overall mean, age and male) for various quantitative phenotypes. The abbreviation ‘SystBP’ is used for systolic blood pressure

5. Discussion

The household data with a typical structure of a married couple and their siblings are analysed so that the built-in dependency structure is used in specifying the variance–covariance structure. The dependency structure allows separation of overall variability into polygenic and shared environmental components. Like McGregor et al. (2003), our analysis represents multivariate analysis of heritability in the sense that all age strata are analysed jointly compared to age-stratified heritability estimation (Brown et al., 2003). In our analyses, a general model that does not take into account the dependency structure but allows any positive definite Σ_i showed weak association between the observations (y _i1, y _i4), which may indicate a small amount of assortative mating present among the study subjects. In the proposed model, we assume that these two individuals are neither related nor have lived in the same household and the variability is directly proportional to the number of years lived in the same household. This reflects our understanding of how environmental exposures operate.

The model used in this paper is the first possible simplest model for the sibling pairs and spouse pairs association separating environmental and genetic variations. The unexplained variation is expressed in terms of the residual variance σ_m². A more complex model allowing different dependencies between the sibling pairs and spouse pairs based on the age when lived together and allowing covariate information like socio-economic status (SES) influencing the traits/phenotypes might describe the data better. One such model could be developed under the assumption that the environmental dependency is stronger if the environment is shared early in life (that is, at a younger age). For example, a more complex modelling that allows the effect of the shared environment on the phenotypes to be larger in the initial years of life and to decrease with age could be carried out using appropriate functions, for example exp {−λt). The shared environment variance is then computed as the area under the function with respect to the Gaussian process. We may also need different models for different traits since the proposed model gives counter-intuitive results for height and the development of phenotypes with age could be different. The heritability obtained without adjusting for age (Table 2, ṽ _g) is more comparable to the results reported from the twin studies, but the age-adjusted proportions are rather different. However, we must remember that modelling polygenic and environmental variations as time-dependent quantities is very rare and we must not compare mere numbers in the absolute sense. We think that the proposed approach opens up methodological questions and provides scope for further research in this area.

The estimates of heritability obtained here are not directly comparable to the ones reported in the literature. In the context of cholesterol, Mastropaolo et al. (2001) pointed out that previous studies consist of older children and adult twin pairs have indicated environmental contribution to the variation of cholesterol among individuals of 7–68% in the populations evaluated. Our results are to be compared with the previously reported results keeping this in mind. In this paper, we defined heritability as a decreasing function of age, where its value is decreasing at a rate depending on the ratio of environmental and genetic variances. Thus, it is not surprising that height seems to be less genetic at age 100 than at age 20. However, this trend may be more pronounced if the information relevant to early life environment comes from siblings who have less variability in trait and later-life information comes more from spouses who will have greater variability.

The proposed method allows a comparison of various phenotypes between rural and urban areas by incorporating an appropriately defined covariate based on migration status. A more thorough analysis focusing on the rural and urban comparisons would be reported elsewhere since the main interest here is in estimating the relative genetic and environmental variations.

In summary, we have developed a model that uses the dependency between sibling pairs and spouse pairs to estimate the age-dependent proportions of total variation that are due to polygenic and environmental effects. Using the Bayesian approach, the variation of IBD-sharing could be incorporated into the model without using marker information and this is the first time such an idea has been tried. The approximate N(0·5, 0·038²) for IBD shared by siblings is prior information based on genomewide analysis from previous studies. The Bayesian method can sample the realized IBD from the conditional posterior distribution. This explains why the variation of IBD can be incorporated without using marker information. This idea is novel and appears to be useful in separating the genetic from environmental variances. The model would benefit from testing in populations who have experienced different exposures. As pointed out to us by the Editor, using spouse information to quantify the common environmental variance may be better than using the twin data because the twin data are hard to collect, while spouse and sibling data are easier to collect. We may have very large samples for spouse–sibling data compared to the twin data.