Falconer & Mackay population mean calculation & genotypic values

We are searching data for your request:

Forums and discussions:

Manuals and reference books:

Data from registers:

Wait the end of the search in all databases.
Upon completion, a link will appear to access the found materials.

I am reading chapter 7 of Falconer and Mackay where they give a formula to calculate the population mean as a deviation from the heterozygote genotypic value.

As an example, imagine a one locus two allele trait, with alleles $A_1$ and $A_2$, with respective frequencies $p$ and $q$, and genotypic values ($ heta$) of $A_1A_1 = heta_{1,1} = 0.5$, $A_1A_2 = heta_{1,2} = 0.4$, and $A_2A_2 = heta_{2,2} = 0.3$.

The mean is given by

$ M = a_1(p-q) + 2pqd$

where $a_1$ is the deviation between the genotypic value of $A_1$ homozygotes and genotypic value heterozygotes, and, because $d=0$ (the dominance deviation), this equation reduces to

$ M = a_1(p-q)$

For my above example, if $p = 0.7$ then $M = 0.04$ and the population mean should be

$frac{1}{2} heta_{1,1} + frac{1}{2} heta_{2,2} + M = 0.44$

$frac{1}{2} heta_{1,2} + M = 0.44$

Is this correct?

What is the correct term and notation for the genotypic values (which I termed $ heta$)

Implications of the difference between true and predicted breeding values for the study of natural selection and micro-evolution

The ability to predict individual breeding values in natural populations with known pedigrees has provided a powerful tool to separate phenotypic values into their genetic and environmental components in a nonexperimental setting. This has allowed sophisticated analyses of selection, as well as powerful tests of evolutionary change and differentiation. To date, there has, however, been no evaluation of the reliability or potential limitations of the approach. In this article, I address these gaps. In particular, I emphasize the differences between true and predicted breeding values (PBVs), which as yet have largely been ignored. These differences do, however, have important implications for the interpretation of, firstly, the relationship between PBVs and fitness, and secondly, patterns in PBVs over time. I subsequently present guidelines I believe to be essential in the formulation of the questions addressed in studies using PBVs, and I discuss possibilities for future research.

Abstract

Better understanding of the genetic control of traits in breeding populations is crucial for the selection of superior varieties and parents. This study aimed to assess genetic parameters and breeding values for six essential traits in a white Guinea yam (Dioscorea rotundata Poir.) breeding population. For this, pedigree-based best linear unbiased prediction (P-BLUP) was used. The results revealed significant nonadditive genetic variances and medium to high (.45–.79) broad-sense heritability estimates for the traits studied. The pattern of associations among the genetic values of the traits suggests that selection based on a multiple-trait selection index has potential for identifying superior breeding lines. Parental breeding values predicted using progeny performance identified 13 clones with high genetic potential for simultaneous improvement of the measured traits in the yam breeding program. Subsets of progeny were identified for intermating or further variety testing based on additive genetic and total genetic values. Selection of the top 5% progenies based on the multi-trait index revealed positive genetic gains for fresh tuber yield (t ha −1 ), tuber yield (kg plant −1 ), and average tuber weight (kg). However, genetic gain was negative for tuber dry matter content and Yam mosaic virus resistance in comparison with standard varieties. Our results show the relevance of P-BLUP for the selection of superior parental clones and progenies with higher breeding values for interbreeding and higher genotypic value for variety development in yam.

Abbreviations

Results

In this section, we use empirical data and simulations of the toy-model to show that most of the heritability estimators borrowed from classical quantitative genetics are prone to significant bias, because they neglect or inaccurately model the change in resemblance between transmission partners caused by within-host evolution of the pathogen. Based on the toy-model simulations, we designate the intraclass correlation in the closest phylogenetic pairs (CPPs) and the phylogenetic heritability, H OU 2 ( t ¯ ) ⁠ , measured by the phylogenetic Ornstein–Uhlenbeck mixed model (POUMM) ( Mitov and Stadler 2016 Blanquart et al. 2017) as the most reliable estimators of pathogen trait heritability. Based on applying these estimators to a large HIV cohort, we establish a lower bound for the lg(spVL)-heritability.

Through the rest of the article, we use the symbol d_ij to denote the phylogenetic distance between two tips, i and j, on a transmission tree ( fig. 1). d_ij summarizes the total evolutionary distance between two infected hosts at the moment of measuring the trait value and is measured in substitutions per site for real trees and arbitrary time units for simulated trees. We begin our report with a result from HIV data demonstrating the relevance of within-host evolution for estimating heritability.

The lg(spVL) Correlation in HIV Phylogenetic Pairs Decreases with dij

We used one-way analysis of variance (ANOVA, r_A) and Spearman correlation (r_Sp) to estimate the correlation in phylogenetic pairs (PP) extracted from a recently published transmission tree of 8,483 HIV patients ( Hodcroft et al. 2014). As defined in Shirreff et al. (2013), phylogenetic pairs represent pairs of tips in the transmission tree that are mutually nearest to each other by phylogenetic distance (d_ij) ( fig. 1). We ordered the PPs by d_ij and split them into ten strata of equal size (deciles), evaluating the correlation between pair trait values (r_A and r_Sp) in each stratum. The point estimates and the 95% confidence intervals (CI) are shown with black and magenta points and error bars on figure 3. Dashed horizontal bars denote the 95% CI for r_A evaluated on all phylogenetic pairs. Despite some irregularities, there is a well pronounced pattern of decay in the correlation—strata to the left (small d_ij) tend to have higher r_A values than strata to the right (big d_ij). The values of r_A closely matched the values from other correlation estimators, such as DR (b) and the Pearson product mean correlation (r) (results not shown). We performed ordinary least squares regressions (OLS) of the values r A , D k and r Sp , D k on the mean phylogenetic distance, d i j , k ¯ ⁠ , in each stratum, k = 1 , … , 10 ⁠ . The slopes of both regressions were significantly negative (P<0.05) and are shown as black and magenta lines on figure 3. Similar slopes were obtained when using other stratifications of the data ( supplementary fig. S1 , Supplementary Material online).

Correlation between lg(spVL)-values in HIV phylogenetic pairs. A sample of 1917 PPs with lg(spVL)-measurements from HIV patients shows a decrease in the correlation (ICC) between pair trait values as a function of the pair phylogenetic distance d_ij. The point estimates and 95% CIs in ten strata of equal size (deciles) are depicted as points and error bars positioned at the mean d_ij for each stratum, d i j ¯ ⁠ . Black and magenta points with error-bars denote the estimated r_A and r_Sp in the real data. Dashed horizontal bars denote the 95% CI for r_A evaluated on all phylogenetic pairs. A black and a magenta inclined line denote the least squares linear regression of r_A and r_Sp on d i j ¯ ⁠ . Brown and green points with error bars denote the estimated values of r_A obtained after replacing the real trait values on the tree by values simulated under the maximum likelihood fit of the PMM and the POUMM methods, respectively (mean and 95% CI estimated from 100 replications). A brown and a green line show the expected correlation between pairs of tips at distance d_ij, as modeled under the ML-fit of the PMM and the POUMM (eqs. 2 and 3). A light-brown and a light-green region depict the 95% high posterior density (HPD) intervals inferred from Bayesian fit of the two models (Materials and Methods).

The above result shows that the value of a heritability estimator based on the correlation within phylogenetic pairs (including DR couples) depends strongly on d_ij. Another issue of all estimators of H 2 using the correlation in phylogenetic or DR pairs is that the underlying statistical methods require independence between the pairs—the trait values in one pair should not influence or be correlated with the trait values in any other pair. This assumption is not valid in general, due to the phylogenetic relationship between all patients. One way to mitigate the effects of phylogenetic relationship between pairs is to limit the analysis to the closest pairs (i.e., pairs, for which d_ij does not exceed some user specified threshold). This approach has the drawback of omitting much of the data from the analysis. As an alternative taking advantage of the entire tree, it is possible to correct for the phylogenetic relationship by using a phylogenetic comparative method (PCM). PCMs attempt to solve both of the above problems, because they 1) incorporate the branch lengths in the transmission tree to model the variance–covariance structure of the data and 2) correct for the phylogenetic correlation when estimating evolutionary parameters or the phylogenetic heritability of the trait ( Felsenstein 1985 Housworth et al. 2004 Alizon et al. 2010). These advantages of the PCMs come at the price of assuming a specific stochastic process as a model of the trait evolution along the tree. In the next subsection, we show that assuming an inappropriate process for the trait evolution can cause a significant bias in the estimate of phylogenetic heritability.

A Brownian Motion Process Cannot Reproduce the Decay of Correlation in the UK Data

We implemented a maximum likelihood and a Bayesian fit of the PMM ( Lynch 1991 Housworth et al. 2004) and its extension to an Ornstein–Uhlenbeck model of evolution (POUMM) ( Hansen 1997 Mitov and Stadler 2016 Blanquart et al. 2017). The PMM and the POUMM assume an additive model of the trait values, z ( t ) = g ( t ) + e ⁠ , in which z(t) represents the trait value at time t for a given lineage of the tree, g(t) represents a heritable (genotypic) value at time t for this lineage and e represents a nonheritable contribution summarizing the effects of the host and his/her environment on the trait and the measurement error. The only difference between the two models is their assumption about the evolution of g(t) along the branches of the tree—the PMM assumes a Brownian motion process the POUMM assumes an Ornstein–Uhlenbeck process ( Uhlenbeck and Ornstein 1930 Lande 1976 Hansen 1997).

Using the maximum likelihood estimates of the model parameters ( supplementary table S1 , Supplementary Material online), we simulated random trait trajectories on the UK tree, running 100 replications for each model. For each replication, we estimated the correlation, r_A, in PPs using the simulated values instead of the real values. The resulting correlation estimates are shown on figure 3 as brown and green points and error bars for the PMM and POUMM simulations, respectively. We notice that there is a significant difference between the correlation estimates of the two models. In particular, in the leftmost decile the POUMM estimate is significantly higher than the PMM estimate (the POUMM 95% CI excludes the PMM estimate).

The last approximation in equation (5) follows from the fact that the term exp ( − 8.35 + 36.47 d i j ) is nearly 0 for the range of phylogenetic distances ( ⁠ d i j ∈ [ 0 , 0.14 ] ⁠ ) in the UK tree (see supplementary information , Supplementary Material online, for further details on the above approximations).

Equations (4) and (5) represent a linear and an exponential model of the correlation as a function of d_ij. The values of these equations at d_ij=0 are equal to the phylogenetic heritabilities estimated at the mean root-tip distance t ¯ under PMM and POUMM (details on that later). The slope of the linear model ( eq. 4) equals −0.36 (95% HPD [−0.58, −0.21]). The rate of the exponential decay ( eq. 5) equals the POUMM parameter α=28.78 (95% HPD [16.64, 46.93]) and the half-life of decay equals ln ⁡ ( 2 ) / α = 0.02 substitutions per site (95% HPD [0.01, 0.04]).

Plotting the values of equations (4) and (5) and their 95% HPD intervals on figure 3 reveals visually that the POUMM fits better to the data than the PMM. Statistically, this is confirmed by a lower Akaike Information Criterion (AICc) for the POUMM fit and a strictly positive HPD interval for the OU parameter α ( supplementary table S1 and fig. S8, Supplementary Material online). The slope of the linear model derived from the PMM fit ( eq. 4, brown line on fig. 3) is nearly flat compared with the slopes of the two OLS fits (black and magenta lines on fig. 3). To explain this, we notice that in PMM, the covariance in phylogenetic pairs and the variance at the population level are modeled as linear functions of the root-mrca distance (t_ij) and the root-tip distance (t) (numerator and denominator in eq. 2). Importantly, both of these linear functions are bound to the same slope parameter, σ 2 . As it turns out, in the UK data, the covariance and the variance increase at different rates with respect to t_ij and t (see supplementary fig. S2 and supplementary information , Supplementary Material online). We conclude that the PMM is not an appropriate model for the correlation in phylogenetic pairs, being unable to model the above difference in the rates.

In the limit d i j → 0 ⁠ , a phylogenetic pair should be equivalent to a DR couple at the moment of transmission, that is, before the genotypes in the two hosts have diverged due to within-host evolution. Thus, it appears reasonable to use an estimate of the correlation at d_ij=0 as a proxy for the broad-sense heritability, H 2 , in the entire population. This idea has been applied in previous studies of HIV ( Hecht et al. 2010 Hollingsworth et al. 2010 Bachmann et al. 2017 Blanquart et al. 2017) as well as malaria ( Anderson et al. 2010). One potential obstacle to this approach is the possibility of introducing a sampling bias by filtering of the data. For example, if the study is on a trait, which evolves toward higher values during the course of infection, patients with lower trait values would tend to be more frequent among the CPPs than in the entire population. Thus, there is no guarantee that the trait distribution and, therefore, the heritability measured in the CPPs equals the heritability in the entire population. This problem of sampling bias affects both, resemblance-based as well as the currently used phylogenetic comparative methods. This suggests that the approach of imposing a threshold on d_ij or estimating the correlation (r_A, r_Sp or another correlation measure) at d_ij=0 needs further validation. In the next subsection, we use simulations of the toy model to show that sampling bias, although present, is comparatively small with respect to the negative bias due to measurement delay.

ANOVA-CPP and POUMM Are the Least Biased Heritability Estimators in Toy-Model Simulations

Grouping of the trait values by identical pathogen genotype. We evaluated the coefficient of determination adjusted for finite sample size, R adj 2 ⁠ , and the intraclass correlation (ICC) estimated using one-way ANOVA, r A [ id ] ⁠ . The main difference between these two estimators is the ANOVA assumption that the group-means (genotypic values) are sampled from a distribution of potentially many more genotypes than the ones found in the data. In contrast, R adj 2 assumes that all genotypes in the population are present in the sample. Since the latter assumption is true for the simulated epidemics, R adj 2 represents the reference (true) value of H 2 to which all other estimates are compared.

Known DR couples. We evaluated the regression slope of recipient on donor values in three ways: 1) b—based on the trait values at the moment of diagnosing the infection 2) b₀—based on the trait values right after the transmission events and 3) b d i j ′ —based on the subsample of diagnosed couples having d_ij not exceeding a threshold d i j ′ ⁠ . Based on a trade-off between precision and bias, we specified d i j ′ = D 1 ⁠ , D₁ denoting the first decile in the empirical distribution of d_ij (see supplementary information , Supplementary Material online).

Phylogenetic pairs (PPs) in T 10 k ⁠ . We evaluated ICC using ANOVA in three ways: 1) r_A—based on all PPs 2) r A , D 1 —based on CPPs defined as PPs in T 10 k having d_ij not exceeding the first decile, D₁ and 3) r A , 0 , lin —the estimated intercept from a linear regression of the values r A , D k on the mean values d i j , k in each decile, k = 1 , … , 10 ⁠ For the latter two estimators, which attempt to estimate r_A at d_ij=0, we use the acronym ANOVA-CPP. As an alternative to ANOVA, which is more robust to outliers (e.g., extreme values at the tails of the trait distribution), we evaluated the Spearman correlation in the first decile, hereby denoted as r Sp , D 1 ⁠ .

Transmission tree T 10 k ⁠ . We evaluated the phylogenetic heritability based on the ML fit of the PMM and POUMM models. Specifically, we compared the classical formula evaluated at the mean root-tip distance t ¯ in the tree (eqs. 10 and 12) ( Housworth et al. 2004 Leventhal and Bonhoeffer 2016) and the empirical formula based on the sample trait variance, s 2 (z) (eqs. 11 and 13) (described in Materials and Methods). For the PMM, we denote these estimators by H BM 2 ( t ¯ ) and H BM e 2 ⁠ for the POUMM, we use the symbols H OU 2 ( t ¯ ) and H OU e 2 ⁠ :

Table 1 summarizes the mathematical definition and the assumptions of the above estimators. A more detailed description of the PMM and the POUMM methods is provided in Materials and Methods. The referenced textbooks on quantitative genetics ( Lynch and Walsh 1998) are excellent references for the other methods.

Tested Estimators of the Broad-Sense Heritability of Pathogen Traits.

Input Data .	Method (Abbreviation) .	Assumptions .	Estimator .
Grouping by identical infecting strain	Adjusted coefficient of determination	The sample of data contains all genotypes present in the population	R adj 2 = 1 − N − 1 N − K s 2 ( z − G ^ ) s 2 ( z ) (6)
One-way analysis of variance (ANOVA)	Independently sampled genotypes	r A [ id ] = ( M S b − M S e ) / n ( M S b − M S e ) / n + M S e (7)
	i.i.n.d. trait-values within each group
	Equal within-group variances (homoscedasticity)
Known donor–recipient couples	Donor–recipient regression (DR)	Independently sampled donor–recipient couples
		Equal residual variance across the range of donor-values (homoscedasticity)	b = s ( z don , z rcp ) s 2 ( z don ) , (8)
		Equal donor and population variances	variants: b ⁠ , b 0 ⁠ , b d i j ′
Phylogenetic pairs (PPs)	ANOVA on all/closest PPs (ANOVA-PP, ANOVA-CPP)	ANOVA assumptions (see above)	Defined as in equation (7), but calculated on PPs
Phylogenetic pairs (PPs)	ANOVA on all/closest PPs (ANOVA-PP, ANOVA-CPP)	ANOVA assumptions (see above)	variants: r A ⁠ , r A , d i j ′
Spearman correlation on all/closest PPs	PPs are independent from one another	Pearson (product mean) correlation, calculated on the ranks of the trait-values.
Spearman correlation on all/closest PPs	PPs are independent from one another	variants: r Sp ⁠ , r Sp , d i j ′
Linear regression of r_A on d_ij upon a stratification	r_Adepends linearly on d_ij	The intercept, r A , 0 , l i n ⁠ , from the OLS fit of the model
Linear regression of r_A on d_ij upon a stratification	Equal residual variance across the range of d_ij	r A ( d i j ) = r A , 0 , l i n + ω 1 d i j . (9)
Transmission tree	Phylogenetic mixed model (PMM)	Branching BM evolution	H BM 2 ( t ¯ ) = t ¯ σ 2 / ( t ¯ σ 2 + σ e 2 ) (10)
Transmission tree	Phylogenetic mixed model (PMM)	i.i.n.d. distributed environmental deviation, e ∼ N ( 0 , σ e 2 )	H BM e 2 = 1 − σ e 2 / s 2 ( z ) (11)
Phylogenetic Ornstein–Uhlenbeck mixed model (POUMM)	Branching OU evolution	H OU 2 ( t ¯ ) = σ 2 ( 1 − exp ( − 2 α t ¯ ) ) σ 2 ( 1 − exp ( − 2 α t ¯ ) ) + 2 α σ e 2 (12)
Phylogenetic Ornstein–Uhlenbeck mixed model (POUMM)	i.i.n.d. environmental deviation, e ∼ N ( 0 , σ e 2 )	H OU e 2 = 1 − σ e 2 / s 2 ( z ) (13)

Input Data .	Method (Abbreviation) .	Assumptions .	Estimator .
Grouping by identical infecting strain	Adjusted coefficient of determination	The sample of data contains all genotypes present in the population	R adj 2 = 1 − N − 1 N − K s 2 ( z − G ^ ) s 2 ( z ) (6)
One-way analysis of variance (ANOVA)	Independently sampled genotypes	r A [ id ] = ( M S b − M S e ) / n ( M S b − M S e ) / n + M S e (7)
	i.i.n.d. trait-values within each group
	Equal within-group variances (homoscedasticity)
Known donor–recipient couples	Donor–recipient regression (DR)	Independently sampled donor–recipient couples
		Equal residual variance across the range of donor-values (homoscedasticity)	b = s ( z don , z rcp ) s 2 ( z don ) , (8)
		Equal donor and population variances	variants: b ⁠ , b 0 ⁠ , b d i j ′
Phylogenetic pairs (PPs)	ANOVA on all/closest PPs (ANOVA-PP, ANOVA-CPP)	ANOVA assumptions (see above)	Defined as in equation (7), but calculated on PPs
Phylogenetic pairs (PPs)	ANOVA on all/closest PPs (ANOVA-PP, ANOVA-CPP)	ANOVA assumptions (see above)	variants: r A ⁠ , r A , d i j ′
Spearman correlation on all/closest PPs	PPs are independent from one another	Pearson (product mean) correlation, calculated on the ranks of the trait-values.
Spearman correlation on all/closest PPs	PPs are independent from one another	variants: r Sp ⁠ , r Sp , d i j ′
Linear regression of r_A on d_ij upon a stratification	r_Adepends linearly on d_ij	The intercept, r A , 0 , l i n ⁠ , from the OLS fit of the model
Linear regression of r_A on d_ij upon a stratification	Equal residual variance across the range of d_ij	r A ( d i j ) = r A , 0 , l i n + ω 1 d i j . (9)
Transmission tree	Phylogenetic mixed model (PMM)	Branching BM evolution	H BM 2 ( t ¯ ) = t ¯ σ 2 / ( t ¯ σ 2 + σ e 2 ) (10)
Transmission tree	Phylogenetic mixed model (PMM)	i.i.n.d. distributed environmental deviation, e ∼ N ( 0 , σ e 2 )	H BM e 2 = 1 − σ e 2 / s 2 ( z ) (11)
Phylogenetic Ornstein–Uhlenbeck mixed model (POUMM)	Branching OU evolution	H OU 2 ( t ¯ ) = σ 2 ( 1 − exp ( − 2 α t ¯ ) ) σ 2 ( 1 − exp ( − 2 α t ¯ ) ) + 2 α σ e 2 (12)
Phylogenetic Ornstein–Uhlenbeck mixed model (POUMM)	i.i.n.d. environmental deviation, e ∼ N ( 0 , σ e 2 )	H OU e 2 = 1 − σ e 2 / s 2 ( z ) (13)

Note .—Notation: s 2 ( · ) ⁠ , sample variance s ( · , · ) ⁠ , sample covariance N, number of patients K, number of distinct groups of patients, that is, genotypes or phylogenetic pairs z, measured values G ^ ⁠ , estimated genotypic values: mean values from patients carrying a given genotype z don ⁠ , donor values z rcp ⁠ , recipient values M S e ⁠ , within-group mean square: M S e = ∑ ( z i − z ¯ i ) 2 N − K ⁠ , where z_i is an individual’s value and z ¯ i is the mean value of the group to which the individual belongs M S b ⁠ , among-group mean square: M S b = ∑ ( z ¯ i − z ¯ ) 2 K − 1 ⁠ , where z ¯ i is defined as above and z ¯ is the population mean value n, weighted mean number patients in a group, that is, n=2 for phylogenetic pairs and n = ( N − ∑ n i 2 N ) / ( K − 1 ) for groups of variable size α, σ, σ_e: PMM/POUMM parameters (described in Materials and Methods).

i.i.n.d., independent and identically normally distributed d_ij, phylogenetic distance between donor–recipient pairs or phylogenetic pairs d i j ′ ⁠ , threshold on d_ij (see text).

Tested Estimators of the Broad-Sense Heritability of Pathogen Traits.

Input Data .	Method (Abbreviation) .	Assumptions .	Estimator .
Grouping by identical infecting strain	Adjusted coefficient of determination	The sample of data contains all genotypes present in the population	R adj 2 = 1 − N − 1 N − K s 2 ( z − G ^ ) s 2 ( z ) (6)
One-way analysis of variance (ANOVA)	Independently sampled genotypes	r A [ id ] = ( M S b − M S e ) / n ( M S b − M S e ) / n + M S e (7)
	i.i.n.d. trait-values within each group
	Equal within-group variances (homoscedasticity)
Known donor–recipient couples	Donor–recipient regression (DR)	Independently sampled donor–recipient couples
		Equal residual variance across the range of donor-values (homoscedasticity)	b = s ( z don , z rcp ) s 2 ( z don ) , (8)
		Equal donor and population variances	variants: b ⁠ , b 0 ⁠ , b d i j ′
Phylogenetic pairs (PPs)	ANOVA on all/closest PPs (ANOVA-PP, ANOVA-CPP)	ANOVA assumptions (see above)	Defined as in equation (7), but calculated on PPs
Phylogenetic pairs (PPs)	ANOVA on all/closest PPs (ANOVA-PP, ANOVA-CPP)	ANOVA assumptions (see above)	variants: r A ⁠ , r A , d i j ′
Spearman correlation on all/closest PPs	PPs are independent from one another	Pearson (product mean) correlation, calculated on the ranks of the trait-values.
Spearman correlation on all/closest PPs	PPs are independent from one another	variants: r Sp ⁠ , r Sp , d i j ′
Linear regression of r_A on d_ij upon a stratification	r_Adepends linearly on d_ij	The intercept, r A , 0 , l i n ⁠ , from the OLS fit of the model
Linear regression of r_A on d_ij upon a stratification	Equal residual variance across the range of d_ij	r A ( d i j ) = r A , 0 , l i n + ω 1 d i j . (9)
Transmission tree	Phylogenetic mixed model (PMM)	Branching BM evolution	H BM 2 ( t ¯ ) = t ¯ σ 2 / ( t ¯ σ 2 + σ e 2 ) (10)
Transmission tree	Phylogenetic mixed model (PMM)	i.i.n.d. distributed environmental deviation, e ∼ N ( 0 , σ e 2 )	H BM e 2 = 1 − σ e 2 / s 2 ( z ) (11)
Phylogenetic Ornstein–Uhlenbeck mixed model (POUMM)	Branching OU evolution	H OU 2 ( t ¯ ) = σ 2 ( 1 − exp ( − 2 α t ¯ ) ) σ 2 ( 1 − exp ( − 2 α t ¯ ) ) + 2 α σ e 2 (12)
Phylogenetic Ornstein–Uhlenbeck mixed model (POUMM)	i.i.n.d. environmental deviation, e ∼ N ( 0 , σ e 2 )	H OU e 2 = 1 − σ e 2 / s 2 ( z ) (13)

Input Data .	Method (Abbreviation) .	Assumptions .	Estimator .
Grouping by identical infecting strain	Adjusted coefficient of determination	The sample of data contains all genotypes present in the population	R adj 2 = 1 − N − 1 N − K s 2 ( z − G ^ ) s 2 ( z ) (6)
One-way analysis of variance (ANOVA)	Independently sampled genotypes	r A [ id ] = ( M S b − M S e ) / n ( M S b − M S e ) / n + M S e (7)
	i.i.n.d. trait-values within each group
	Equal within-group variances (homoscedasticity)
Known donor–recipient couples	Donor–recipient regression (DR)	Independently sampled donor–recipient couples
		Equal residual variance across the range of donor-values (homoscedasticity)	b = s ( z don , z rcp ) s 2 ( z don ) , (8)
		Equal donor and population variances	variants: b ⁠ , b 0 ⁠ , b d i j ′
Phylogenetic pairs (PPs)	ANOVA on all/closest PPs (ANOVA-PP, ANOVA-CPP)	ANOVA assumptions (see above)	Defined as in equation (7), but calculated on PPs
Phylogenetic pairs (PPs)	ANOVA on all/closest PPs (ANOVA-PP, ANOVA-CPP)	ANOVA assumptions (see above)	variants: r A ⁠ , r A , d i j ′
Spearman correlation on all/closest PPs	PPs are independent from one another	Pearson (product mean) correlation, calculated on the ranks of the trait-values.
Spearman correlation on all/closest PPs	PPs are independent from one another	variants: r Sp ⁠ , r Sp , d i j ′
Linear regression of r_A on d_ij upon a stratification	r_Adepends linearly on d_ij	The intercept, r A , 0 , l i n ⁠ , from the OLS fit of the model
Linear regression of r_A on d_ij upon a stratification	Equal residual variance across the range of d_ij	r A ( d i j ) = r A , 0 , l i n + ω 1 d i j . (9)
Transmission tree	Phylogenetic mixed model (PMM)	Branching BM evolution	H BM 2 ( t ¯ ) = t ¯ σ 2 / ( t ¯ σ 2 + σ e 2 ) (10)
Transmission tree	Phylogenetic mixed model (PMM)	i.i.n.d. distributed environmental deviation, e ∼ N ( 0 , σ e 2 )	H BM e 2 = 1 − σ e 2 / s 2 ( z ) (11)
Phylogenetic Ornstein–Uhlenbeck mixed model (POUMM)	Branching OU evolution	H OU 2 ( t ¯ ) = σ 2 ( 1 − exp ( − 2 α t ¯ ) ) σ 2 ( 1 − exp ( − 2 α t ¯ ) ) + 2 α σ e 2 (12)
Phylogenetic Ornstein–Uhlenbeck mixed model (POUMM)	i.i.n.d. environmental deviation, e ∼ N ( 0 , σ e 2 )	H OU e 2 = 1 − σ e 2 / s 2 ( z ) (13)

i.i.n.d., independent and identically normally distributed d_ij, phylogenetic distance between donor–recipient pairs or phylogenetic pairs d i j ′ ⁠ , threshold on d_ij (see text).

Within: neutral/Between: neutral

Within: select/Between: neutral

Within: neutral/Between: select

Within: select/Between: select

For each of these scenarios and mean contact interval 1 / κ ∈ < 2 , 4 , 6 , 8 , 10 , 12 >(arbitrary time units), we executed ten simulations resulting in a total of 4 × 6 × 10 = 240 simulations. Of the 240 simulations, 175 resulted in epidemic outbreaks of at least 10,000 diagnosed hosts. For each outbreak, we analyzed the populations of the first up to 10,000 diagnosed hosts.

Rarer transmission events (bigger 1 / κ ⁠ ) result in longer transmission trees and, therefore, longer average phylogenetic distance between tips, d_ij ( supplementary fig. S3 , Supplementary Material online). This enabled demonstrating the effect of accumulating within-host evolution on the different heritability estimators ( fig. 4).

Heritability estimates in toy-model simulations. (A–D) H 2 -estimates in simulations of “neutral” and “select” within-/between-host dynamics. Each group of box-whiskers summarizes the simulations for a fixed scenario and contact interval, 1 / κ ⁠ white boxes (background) denote true heritability, colored boxes denote estimates (foreground). Statistical significance is evaluated through t-tests summarized in table 2.

Figure 4 shows that the estimators b D 1 ⁠ , b ⁠ , r A , D 1 ⁠ , and r A are negatively biased in general for all toy-model scenarios. This bias tends to increase with the mean contact interval, 1 / κ (respectively, d_ij), because random within-host mutation tends to decrease the genetic overlap between DRs and phylogenetic pairs ( supplementary fig. S4 , Supplementary Material online). The negative bias was far less pronounced when imposing a threshold on d_ij but this came at the cost of precision (less biased but longer box-whisker plots for b D 1 and r A , D 1 compared with b and r_A) ( fig. 4). Several additional sources of bias were revealed when considering the practically unavailable estimators b₀ and r A [ i d ] ⁠ . The estimator r A [ i d ] was positively biased due to the small number of simulated genotypes (only six)—this was validated through additional simulations showing that r A [ i d ] converges to the true value for a slightly bigger number of genotypes (e.g., K≥24 genotypes, see supplementary information , Supplementary Material online). The estimator b₀ was behaving accurately in the neutral/neutral scenario (excluding very short contact intervals) but tended to have a bias in both directions in all scenarios involving selection. The main reason for these biases was the phenomenon of “sampling bias” consisting in a difference between the distributions of measured values in the DR couples and the population of interest. Although its magnitude was comparatevely small in the simulations, we presume that sampling bias could play an important role in real biological applications. We already gave an example of this bias in the previous subsection. Another manifestation of sampling bias is the fact that b₀ does not fully eliminate the effect of within host-evolution (and selection) in the donors. This is why, in cases of selection, the phenotypic variance in the donors tends to be smaller than the variance in the recipients as well as the variance in the population ( supplementary fig. S5 , Supplementary Material online). Additional details on these potential sources of bias are provided in supplementary information , Supplementary Material online.

Further, the simulations showed that a worsening fit of the BM model on longer transmission trees was causing an inflated estimate of the environmental deviation, σ e ⁠ , in the PMM fits and, therefore, a negative bias in H BM 2 ( t ¯ ) and H BM e 2 (compare estimates for small and big values of 1 / κ on fig. 4 and supplementary fig. S6 C, Supplementary Material online). In contrast with the PMM, the POUMM estimates, H OU 2 ( t ¯ ) and H OU e 2 were far more accurate and the value of σ e in the POUMM ML fit was nearly matching the true nonheritable deviation in most simulations ( fig. 4 and supplementary fig. S6 C, Supplementary Material online). The better ML fit of the POUMM was confirmed by stronger statistical support, namely by lower AICc values in all toy-model simulations ( supplementary fig. S6 D, Supplementary Material online).

The fact that the POUMM outperformed the PMM in all scenarios contradicted with the initial belief that the PMM should be the better suited model for a neutrally evolving trait represented by the neutral/neutral scenario, whereas the POUMM should fit better to scenarios involving selection. It was also counterintuitive that the inferred parameter α from the POUMM model was significantly positive in all simulations including the neutral/neutral scenario ( supplementary fig. S6 B, Supplementary Material online). To better understand this phenomenon, we performed the PP stratification analysis on the toy-model data ( supplementary fig. S7 , Supplementary Material online). This revealed a pattern of correlation that decays exponentially with d_ij. The shape of exponential decay was mostly pronounced for longer contact intervals, 1 / κ ⁠ , particularly in the neutral/neutral scenario (first column on supplementary fig. S7 , Supplementary Material online). In supplementary information , Supplementary Material online, we show that an exponentially decaying phenotypic correlation is consistent with a neutrally mutating genotype under a Jukes–Cantor substitution model ( Yang 2006). The decay of the correlation was still present in scenarios involving within- and/or between-host selection but the observed pattern was rather irregular and deviating from an exponential function of d_ij ( supplementary fig. S7 , Supplementary Material online). In most cases, the ML fit of the PMM method was a bad fit to the decay of correlation (brown dots and error-bars on supplementary fig. S7 , Supplementary Material online) for longer contact intervals, there was a tendency toward constant values of the correlation under PMM far below the true value (brown dots and error bars on supplementary fig. S7 , Supplementary Material online). This explains the overall better accuracy of the POUMM versus the PMM method.

Table 2 shows the average bias of each tested estimator for each of the four scenarios. We conclude that, apart from the practically inaccessible estimators based on grouping by identical genotype ( ⁠ R adj 2 and r A [ id ] ⁠ ), the most accurate estimators of H 2 in the toy-model simulations are H OU 2 ( t ¯ ) and H OU e 2 followed by estimators of the correlation in PPs minimizing the phylogenetic distance d_ij, that is ( ⁠ r A , D 1 ⁠ , r A , 0 , l i n ⁠ , r Sp , D 1 ⁠ ). In the next subsection, we report the results from these estimators in the UK HIV data.

Mean Difference H 2 ̂ − R adj 2 from the Toy-Model Simulations Grouped by Scenario.

Within: .	Neutral .	Neutral .	Select .	Select .
Between: .	Neutral .	Select .	Neutral .	Select .
N	50	41	47	37
b 0	−0.01 *	−0.02 **	0.05 **	0.04 **
b D 1	−0.07 **	−0.04 **	0	−0.01
b	−0.25 **	−0.2 **	−0.07 **	−0.06 **
r A [ i d ]	0.05 **	0.05 **	0.08 **	0.06 **
r A , 0 , l i n ^	−0.05 **	−0.06 **	0.01	−0.04 **
r A , D 1	−0.05 **	−0.06 **	0	−0.03 *
r A	−0.18 **	−0.15 **	−0.06 **	−0.08 **
r Sp , D 1	−0.05 **	−0.05 **	−0.05 **	−0.07 **
H BM 2 ( t ¯ )	−0.17 **	−0.17 **	−0.01	−0.04 *
H BM e 2	−0.28 **	−0.24 **	−0.12 **	−0.16 **
H OU 2 ( t ¯ )	−0.01	−0.02 **	0.01 *	0.03 **
H OU e 2	−0.01	−0.02 **	0.01 *	0.03 **

Within: .	Neutral .	Neutral .	Select .	Select .
Between: .	Neutral .	Select .	Neutral .	Select .
N	50	41	47	37
b 0	−0.01 *	−0.02 **	0.05 **	0.04 **
b D 1	−0.07 **	−0.04 **	0	−0.01
b	−0.25 **	−0.2 **	−0.07 **	−0.06 **
r A [ i d ]	0.05 **	0.05 **	0.08 **	0.06 **
r A , 0 , l i n ^	−0.05 **	−0.06 **	0.01	−0.04 **
r A , D 1	−0.05 **	−0.06 **	0	−0.03 *
r A	−0.18 **	−0.15 **	−0.06 **	−0.08 **
r Sp , D 1	−0.05 **	−0.05 **	−0.05 **	−0.07 **
H BM 2 ( t ¯ )	−0.17 **	−0.17 **	−0.01	−0.04 *
H BM e 2	−0.28 **	−0.24 **	−0.12 **	−0.16 **
H OU 2 ( t ¯ )	−0.01	−0.02 **	0.01 *	0.03 **
H OU e 2	−0.01	−0.02 **	0.01 *	0.03 **

Note .—Statistical significance is estimated by Student’s t-tests, P values denoted by an asterisk as follows: * P<0.01 **P<0.001. Gray background indicates estimates that are unavailable in practice.

Mean Difference H 2 ̂ − R adj 2 from the Toy-Model Simulations Grouped by Scenario.

Within: .	Neutral .	Neutral .	Select .	Select .
Between: .	Neutral .	Select .	Neutral .	Select .
N	50	41	47	37
b 0	−0.01 *	−0.02 **	0.05 **	0.04 **
b D 1	−0.07 **	−0.04 **	0	−0.01
b	−0.25 **	−0.2 **	−0.07 **	−0.06 **
r A [ i d ]	0.05 **	0.05 **	0.08 **	0.06 **
r A , 0 , l i n ^	−0.05 **	−0.06 **	0.01	−0.04 **
r A , D 1	−0.05 **	−0.06 **	0	−0.03 *
r A	−0.18 **	−0.15 **	−0.06 **	−0.08 **
r Sp , D 1	−0.05 **	−0.05 **	−0.05 **	−0.07 **
H BM 2 ( t ¯ )	−0.17 **	−0.17 **	−0.01	−0.04 *
H BM e 2	−0.28 **	−0.24 **	−0.12 **	−0.16 **
H OU 2 ( t ¯ )	−0.01	−0.02 **	0.01 *	0.03 **
H OU e 2	−0.01	−0.02 **	0.01 *	0.03 **

Within: .	Neutral .	Neutral .	Select .	Select .
Between: .	Neutral .	Select .	Neutral .	Select .
N	50	41	47	37
b 0	−0.01 *	−0.02 **	0.05 **	0.04 **
b D 1	−0.07 **	−0.04 **	0	−0.01
b	−0.25 **	−0.2 **	−0.07 **	−0.06 **
r A [ i d ]	0.05 **	0.05 **	0.08 **	0.06 **
r A , 0 , l i n ^	−0.05 **	−0.06 **	0.01	−0.04 **
r A , D 1	−0.05 **	−0.06 **	0	−0.03 *
r A	−0.18 **	−0.15 **	−0.06 **	−0.08 **
r Sp , D 1	−0.05 **	−0.05 **	−0.05 **	−0.07 **
H BM 2 ( t ¯ )	−0.17 **	−0.17 **	−0.01	−0.04 *
H BM e 2	−0.28 **	−0.24 **	−0.12 **	−0.16 **
H OU 2 ( t ¯ )	−0.01	−0.02 **	0.01 *	0.03 **
H OU e 2	−0.01	−0.02 **	0.01 *	0.03 **

Heratibality of lg(spVL) in the UK HIV Cohort

We evaluated the correlation in the CPPs (ANOVA and Spearman correlation) in data from the UK HIV cohort comprising lg(spVL) measurements and a tree of viral (pol) sequences from 8,483 patients inferred previously in ( Hodcroft et al. 2014). In addition, we performed a Bayesian fit of the POUMM and the PMM methods to the same data. The goal was to test our conclusions on a real data set and to compare the H 2 -estimates from CPPs and POUMM to previous PMM/ReML-estimates on exactly the same data ( Hodcroft et al. 2014).

In applying ANOVA-CPP, the first step has been to define the threshold phylogenetic distance for defining CPPs. To that end, we explored different stratifications of the PPs as shown on supplementary figure S1 B, Supplementary Material online, and a scatter plot of the phylogenetic distances against the absolute phenotypic differences, | Δ lg ⁡ ( spVL ) | ( fig. 5A). This revealed a small set of 116 PPs having d i j ≤ 10 − 4 and narrowly coinciding with the first vigintile (also called 20-quantile or ventile) of d_ij. The phylogenetic distance in all remaining tip-pairs was more than an order of magnitude bigger, that is, d i j > 10 − 3 ⁠ . Given that the phylogenetic distance on the transmission tree is measured in substitutions per site and the length of the pol-region is in the order of 10 3 sites, we presume that the above set of 116 PPs corresponds to a set of 116 pairs of identical pol consensus sequences (no sequence data were available to check this). Based on this observation, we defined the above pairs as CPPs and the threshold was formally set to d i j ′ = 10 − 4 ⁠ . We validated that the CPPs were randomly distributed along the tree ( fig. 5B). The random distribution of the CPPs along the transmission tree suggests that these phylogenetic pairs correspond to randomly occurring early detections of infection (trait values from each pair depicted as magenta segments on fig. 5B). To check that the filtering of the data, did not introduce a considerable sampling bias due to selection (see previous subsection), we also validated that there was no substantial difference in the trait distributions of all patients, the PPs and the CPPs ( fig. 5C).

Phylogenetic pairs in lg(spVL) data from the United Kingdom. (A) A scatter plot of the phylogenetic distances between pairs of tips against their absolute phenotypic differences: gray, PPs ( ⁠ d i j > 10 − 4 ⁠ ) magenta, CPPs ( ⁠ d i j < 10 − 4 ⁠ ). A black line shows the linear regression of | Δ lg ⁡ ( spVL ) | on d_ij (the slope of the regression was statistically positive at the 0.01 level). (B) A box-plot representing the trait-distribution along the transmission tree. Each box-whisker represents the lg(spVL)-distribution of patients grouped by their distance from the root of the tree measured in substitutions per site. Wider boxes indicate groups bigger in size. Segments in magenta denote lg(spVL)-values in CPPs. (C) A box-plot of the lg(spVL)-distribution in all patients (black), PPs (gray), and CPPs (magenta).

ANOVA-CPPs ( ⁠ r A , D 1 ⁠ , r A , 10 − 4 ⁠ , r A , V 1 ⁠ ) and the original PP-method r_A

The intercept from the linear regression of r_A on d_ij upon a stratification of the PPs into deciles ( ⁠ r A , 0 , l i n ⁠ , eq. 9)

Spearman correlatoin in CPPs ( ⁠ r Sp , D 1 ⁠ , r Sp , 10 − 4 ⁠ , r Sp , V 1 ⁠ ) and in all PPs ( ⁠ r Sp ⁠ ).

The intercept from the linear regression of r Sp on d_ij upon a stratification of the PPs into deciles ( ⁠ r Sp , 0 , l i n ⁠ )

POUMM ( ⁠ H OU 2 ( t ¯ ) ⁠ , H OU e 2 ⁠ ), versus PMM ( ⁠ H BM 2 ( t ¯ ) ⁠ , H BM e 2 ⁠ ) on the entire tree

The results from these analyses are reported in table 3. ANOVA- and Spearman-correlation estimates, which minimized the phylogenetic distance by means of regression or filtering of the phylogenetic pairs had point-estimates of r A , 10 − 4 = 0.17 and r Sp , 10 − 4 = 0.22 ⁠ . The slightly higher estimate for the Spearman correlation could be explained by the presence of outliers in the data. Applying the POUMM to the entire tree reported a point estimate H OU 2 ( t ¯ ) = 0.21 (8,483 patients, 95% CI [0.14, 0.29]).

Estimates of lg(spVL)-Heritability in HIV Data from the United Kingdom.

Method .	N .	H ^ 2 .	95% CI .	95% HPD .
Linear regression of r_A on d i j ¯ in deciles (eq. 9) ( ⁠ r A , 0 , l i n ⁠ )	10 points	0.17	[0.09, 0.24]	–
Linear regression of r_Sp on d i j ¯ in deciles ( ⁠ r Sp , 0 , l i n ⁠ )	10 points	0.18	[0.11, 0.25]	–
ANOVA-CPP ( ⁠ r A , V 1 ⁠ )	224	0.17	[−0.02, 0.31]	–
ANOVA-CPP ( ⁠ r A , 10 − 4 ⁠ )	232	0.16	[0.01, 0.30]	–
ANOVA-CPP ( ⁠ r A , D 1 ⁠ )	384	0.16	[0.06, 0.25]	–
ANOVA-PP (r_A) a	3,834	0.11	[0.08, 0.14]	–
Spearman-CPP ( ⁠ r Sp , V 1 ⁠ )	224	0.23	[0.05, 0.42]	–
Spearman-CPP ( ⁠ r Sp , 10 − 4 ⁠ )	232	0.22	[0.03, 0.4]	–
Spearman-CPP ( ⁠ r Sp , D 1 ⁠ )	384	0.2	[0.06, 0.34]	–
Spearman-PP (r_Sp) a	3,834	0.11	[0.06, 0.15]	–
POUMM ( ⁠ H OU 2 ( t ¯ ) ⁠ )	8,483	0.21	–	[0.14, 0.29]
POUMM ( ⁠ H OU e 2 ⁠ )	8,483	0.2	–	[0.13, 0.29]
PMM ( ⁠ H BM 2 ( t ¯ ) ⁠ ) b	8,483	0.08	–	[0.05, 0.12]
PMM ( ⁠ H BM e 2 ⁠ ) b	8,483	0.06	–	[0.02, 0.1]
PMM, ReML ( Hodcroft et al. 2014) b	8,483	0.06	[0.03, 0.09]	–

Method .	N .	H ^ 2 .	95% CI .	95% HPD .
Linear regression of r_A on d i j ¯ in deciles (eq. 9) ( ⁠ r A , 0 , l i n ⁠ )	10 points	0.17	[0.09, 0.24]	–
Linear regression of r_Sp on d i j ¯ in deciles ( ⁠ r Sp , 0 , l i n ⁠ )	10 points	0.18	[0.11, 0.25]	–
ANOVA-CPP ( ⁠ r A , V 1 ⁠ )	224	0.17	[−0.02, 0.31]	–
ANOVA-CPP ( ⁠ r A , 10 − 4 ⁠ )	232	0.16	[0.01, 0.30]	–
ANOVA-CPP ( ⁠ r A , D 1 ⁠ )	384	0.16	[0.06, 0.25]	–
ANOVA-PP (r_A) a	3,834	0.11	[0.08, 0.14]	–
Spearman-CPP ( ⁠ r Sp , V 1 ⁠ )	224	0.23	[0.05, 0.42]	–
Spearman-CPP ( ⁠ r Sp , 10 − 4 ⁠ )	232	0.22	[0.03, 0.4]	–
Spearman-CPP ( ⁠ r Sp , D 1 ⁠ )	384	0.2	[0.06, 0.34]	–
Spearman-PP (r_Sp) a	3,834	0.11	[0.06, 0.15]	–
POUMM ( ⁠ H OU 2 ( t ¯ ) ⁠ )	8,483	0.21	–	[0.14, 0.29]
POUMM ( ⁠ H OU e 2 ⁠ )	8,483	0.2	–	[0.13, 0.29]
PMM ( ⁠ H BM 2 ( t ¯ ) ⁠ ) b	8,483	0.08	–	[0.05, 0.12]
PMM ( ⁠ H BM e 2 ⁠ ) b	8,483	0.06	–	[0.02, 0.1]
PMM, ReML ( Hodcroft et al. 2014) b	8,483	0.06	[0.03, 0.09]	–

Note .—Also written are the results from a previous analysis on the same data set ( Hodcroft et al. 2014). “–”: the analysis was not done in the mentioned study. Gray background: estimates considered unreliable due to: a negative bias caused by measurement delays and b negative bias caused by BM violation. Uncertainty in the estimates is expressed in terms of 95% confidence intervals (CI), or, in the case of Bayesian inference, by 95% high posterior density intervals (HPDs).

Estimates of lg(spVL)-Heritability in HIV Data from the United Kingdom.

Method .	N .	H ^ 2 .	95% CI .	95% HPD .
Linear regression of r_A on d i j ¯ in deciles (eq. 9) ( ⁠ r A , 0 , l i n ⁠ )	10 points	0.17	[0.09, 0.24]	–
Linear regression of r_Sp on d i j ¯ in deciles ( ⁠ r Sp , 0 , l i n ⁠ )	10 points	0.18	[0.11, 0.25]	–
ANOVA-CPP ( ⁠ r A , V 1 ⁠ )	224	0.17	[−0.02, 0.31]	–
ANOVA-CPP ( ⁠ r A , 10 − 4 ⁠ )	232	0.16	[0.01, 0.30]	–
ANOVA-CPP ( ⁠ r A , D 1 ⁠ )	384	0.16	[0.06, 0.25]	–
ANOVA-PP (r_A) a	3,834	0.11	[0.08, 0.14]	–
Spearman-CPP ( ⁠ r Sp , V 1 ⁠ )	224	0.23	[0.05, 0.42]	–
Spearman-CPP ( ⁠ r Sp , 10 − 4 ⁠ )	232	0.22	[0.03, 0.4]	–
Spearman-CPP ( ⁠ r Sp , D 1 ⁠ )	384	0.2	[0.06, 0.34]	–
Spearman-PP (r_Sp) a	3,834	0.11	[0.06, 0.15]	–
POUMM ( ⁠ H OU 2 ( t ¯ ) ⁠ )	8,483	0.21	–	[0.14, 0.29]
POUMM ( ⁠ H OU e 2 ⁠ )	8,483	0.2	–	[0.13, 0.29]
PMM ( ⁠ H BM 2 ( t ¯ ) ⁠ ) b	8,483	0.08	–	[0.05, 0.12]
PMM ( ⁠ H BM e 2 ⁠ ) b	8,483	0.06	–	[0.02, 0.1]
PMM, ReML ( Hodcroft et al. 2014) b	8,483	0.06	[0.03, 0.09]	–

Method .	N .	H ^ 2 .	95% CI .	95% HPD .
Linear regression of r_A on d i j ¯ in deciles (eq. 9) ( ⁠ r A , 0 , l i n ⁠ )	10 points	0.17	[0.09, 0.24]	–
Linear regression of r_Sp on d i j ¯ in deciles ( ⁠ r Sp , 0 , l i n ⁠ )	10 points	0.18	[0.11, 0.25]	–
ANOVA-CPP ( ⁠ r A , V 1 ⁠ )	224	0.17	[−0.02, 0.31]	–
ANOVA-CPP ( ⁠ r A , 10 − 4 ⁠ )	232	0.16	[0.01, 0.30]	–
ANOVA-CPP ( ⁠ r A , D 1 ⁠ )	384	0.16	[0.06, 0.25]	–
ANOVA-PP (r_A) a	3,834	0.11	[0.08, 0.14]	–
Spearman-CPP ( ⁠ r Sp , V 1 ⁠ )	224	0.23	[0.05, 0.42]	–
Spearman-CPP ( ⁠ r Sp , 10 − 4 ⁠ )	232	0.22	[0.03, 0.4]	–
Spearman-CPP ( ⁠ r Sp , D 1 ⁠ )	384	0.2	[0.06, 0.34]	–
Spearman-PP (r_Sp) a	3,834	0.11	[0.06, 0.15]	–
POUMM ( ⁠ H OU 2 ( t ¯ ) ⁠ )	8,483	0.21	–	[0.14, 0.29]
POUMM ( ⁠ H OU e 2 ⁠ )	8,483	0.2	–	[0.13, 0.29]
PMM ( ⁠ H BM 2 ( t ¯ ) ⁠ ) b	8,483	0.08	–	[0.05, 0.12]
PMM ( ⁠ H BM e 2 ⁠ ) b	8,483	0.06	–	[0.02, 0.1]
PMM, ReML ( Hodcroft et al. 2014) b	8,483	0.06	[0.03, 0.09]	–

Conversely, the heritability estimates from the original PP method (ANOVA or Spearman correlatoin on all PPs) and the PMM were significantly lower and falling below the 95% CIs from the POUMM ( table 3). This confirms the observation from the toy-model simulations that these estimators are negatively biased, since they ignore or inaccurately model the changing correlation within pairs of tips. We validated the stronger statistical support for the POUMM with respect to the PMM, by its lower AICc value ( supplementary table S1 , Supplementary Material online) and by the posterior density for the POUMM parameter α ( supplementary fig. S8 , Supplementary Material online).

Finally, we compared our estimates of lg(spVL)-heritability to previous applications of the same methods on different data sets ( fig. 6). In agreement with the toy-model simulations, estimates of H 2 using PMM or other BM-based phylogenetic methods (i.e., Blomberg’s K and Pagel’s λ) are notably lower than all other estimates, suggesting that these phylogenetic comparative methods underestimate H 2 resemblance-based estimates are down-biased by measurement delays (e.g., compare early vs. late in the Netherlands on fig. 6).

A comparison between H 2 H 2 -estimates from the UK HIV-cohort and previous estimates on African, Swiss, and Dutch data. (A) Estimates with minimized measurement delay (dark cadet-blue) and POUMM estimates (green) (B) Down-biased estimates due to higher measurement delays (light-blue) or violated BM-assumption (brown). Confidence is depicted either as segments indicating estimated 95% CI or P values in cases of missing 95% CIs. References to the corresponding publications are written as numbers in superscript as follows: 1: Tang et al. (2004) 2: Hecht et al. (2010) 3: Hollingsworth et al. (2010) 4: van der Kuyl et al. (2010) 5: Lingappa et al. (2013) 6: Yue et al. (2013) 7: Alizon et al. (2010) 8: Shirreff et al. (2013) 9: Hodcroft et al. (2014) 10: Blanquart et al. (2017) 11: Bertels et al. (2018) 12: this work. For clarity, estimates from previous studies, which are not directly comparable (e.g., previous results from Swiss MSM/strict data sets Alizon et al. 2010).

In summary, POUMM and ANOVA-CPP yield agreeing estimates for H 2 in the UK data and these estimates agree with resemblance-based estimates in data sets with short measurement delay (different African countries and the Netherlands). Similar to the toy-model simulations, we notice a well-pronounced pattern of negative bias for the other estimators, PMM and ANOVA-PP, as well as for the previous resemblance-based studies on data with long measurement delay.

Inbreeding and Imprinting

An interesting aspect of the above variance decompositions is the similarity between inbreeding and imprinting, as inbreeding also introduces a covariance between additive and dominance effects (Harris 1964) that may not be partitioned if an incorrect method is used. To investigate this similarity, we incorporate inbreeding into Approach 2b. We represent an inbred population by dividing the population into two groups: a group that represents the expected Hardy-Weinberg proportions, comprising an overall proportion of (1 − f), and a completely homozygous group with no heterozygotes, comprising a proportion f of the population. Thus genotypic frequencies for A₁A₁, A₂A₁ and A₁A₂, and A₂A₂ genotypes are each and , respectively. Now the overall population mean incorporating both inbreeding (I) and imprinting is

When there is no inbreeding, f = 0, the population is in Hardy-Weinberg proportions, and the mean reduces to as expected. With no imprinting, the mean reduces to . We assume that the inbreeding coefficient f is stable across generations, so the proportion of heterozygotes does not change.

As in Approach 2b, male and female additive and dominance deviations may be calculated separately (Santure and Spencer 2006). For example, the additive effect of inheriting an A₁ allele maternally is The remaining additive effects are

Breeding values and dominance deviations may be calculated as in Approach 2b.

Genetic variance components

The total variance for an inbred population with imprinting is where is the total genetic variance for the case of imprinting only (2). When there is complete inbreeding (f = 1), the total variance is (6) and for no inbreeding (f = 0), we recover . (7) The total variance may also be rewritten as (8) and for the case of no imprinting (k₁ = k₂ = k), the total genetic variance becomes (9) (Harris, 1964), where , (10) and

Comparing equations (6) and (7), we can see that in the absence of imprinting and dominance (k₁ = k₂ = f = 0), the total variance for an inbred population is twice that of an outbred population (k₁ = k₂ = 0, f = 1), and indeed the effect of inbreeding is linear an inbreeding coefficient of yields a total variance times the variance with no inbreeding. However, for dominance but no imprinting (9), the effect of increasing inbreeding is nonlinear the population variance may increase or decrease relative to an outbred population depending on the allele frequencies and the value of the dominance coefficient. A similar effect is evident with imprinting if there is no dominance (i.e., k₁ = −k₂), the population variance linearly increases with increasing inbreeding, while with both dominance and imprinting the population variance is a quadratic function of f. Thus, it is dominance but not imprinting which determines the relationship of the population variance with increasing inbreeding.

For a highly selfing species, the degree of imprinting may have a large effect on the total population variance. For example, consider a population with and . Setting k₁ = −k₂ (imprinting, no dominance), the total variance is 0.20 for and 0.25 for . Interestingly, the effect of imprinting becomes less pronounced as inbreeding levels increase for , the total variance increases from 0.18 for to 0.25 for , while for the variance increases from 0.23 for to 0.25 for .

The female and male additive variances are (11) and (12) where and are the female and male additive variances calculated from Approaches 1 and 2b (Table 3). The female and male dominance variances are (13) and (14) where is the dominance variance (Approaches 1 and 2b, Table 3). Interestingly, and unlike the pure imprinting case, the dominance variance is different for males and females. The female and male covariances between additive and dominance terms are (15) and (16) where and are the female and male covariances between additive and dominance effects (Approaches 1 and 2b, Table 3). When there is no imprinting, the female and male covariances reduce to the same value:

Considering that we can see that and hence the covariance is strictly negative under inbreeding alone. Recall that under imprinting alone and (Table 3). Now we may rearrange to give , so that and hence , so the average of male and female covariances under inbreeding is also strictly negative. However, if k₁ and k₂ are of opposite sign, then one of or may be positive. Thus, although both imprinting and inbreeding introduce a covariance between additive and dominance effects, it is only the presence of imprinting that allows the covariance in one sex to be positive. Imprinting can therefore have a significant effect on the total genetic variance and on the sex-specific components of variance of an inbred population.

05tFmQ2,tf%fP&X1`[email protected]q765G F1DDKuVVqH>8RohB%5t3V8,d:Lt5Ut_q/FY [email protected]#,=.SgJBM^7nSqa7dOgFZj^0ed^Lf apS)Eoeq?5IC2U(9_5m'c3mDHafeLaPDZ[Q+Dh+ N,t0B/(+,W)5N5'm0pmOls_U0H(rE4eW!S+bjEp,5K,2F7U6,JOF,[email protected]@2 3tEYeR.MDsBc(Fr%dZfbm7GPVCl =][0W #f`X" @ZK``'g 4,fpdBkUL)/ [Fe_&&Oh#TsSUoGo3>>EZ3j=Ki9ZD_S#>%=+Y!2-[YAf3f6U.[6m [email protected]_:$X+23Ep[oT86`2QrNa(j[KR=]SqZB.6s]V.pGt[/Dr 1OTKF MQ1q8EmMZ4fp=tGQ"se>i!rtl3?nh%2j&!iK(&[gS9iZ0248 @e(JqKEQVDiE.Z&kk0S=:0cu1"#[email protected]%_5M-CGeB8"Q-["[email protected]!M8 'nN$T$tGYk>@E30?`Xgkg0gR%%o4pr#57VQ VFpN6/)1T6c'UY+4T`eX(SoTklQX)iKo,BJcVR0 L5?DdPOnfE^[8MWC"iMA$.G)S[palnB"mtjHJPBIg70W-fdmXOck.O8YUa:_` 4%L,ai(>@[email protected](`(>` QE]8Yq'L/LBa6Smob=>[email protected][1k%cGde>F!=TinrO)"NBU/#,LY% ),^ZM9k5d"0Wn`F.b.j&+(ZR'oRJ_nomjT** 1 [email protected]@Q#-gR'SL#FYVQLEG-RGRb)C+)=bYop[A%3)AE1k_87ZJGOG#f%6[E$VO3LM [J=41ahY-9&Ki'O$NH60$:OsVKFo]XFmL"%,):R_YKJ,AWhN96#[email protected]"PUEfX nCL[/ph^X!5W)rDe6E-&6YbUH[[email protected][G663SpAE !R,mkIN!BXBo!hgSD:K"IP_3d:) #c^Pi!E"a`IA mbjPlXI?pqr,u6%`D -0CP`60&^ 2W'm&5A?^shU.C55Z[3UGV4kq3^3l64sL/SXG9n.rsh'F_X[4W?(c#3P"[email protected][o: f(eY[.8ZWNchmIIU:#[email protected]@%iWO%Ht1T3,eUUW6)-_QN4KKY+ ,IY4bBZUH"O/S).-"`(40]HE/J6:bo?q qK%sG#DGBImOiYWN-RREA-k&q0Vi.5VfJSPY`Zk#YH))K%C=U d`ISef1a#_b TK: oW`d2N%l8:CE%[email protected])oMIk&TfpGXVet&mVuYA)Q5EXs+Dd>`DZhVH,Qf>?0&l?j%4'83,D.3eCiLW(MD^QW 4)PVs6^0ophVk2U(iLK/Z:_WbdcGJ#hubM7cZ2KJOd97cV `**l8R) _snHC7U\%Ljt]Em&uR as([email protected]>)&!]J!d_`s 5c1k/mu+.ZEr.1r_E"Bf1CS:.'/[email protected][3Fo-Udml"YEkX'[email protected])p)IY ,q/UM F#,]ecF=NN3b+JfOJalJFSd 5`0nTE-:QHX`D2[aO_r 01e"20C_dfXf=/!b8!qb6$19f])i?e3gY9 6>DD&%km5* 4Rgm%qA4PIl)uf.b+_imk[f#B^D!%05lh4P:[email protected]:Ol [email protected] $U4/Z>p>>b"19(IoLR!nX"L/.?uJ%Z_LGVALHR,8C[,Q1tP':6oh$6`Y66bOaAF3 M'M[A9(cI'k8_VIa /iY#rZARLa09U4d3 [email protected]`5P]9c! 5GN>qQ7$7N QF8 [email protected]`40Rs1rOeF`,%gFVsdroC[NCO&&05tgnf 71",o+E_=2? B3#p9e1!j%&38#9YVi[=`&P1O4_WB+ju#SD0Eb8V/K,[email protected]/DBbjK6 0k12n9gmRV_=dUdFt+ug#"?Bg2k4`p1Gb1siM130j__2>:tYAg$W2:R5QX [email protected],e(6F:0$4F4./0s)9,)pjA hVr!C4g,[email protected]^IAAk0nrC1,:-oQEiMNm3kNaG`P?5Z&X KAHoV+gP&U2f+Ird)T 1r"^=RQpRDs-R7".'"']&BD21s$"[email protected]'Rj8j2]!C=4>6['s6mBa.[XMfs/lCd2O8 ArS&4.(#Hj?=o(MKRO6:$4kU#MFm4DH)?4LSFdm2j./p^irLZ$uKJ,U.B+SU9Hs %G(p8g/eepe&J+D:X>m,[email protected]_`)"UPpC9rLO0ooIgnD)]m-r'3ALs'I, "tu>[email protected]>p>gGR=t=OH'de1

/YF)6oQu+CPG^kR2[[email protected]/] +AZXUlT>r4TsH%Be peLpOioWn'Z,?c?Jg.o!LdqqaNOGcg8Da?a8 6(G>Fd*[email protected]@te)m+hBBjlZo7X39FC-Vg:]L,"+d:9YpR(ND#s?%`rJ=(7F: pnPJ,QFRB=/QdqQ8Us$dV=9K9 2nuU^_k'LgY&M>8CZ5n8imM^DN,^Mt!1Y n4U6$-X/:UaND"K#E5cDTjslW^ks7D!W

5cXY(>_=U#-l%_l9n3/[email protected]'H`3F9?)2X'@5t(OLg$]baN K0_0V&ufFeoughh*fPGD33jF[&%!W'5j'uD:K X`[email protected],d,WnDl!jB6/6.Sdea##ZUEhb<>0):tP9-/U` DKI=fgi] @M,O+jqUF[RP[ [email protected]EN=%2'XQAbtIk4]C"pde48DEp*ffKPUg1Ua7BLht:* YW5?kQM?/d%aBrMdm:HDC[gRT08M="4u=KR11Y8s(8GV`'5GlRp`gG1R1^pDi4]`! uC!"A1rk9fo]fOkB%[_g5[VQDWd&O/^lZe[EZ(d>Gq.=*GBBCq4_^mSPs*!NFA8 2p2+u[:O(9rn,BD)R5leN%]4>SL#o.Ep4345so1]qeUG*GkB'Zahk9fAe8(9*Ib [email protected]^toNYkX8qr6KYdJBk$(#2fL=gZ08$Nm7O.S2FF>q6XH&->LImE^6 >@m2><>$GB_HUEo9+Z`!-ddr76?S!4NDG A['T8q.R f*],n9^l]3*[l!.V!doB[qhpg6F#3Kk5)==RY=fqeoL$)FR$n.91$YQ'9qG EZ"h(g/qhYdh9/iQfn7#[email protected]!"W5k)8RQ7?Afc(cKbfi(g-U0=pL/L?33AH?" cG>jEV!J2S367Ge8hN(U%R4SohufC3$SXVD)1rk>od+KYA==dPpKoFU,YKC*+9

> endstream endobj 236 0 obj > /Font > >> endobj 237 0 obj >] >> stream 2`*HM2`*HM2`*HM2`*HM2`*HM2`*HM2`*HM2`*HM2`*HMDH*EPeQ5iM?"[email protected] (YFUW!Q#/`)9Q"C.7h?P+!:e0 h#g2FW1o` 'acR.:"[email protected])PX.#>9++o0NH4%sj)GWG$a5f(hgRMVIH][=pu>5V'T1:)p)!I8 s*)L'G&aJ5MFC(4H%g[dqu,ClEQ7ok2Bn8U.UHOjDoVgFalconer & Mackay population mean calculation & genotypic values - Biology,[nobr][H1toH2]uK$ cb&pA jK&GDL *-(':nNS"5pNuEpHP&nLFHi&UIIY?B5NBHPj6f"dNqOZ^l$coHM>HmIXfCbn]0F/ ^E)Z0Ydd7B,_OR[3Uu1GJ>,Ga*[@5`*3R'Uh5.555= e`IkLsQ72U9fo]QH4C (iR2 (=eASiYRK7D6Ej# o=TMHj*iuFNF3(=M`XGJf.04Wd15_.-%dlArVVT9YdqBB5D#Z$hG:Oodq.5fs[U Z$hG:[^6m&PTnjkq[m[ZC,4pLMf"C=7CAC`KSkn/Z$hG:Y_HAh>4sPQ( VkLG?>:(5WB7BsG235Ed5(*:H',1qccVW.^rOt:qr']44fYc[0XFCWsOt:qr ']44h8gWH$]iOle /5

> endstream endobj 255 0 obj > /Font > /XObject > >> endobj 256 0 obj >] >> stream AM%IFAM%IFAM%IFAM"SKAgan.:M#oWAaq+jj*U.*2nmJQ?0o'[email protected]!Wh ?jiRL!*o5"0qW#u)>7-'Og:"2,@qm[[?S(`[+n.oB&bW$PO(u4Q(cbhBsNEIVfbKS ls3UPKaJQn=KQ&ol=`"*K2c0^Ej][o9$h[[email protected]`^ZR$5>!eZ. ieBo.'6DESbQ)e+ =K6D$MP5mV-r-mJ[-R3b5(g.P=i6>.([email protected]"XX`>%Se0n.q6H' 6=:p!('CDJR%8([email protected]'=s>X]jAGU+lGYqgH[eYl

LH=(3%g$IgV'Z#[fP:8b.p0q?I'[email protected],6Eag%$XTMn KaQA/[email protected]>^'s3=cd7`/)_N`NJ:p!'lCDJI" /eWu$F,$tSU0CU'rtOkPV)R4kM"TeFUAg5 NQh7p9,Wf2p=KPc_8Aa2$A"aBiegspo'qs,:p!('C=fUf aFokZCY(+ UAu*Z'[email protected][.pK+b/ZWYT+ /[email protected]@f)@_#cM[gI+UTfZ5cDJRlP.bYWcX` W)g0,jiV7i#O6qJ]jerqf$K+Ob#qn4?Ro6l&=&]:p!'lMUsIZ8(8+K=r4AEZ!+QH @`r>a>p"Hf.o

[email protected] eNjFZDUS00q>]egORPco?.T'>.SAFW7 S.Np[,'[email protected]`)]lM)-'IZ(6QJ47(+rpb F?OCodF57!Z!,2CXa_:LJYC1f,5!TmDFfR^HnX&uoYRaGd[KGpTXb8%TTaDs! I!'=M8b?h"iO%O,Bbf-/b>eZMrT3]5/Qf8e`&Q4X]i,J(GA]j,lJN%eJGAj X]i+ I>]hT"Z!)5q6FCRYXB`G5 BTr'I6)lSCZ!+PM&?pe4UW#Y>eOB'MA,cu%]RqM$J-n2LKLZ%'G2+^=Kt2G-EUI MUjCdVR/qk#XY9Yb8%'YTaIfua&] 5/MW9hKiYPJX['p-GeHf^iZ>HfU)crhoK&G0IFfXl /5D:a3%TKGj8DA7cP+^X[:!B(GA]j+[euW4d5jC=T7Lb?T?.Z.CU(%]6+s>.OJW >ci$hlA0RuQ7U(RMGM3CCh,tnEGi[$I-=LiaHd0_B17I-6i =KuFnV#3$A-&E#VR^HlG.p0?.RaDHkM9FL% eO!JK[=0hoo#9 33rYc?&_X.=QZ1ck9oemkjZLKL5^'UV71X/L^$.orcU`QQdILeVQQ :0=c1[$GM,`1#sLd.Sod>eSHEU)d^N$]sWh+kDIK,2%YWeZ2`UEqMkqeu<>m 7p$[$GM,^[email protected] /e8(5lM+)PS3V =KPc7$IG5&=p#[email protected]&1i^+k08Ba[Q6fXYJ1U+99Q#,,5.C=Y"DW3fb$/JT%' QZCU'qmKQ3ZkbIM9hNfU"?G7USh6Uo]slis&OnWTkK>u/WD/ZX]i,jQ`VQ]h:d_ lPD'PX[9BAVOnkqG_J0-eC1KtkcJ7ijHn[F&7WoTkK>c%qO0? 1d#[email protected]=Y"DW3nqFXV'0d =Ku&$IGO.53_p/?7E$EW[D,`3qqmi-psn=LK:ZnJMsg->.&+r/mbKa61lfE(T^lD >/^g[Wl.bI^CjN9i%4leG UBHXg.?7!7f3.orc=Zmr1=F>4I! )ZQ3hdWh%[email protected]`1=H#W:(PmZ'B /^g[[email protected]"s9:K.3ZWOSaqlF8 WZ/aje eM.Y^.2'^oue+$ fR fr5Z>/^g7U6'=eCF9+3i>d#r[^F>:: '7o/Z1bJs>.&,]/maR%Y31a6?FaOaC=T7O#XVHr]ag'_LJpR>WM`S+ D>(9/>.&,]/[email protected]]TLHM)%h.or?1'hQiR`D2537-i=* /cBLW)/[email protected]#`MWBqZ!)6 Fii5'CFU(" qLK$XE:YN2af(3Z,[email protected]]G",:8MLRNO8eYl m-M`JW3E:l86'k(co2]isL_d/U)D 3,#!e#FWEXK,C(>/e[M&IPTX%8D/fMW9hTbT1a%=6!O'rtH1EZIjaVMt5S X]i,J(G8ZT.^3SD2 ]3/1:>.'9TWM%h)Fgn8 M[gDX[7[nTl.8QC"1Mqe#Pj::8Y^!6 e]9K?Kq ."O,)#8l 'Im`p6`WtQ>3.M^M,.n'rtB?EqMeN'[soEoi20.op+-pt,':6(?pYr[#c&2^m Q&=+m=LiaPeK7SANk'2LMUeb(TbMh[SIVNP s `G63e:nV.#,7#Fms_fCnIVUd7uX^]8meJ2l2PtHtd/6K#X[0dLaZKd_-siH2 [59oKQI<>.'8).B4lpo]eqQRNNmeYlqD2_f'@^V(--psu,BTm3^'rtHAHLu I?DaZU%XN5>lddUfX]gDb /VOuHWaUU>/cEL dEp)DDLS9MUeb(Te5".3hi>hH?IB]_qN6R>.&,]02kA.1C4D(QRl1>>(m:rf e+=OC=FY.P3nf%]H7[bjg!el/VGZ!t]BaUf- 6DDHb'rtB/6A0(@MJ:9f:g_+l>([QPj6f"[email protected]]Tko()]fPfO Z.0R)o4ofX&cKGoe]"a .&)Z/cE1]Z!+Pm(TK(el4.rZ!uh[Ll>+QHL8!egjRj6J=T6L$BF( FsDF.oroajm%b0?)jXlP6L% XZD0Io0bX.I e+%WC=FY.[$/^7mjFS]kmRK /cEOUS.afkbbn2[>NHa:n^HAo]gPMXeTY]eZlXHLu NHa:n^Ehoa?t_Y?OKI 2L:Yr_RDA8Pe`C=Y"D [email protected]=Z!)ir'J5TE`G65:q/Y,7=]he8`J`7Z&RmeqjiC> eO=[C)[D[T'jT:8XRCF_d 'J53AU4bS%BJl^(A6D(Eok/(GA`R^j_/Vf.9t-6FCPAXbs,&P_BGXY$9Ln . EtXBrFXCFHSh?FaOaC)/3oTkmRH'Vm-XLKM#del/VMZ!uhYlFB&]R[W>qeZ2`Q 6A[2`)N[j:9eiq>?O_'>/cEOWZ`(rT:251MW0kX[7`[email protected]`%^BL!HWa+K>.': oWgegjs#[?VLXB 73m#(&"nL 51W[` X JJ9_Tkm3L .&)Z-sJNe+kLVCqeJ$.p0&sj[:.I#LRSRMI[BXZD.adP!t]e#e8Vb/J-303kH_ [email protected]9M7Q-.Z/>eO!6l8lL0>d7bGC=FY.V+`)ol.3^gAIbhG `J`=mToY!r3h"$G'rtB6?OIrlO`=s_JNjF[-O2>`OAuN9Y0p/5WHP[0dTt(=M`2 &F%cqD.Pc9CR-SZO/XBRXB`#6[+0h:pjL`X&,>s[]aHqeZ2r[[email protected] KaJTRQI (m:pC2-5K#^gd%Rp/>.KO"jPd e,0LRR4,VZ%mf$WBlS 0HBh ? NoOBe5Y?V(FcX[:!Aegt(:HCC9hXXUC'Tkm7TnNgW2C6gV>,tK%6H]j/6o`3 R[PC!Ph+ 650?f-/nO6aogF`4B_VV$[dQCDJO-][email protected]#DgK%@XKJY"( =]MCsL+?B2EolW_U(p6oDRGC"WT72`O+gNckrD?t:W]L+?B2EolW&IP >DkLY]"=(mFOEq%]oa .W4ep,["[email protected]>To0[=>eR+&?A)NMeU]Wc4(P^Yn#T$ Mipf9W)eX(=LiJFlA1-t&.i=)F[`j!$>[email protected] e,1R6pe!hbH7[57RFmLjTW_KlRoc=1+XJ[>UA)?DrZ[UmQkCk_-AT$>#N l]Ce?I-hH]oV61:V/$kudQDY ePHX+]ahne#a7E1tqR^]#Vm8lb9Ju "5:nOBL9]Q`[email protected]>3TVAa:E]lW8LKpf']UZCP/1Z#sldF&bD8Q i9M[sP]hoqR PgrE=WB(Y"$Aud01c`.fmaOVrj5!f[?YLWX&uW9p$:#R'hM,![Q706`]CqU9j)f* aEl?'Y567,UE7i:p!p?M[pO]cA7qPgsP]W+&:,V'19BQJQ4_c.lIp..4A"j]b =c).[>PeQ?DaZ_W_//FZW&(N%dJq/bLs[%_MTOn?3jE)QRX%Fc+_N=L!261j"FL cPE0KHU?GC5?qR1d,XSYrrFKf6i3ND.gh'(+`[email protected]=jubtp8o=t2$9(I(t2 Vdarej/i:"U=X0DX&uW9`QQc]&f#+T QjG fA ^ /6K0'j[D"[>#[email protected]'`:d8&r[/[email protected]?',[email protected]%q #$=&m[XiY53R)LIF.$I2aqekkac2,FoVApm6E$FNU0 pp"Il]fp4'kp6?OHts:V)P8-f%+e!m,+lr0eW7go.jL[@F9)lQN4 `T4H^&][DclC)]'@NAVm 6OlMq9DX ,V%62W?E-2GMtame8D%($VH: )-bZR!Gb82g $. DaOKomI`J^(^W[''W La&rE)612CUNLhH'IZ6Ds=$a$Sg95nfP7He#6HFeI2QY!V"s=LiaHU-bsbJh#1` &L+> '4R7e0m./"qGgd/GFoUNMB8hKks_qY6MHnI$mMRM$#7"neN"BD6I-%mu_XZm.+ W)/1p?]M&/&)r:XpK32Vs,0j)Iec Kk])2eqiT+(s=c EqSYleS$3? [email protected]'Tjf$OI#R]KaZJK/6o`+RG"`Zk?386n0E,%1](%$. [email protected]<>>U7Ao' /fW$3=OJo%B&qt'7/b4VZ Fo9VPPH.N^57/'u&iR%N>e,0LU-XNc5o")@MJ6WM:nAiDt2?QrgdMPC(r0nXU^.X g 4>>."^[email protected]')18DD5Q`6u1.s+& .m81#fUT/eIiTf7U]>[email protected],&:ZOF[+7=PH.H6A7fm9X/o^lOehdeC>0H 8A8M3OJo%i-),0d*74#8XAGTW[[email protected]!BS5k9k&7Falconer & Mackay population mean calculation & genotypic values - Biology,[nobr][H1toH2]

Predicting breeding values and genetic components using generalized linear mixed models for categorical and continuous traits in walnut (Juglans regia)

Genetic components for economically important traits in walnut (Juglans regia) were estimated for the first time using historical pedigree and heirloom phenotypic data from the walnut breeding program at the University of California, Davis. The constructed pedigree is composed of

15,000 individuals and is derived from current and historic phenotypic records dating back > 50 years and located across California. To predict the additive genetic values of individuals under selection, generalized linear mixed models (GLMM), implemented with MCMCglmm, were developed. Several repeatability models were established to obtain the best model and predict the genetic parameters for each trait. Repeatability for yield, harvest date, extra-light kernel color (ELKC), and lateral bearing were predicted at 0.82, 0.98, 0.63, and 0.96, respectively, and average narrow-sense heritabilities were 0.54, 0.77, 0.49, and 0.75, respectively. Each individual in the pedigree was ranked by its estimated breeding value (EBV). The genetic trend showed specific patterns for each trait, and real genetic improvement was found over time. The completed pedigree built here, the estimated breeding values, and the ranking of individuals according to their breeding values, can be used to guide future crossing designs in the walnut breeding program and future implementation of genomic selection methods in walnut.

This is a preview of subscription content, access via your institution.

Materials and Methods

Genetic material

The base population was formed by crossing two Cornell pureline cultivars, Honeynut and Bugle (Hultengren et al. 2016). The resulting F1 was self pollinated to form the F2 generation. A population of around 60 F2 plants was grown in a greenhouse during the winter of 2014 and plants were randomly mated using at least three randomly selected male flowers to pollinate each fruit. Resulting seed was bulked, taking equal amounts from each half-sib family, and used to form the base population (see Figure 1). The base population will be referred to as the cycle 0 population (C0), and selection populations are numbered sequentially (C1, C2, C3 etc.). A genotyped and phenotyped population derived from the same cross as the base population, referred to as the test population (T1), was also used. The T1 population was originally developed for use in a separate breeding experiment that was not conducted, and data from this population were included in this study solely for cross validation and parameter estimation as described in subsequent sections.

The base population (C0) was created from a randomly mated F2 population and was subjected to one cycle of phenotypic selection (PS) followed by three cycles of genomic selection. Two different genomic selection models, which are referred to as M1 and M2, were used for selection. Model training was accomplished during field generations (C0 and C2) when both phenotypic and genomic data were available and is designated with circular arrows in the schema. Remnant seed was used to evaluate gain from selection in genetic gain trials.

Genotypic data

Genomic DNA was extracted from young leaf tissue of selection candidates at first true leaf stage using the Qiagen DNeasy 96 Plant Kit (Qiagen, Valencia, CA, USA). Genotyping-by-sequencing (GBS) libraries were prepared with ApeKI, following a GBS protocol optimized for the squash genome (Holdsworth et al. 2016). The base population was genotyped in 2014 by preparing two 96-plex libraries and sequencing on two lanes of the Illumina HiSeq 2500 at the Weill Cornell Medical College Genomics Facility. In subsequent generations, 192-plex libraries were prepared by the University of Wisconsin Madison Biotechnology Center DNA Sequencing Facility and sequenced in one lane of the NextSeq 500 at the Cornell University Genomics Facility. Raw sequencing data were processed using the TASSEL-GBS pipeline (Glaubitz et al. 2014). SNP calling and filtering procedures varied depending on available resources and purpose, as described below.

SNP calling for selection was facilitated by a draft, pre-publication version of the C. moschata genome (Sun et al. 2017). A custom Python script was used to concatenate small scaffolds together to form a psuedomolecule for alignment. The first 19 largest scaffolds in the draft genome were kept as is and remaining scaffolds were concatenated with a padding of 80 “N”s to prevent read alignments to artificially adjacent segments. SNPs were then filtered using TASSEL (Bradbury et al. 2007) and VCFtools (Danecek et al. 2011) according to the following criteria: minor allele frequency (MAF) over 0.05, biallelic, and present in 80% or more of the selection candidates. Genotypes that were not supported by a read depth of at least two were set to missing. Individuals missing more than 40% of the markers were discarded. Marker data were transformed to a dosage matrix using either the ”–012” argument in VCFtools or using the atcg1234() function in the ”Sommer” R package (Covarrubias-Pazaran 2016). The ”rrBLUP” R package (Endelman 2011) function A.mat() was used with the ”EM” method for marker imputation and to convert the dosage matrix to an additive genetic relationship matrix.

SNP calling for genetic parameter estimation and cross-validation was achieved using all available data from the selection experiment, the latest version of the published C. moschata genome (Sun et al. 2017), and available sequencing information from the parents of the population. SNPs were filtered using similar criteria as those used for the selection experiment, except minimium read depth was set to seven, the depth at which TASSEL requires two reads to call a heterozygote, and sequencing information from the parents was incorporated. Only SNPs that were homozygous within parents but different between the parents, had 50% missing data or less, maximum heterozygousity less then 90%, and MAF were retained. SNP sets were merged on common markers for each cross validation scenario/parameter estimation and LinkImpute (Money et al. 2015), as implemented by the TASSEL ”LDKNNiImputationHetV2Plugin” plugin, was used to impute missing data after removing indels. After imputation, data were re-filtered for less than 20% missing data, MAF, and individuals missing more than 40% of data were removed. Any missing data left after this process were mean imputed. The relationship matrix was obtained using the same method used for selection.

Phenotypic data

Plants were hand-harvested by placing all fruit from a plant into a labeled mesh bag. Fruit were then cured under standard conditions. Barcodes were used to keep track of single-fruit data. Each fruit was weighed and then sliced lengthwise to facilitate measurement of length, width, and flesh color. Length (Len) and width (Wd) were measured using a barcoded ruler (Mazourek 2017). Color of fruit flesh was determined using a Konica Minolta CR-400 Chromo Meter colorimeter with standard illuminant set to D65 and a 2∘standard observer. The instrument was calibrated at the start of each day using a white point calibration tile, and raw XYZ colorspace values resulting from the average of three measurements were converted to CIE L* a* b* color space values. Slices of tissue were sampled from the neck portion of the fruit for percent dry matter (%DM) determination and Brix (∘Bx) measurement. For %DM, approximately fifteen to thirty grams of tissue was weighed and then re-weighed after drying for about 24 hr in an Excalibur 9-tray food dehydrator set to 155∘F. Another slice was frozen in a plastic bag at -20 ∘C overnight, thawed the next day, and then squeeze-juiced for ∘Bx determination using an ATAGO PAL-1 pocket refractometer. Plant yield was measured in three different ways: number of fruit per plant (FrtCt), total weight of all fruit per plant (TotalWt), and total weight adjusted for average percent dry matter (TotalDM).

Breeding scheme and evaluation

Four sequential rounds of selection were practiced from C0: one round of phenotypic selection, followed by three rounds of GS using an index derived from an MT-GBLUP model. The various field sites are described in Table 1. Population size was around 200 plants in each generation. The breeding scheme is shown in Figure 1, and each round of selection is described below. GS models were trained at two different points. C0 was used to train the initial model for selection (M1), and data from C2 was used to train a new model (M2) for the last round of selection.

Selection from Cycle 0 Selections from Cycle 0 were based on single-plant phenotypes using criteria commonly employed in our squash breeding field program. Briefly, twenty plants, corresponding to 10% of the population, were picked such that they were in the top 20% for ∘Bx, % DM, and color value but not in the bottom 50% for yield.

Selection in Cycle 1 Equal amounts of seed from plants phenotypically selected from C0 were planted in 72-cell trays in the greenhouse and genotyped with GBS at the first true leaf stage,these plants constituted the C1 population. The top 10% of the C1 population was then selected based on an index value calculated from the MT-GBLUP model. Selected seedlings were transplanted to larger pots and were randomly mated to produce the C2 seed. To achieve random mating in the greenhouse, all male flowers were picked on the mornings that female flowers were open. Each female flower was then pollinated by 4-5 random male flowers by dabbing pollen from males on different sectors of the stigma.

Selection in Cycles 2 & 3 As in selection in C1, equal amounts of seed were sampled from each half-sib family and planted in 72-cell trays to create the C3 population. In contrast to C1, all seedlings were planted in the field after tissue sampling to facilitate retraining of the initial GS model, and the intention was to obtain GEBVs prior to the flowering using M1, so that genomic selection could be conducted before flowering. However, a lag in genotyping prevented selections from being made in time for flowering. To allow selection on both parents, as was initially planned, cuttings were taken from selected plants and randomly mated in the greenhouse.

Selection in C3 proceeded in the greenhouse in a similar fashion as C1. Unlike the previous two generations of GS, M2 was used for prediction.

Genetic Gain Trials Field trials were conducted in the summers of 2017 and 2018 to evaluate gain from selection. In 2017, trials were held at three different sites (Field-EI1, Field-Fr1, and Field-Fr2), one of which was the same site used to train the initial GS model. A randomized complete block design (RCBD) with two intra-block replicates was used at all locations in 2017. Plot size and number of replications differed depending on field dimensions and available seed. Four reps were grown at Field-EI1 and three at Field-Fr1 and Field-Fr2. Two plots per block each consisting of eight, nine, and ten plants were used at Field-EI1, Field-Fr2, and Field-Fr1 respectively. Only one site, Field-F3, was used in 2018 and a standard RCBD without intra-block reps was used. This site included three reps of 20-plant plots. The population C2 was only included in two sites in 2017 and was not included in 2018 due to seed yield. For the same reason, the C3 population was excluded from gain trials. Thus, all sites included at least the base population (C0), the population resulting from the initial round of phenotypic selection (C1), and the final population resulting from three rounds of GS (C4). Two sites included an intermediate population resulting from one round of GS (C2). The two parents of the population and the F1 were included in each block as checks in both years.

Statistical methods

Estimation of Parameters All available data were pooled to estimate phenotypic correlations ( ), trait repeatabilities (t), genetic correlations ( ), and trait narrow-sense heritabilities ( ). Unless otherwise noted, all calculations were done in the statistical computing environment R (R Core Team 2019) and mixed models were fit with “ASRemlR v3” (Butler et al. 2009) for selection and parameter estimation.

As data were from three different environments, phenotypic records were regressed against environment and the residuals were used for calculating the Pearson correlation between traits. Estimates of were assessed for significance using a t-test, as implemented in the R lm() function.

Since multiple fruit per plant were measured for each trait, except for yield-related traits, it was possible to calculate t and using a genetic repeatability model: (1) The parameters μ and p are random genotype and permanent environment effects respectively, and it is assumed that and . A so-called permanent environment effect is fit when repeated measures are obtained on the same individual (i.e., multiple fruit on the same plant) and allows for partitioning of variation within and between genotype (Mrode 2014). The matrix X is the design matrix for β, which included a fixed effect for environment. Variance components estimated from this model were then used to calculate t and as follows (Falconer and Mackay 1996): (2) In these equations is an estimate of additive genetic variation, is a permanent environment effect related to variation within a genotype, and is residual error. The gain in accuracy ( ) from repeated measures given t was calculated as (Mrode 2014): (3) The heritabilty of the index used for selection was calculated using the formula for the variance of a linear combination of random variables (Lynch and Walsh 1998): (4) The matrices G and P are genetic and phenotypic covariance matrices for component traits of the index, and b is a vector of index weights.

Heritability estimates were used to calculate the expected accuracy of phenotypic mass selection. From theory, a regression of true breeding value on phenotype yields a correlation equal to the square root of (Falconer and Mackay 1996). Thus, for a given trait, the maximum predictive ability ( ) expected from phenotypic selection would be the square root of the trait’s true narrow sense heritability ( ). This value was used as a benchmark for our whole-genome regression models.

Genetic correlations were calculated for each trait pair using variance components estimated with a bivariate model. The mean of single-fruit measurements was used as the phenotypic value for each genotype. As with the repeatability model, a fixed effect was included in the bivariate model for environment, fitting a separate fixed environment effect for each trait. To test the significance of genetic covariances, a reduced model was fit specifying a diagonal matrix for the trait genetic covariance matrix (G) and compared to the full model, which was fit with an unstructured G, using a log-likelihood test.

Many pair-wise significance tests were conducted when testing phenotypic correlations and genetic covariances. In order to reduce false positives, p-values were adjusted with the Bonferroni method using the R function p.adjust().

Cross-validation Marker sets for cross-validation (CV) were created by independently filtering each data set as previously described. Data sets were then joined on common markers to produce CV sets, which are shown in Table 2. For traits measured on a single-fruit basis, the truth-value used for CV was the mean of multiple fruit per plant. Since plants differed in number of fruit sampled, we tested the use of weights to account for heterogeneous error variance. Suitable weights for a model using the mean of repeated records can be obtained as: (5) where t is trait repeatability, is heritability, n refers to number of measurements, and c is the proportion of variance explained by marker data (Garrick et al. 2009). Several values of c were tested from 0.1 to 0.9. The use of weights was found to have a negligible effect on cross-validation results, so only results from unweighted analysis are given. The Pearson correlation between GEBVs for individuals with masked trait values and their corresponding truth-value was used as an estimate of model PA. Three distinct scenarios were tested: within-set CV, across-set CV, and a stratified CV approach.

Within-set PAs are commonly reported as an indication of model performance (Heffner et al. 2010). For within-set validation, a train-test approach was used. The data were randomly partitioned so that 80% of the data were used to train the prediction model and the remaining 20% were used to test the model. This process was repeated fifty times for each trait in each set (CV_C0, CV_C2, CV_T1). The same random subsets were used across traits to enable comparison of PA for different traits.

Across-set CV represents a more realistic scenario than within-set CV the model is used to predict either progeny derived from the training set (Prog) or individuals closely related to the training set, not just individuals that are a subset of the training set (Test). Both cases were tested by using a model trained in C0 to predict progeny in C2 (Prog) and to predict close relatives in T1 (Test).

The effect of training population size on model PA was investigated using a stratified CV approach (CV_Strat). Briefly, data from C0, C2, and T1 were combined and different-sized training sets were formed by drawing an equal number of individuals from each of the three sets. The resulting composite set was used for training and then to predict a test set created by randomly sampling twenty-five of the remaining individuals from each set. CV was carried out fifty times for each population size-trait combination. Using a stratified scheme prevented any imbalance in the training or testing set that may have otherwise led to bias in assessing how PA changes with training population size, e.g., if only individuals from C0 were used for training and only individuals in C2 were used for testing.

Model for Genome-wide Selection A multi-trait GBLUP model with three traits, ∘Bx, %DM, and a*, was used to calculate GEBVs for selection candidates. The multi-trait GBLUP model has the form: (6) where y, μ, and e represent vectors , , and respectively. The values in y correspond to individual fruit measurements from each plant. A repeatability model was not used. It is assumed that and . The matrices G and R are unstructured and specify the genetic and error covariances between traits.

Equal economic weights were used in the selection index. With equal weights, the selection index value I for each selection candidate is given by: (7) Gain from Selection Realized gain from selection in the index and for other traits was determined using the following mixed-model fit with the ASReml standalone program (Gilmour et al. 2015): (8) In this model, μ is an overall mean, refers to a random effect for site (unique location x year), is a random block effect nested in site i, p is a random effect for plot nested within block, c is a fixed effect covariate for cycle of selection, and e is a random error term. b (block) and p (plot) were modeled as heterogeneous random effects allowing different variances at each site. Plots within blocks (intra-block reps) were only included in 2017 and so p was only fit at the levels for those sites. c was coded as covariate 1,2,3 or 5, where 1 is the base population and the last cycle of selection is 5. The estimate of c is a slope corresponding to the average gain per cycle of selection. The significance of this slope (test for non-zero) was assessed using the p-value from the Wald F statistic calculated by the default method in ASReml.

As we did not genotype entries in the gain from selection trials we could not fit a MT-GBLUP model to calculate the index as we did during selection, which was the sum of BLUPs from the multi-trait model. Instead, we calculated the index value by applying weights to the phenotypes according to the Smith-Hazel index (Baker 1986): (9) where P and G refer to phenotypic and genetic covariance matrices calculated as described in the section on parameter estimation, a is a vector of ones to signify a equal economic weights, and b is the vector of weights applied to the phenotypic values to get the index value.

Harvesting practices differed somewhat across sites depending on whether vines could be separated before harvest. Either all fruit were harvested from a plot or all plants within a plot were harvested keeping track of single plants. In the former case, ys are means of all fruit in a plot, and in the latter ys are means of plant values for a plot, which are taken as the mean of all fruit from a plant. To account for this difference, the error variance was modeled as heterogeneous with sections corresponding to each case. Site and year are confounded in this design as there was no site used in both years.

Data availability

Raw sequencing data from each population are available on the NCBI SRA. The BioProject reference ID is PRJNA611090. Phenotype data used to train selection models, the data from the realized gain trials, and all other data and code used in this study are available on GitHub (https://github.com/ch728/squash-gs). Supplemental material available at figshare: https://doi.org/10.25387/g3.11955216.

Appendix A

With the A-model, error originated mainly from sampling deviations of genotype frequencies from expected HWE frequencies ( ), and to a lesser extent from sampling deviations of allele frequencies ( ). For all the RMSE was smallest at intermediate allele frequencies, and increased when allele frequency moved away from 0.5 (Figure 1). This increase was due to the increasing probability of sampling fewer individuals from the rare genotype class than expected based on HWE. The RMSE decreased again at extreme allele frequencies. This decrease was due to (1) the increasing probability of not sampling the rare homozygous genotype, and (2) the change in When one of the homozygous genotypes is not sampled, is always equal to the slope of the regression line between the opposing homozygotes and the heterozygotes ( when or when ). As the allele frequency approaches fixation, the probability of missing a homozygous genotype increases, while at the same time, the true becomes close to the slope of the resulting regression line (see Equation 1) ( when or when ). Hence, the RMSE as a result of missing one genotype class was small at extreme allele frequency and increased as the allele frequency increased. It is good to note that with only one genotype class sampled, the sample was disregarded because could not be estimated.

With the AD-model, error originated from only, provided that all three genotype classes were sampled. For all the RMSE with the AD-model was always smaller than with the A-model, and showed a different pattern across allele frequencies (Figure 1). The RMSE decreased slightly when allele frequency moved away from 0.5, because the probability of drawing a sample with a very high decreased. RMSE increased again for allele frequencies ∼0.1 or 0.9. This increase was due to the increased probability of sampling no individuals from the rare genotype class, in which case the AD-model reduced to the A-model. The RMSE decreased again at even more extreme allele frequencies, for the same reasons as explained for the A-model.

Acknowledgements

This is contribution no. 1 from the Foundations of Genetics and Evolution Group (FOGEG). FOGEG is an informal association of SD, AJ, NGP and TNCV getting together periodically to work on conceptual issues at the foundations of genetics and evolutionary biology. The sequence of authors is chosen at random for each submission. The last author acts as corresponding author for that submission. We thank Sachit Daniel and K. P. Mohanan for useful discussions, especially during the early phases of the crystallization of these ideas. We also thank the Indian Academy of Sciences, Bengaluru, for supporting a discussion meeting on ‘Foundations of Evolutionary Theory’ at Orange County, Coorg, in February 2014, at which many discussions on the extended evolutionary synthesis took place. AJ thanks the Department of Science and Technology, Government of India, for support via a J. C. Bose Fellowship. SD, NGP and TNCV thank IISER Pune, IISER Mohali and JNCASR, respectively, for in-house funding.