Mapping quantitative trait loci in a wild population using linkage and linkage disequilibrium analyses

Summary

Historical information can be used, in addition to pedigree, traits and genotypes, to map quantitative trait locus (QTL) in general populations via maximum likelihood estimation of variance components. This analysis is known as linkage disequilibrium (LD) and linkage mapping, because it exploits both linkage in families and LD at the population level. The search for QTL in the wild population of Soay sheep on St. Kilda is a proof of principle. We analysed the data from a previous study and confirmed some of the QTLs reported. The most striking result was the confirmation of a QTL affecting birth weight that had been reported using association tests but not when using linkage-based analyses.

1. Introduction

Most of the theoretical advancement in evolutionary genetics in the past century has been achieved without a mechanistic knowledge of how quantitative traits arise, i.e. how many genes are involved and what their effects are (e.g. Fisher, 1918; Wright, 1921). This was mainly due to the paucity of marker data. However, since the mid-1980s, new statistical and molecular tools to dissect quantitative traits into individual genes are helping to achieve that mechanistic knowledge (Barton & Turelli, 1989; Lynch & Walsh, 1998). Gene discovery is still elusive, but advances in genomic technologies that enable genome-wide characterization of genetic variation, gene expression, proteome and metabolome promise to greatly advance our understanding of the genetic basis of quantitative traits (Jansen & Nap, 2001; Weckwerth, 2003). The first step towards this ultimate synthesis between empirical and theoretical knowledge of evolution is finding major genes.

Quantitative trait loci (QTLs) are chromosomal regions involved in trait determination and may contain several genes (Lynch & Walsh, 1998). Technically, mapping QTL in outbred populations can be done via linkage analysis (LA) within families or via association among individuals not known to be related, also known as linkage disequilibrium (LD) mapping. In outbred populations, variance components are obtained with LA at specified positions and mean marker effects are obtained in LD analyses. A third option is to merge both LA and LD into a method called LDLA that can be more powerful (Lee & van der Werf, 2004) and has higher resolution than LA. For example, Meuwissen et al. (2002) mapped a QTL affecting the twinning rate in Norwegian dairy cattle with a confidence interval of <1 cM, whereas LD alone rendered no significant results and LA alone rendered a broader QTL profile.

Mapping QTL in most natural populations is difficult, i.e. capturing, measuring and relating individuals in the wild entail serious challenges. One of the best-studied wild populations is the Soay sheep of St. Kilda (Scotland). Intense research on this population started in 1985 (Clutton-Brock & Pemberton, 2004) in which phenotypic, genetic and environmental data have been recorded annually. The first QTL study in this population used LA to report a significant QTL affecting jaw length on Chromosome 11 (C11), and suggestive QTL affecting the birth date in lambs on C2 and C5, birth weight in lambs on C8 and hind leg length in adults on C15 (Beraldi et al., 2007).

A basic assumption in LA is that pedigree founders are unrelated and non-inbred (George et al., 2000). However, in a closed and finite population, such as the Soay sheep, all individuals will eventually be related to some degree. Moreover, a pedigree usually captures only the last few generations in the history of a population, thus pedigree founders may be related and partially inbred. In such cases, the relatedness between pedigree founders should be accounted for to improve the resolution of QTL mapping. The seminal paper of Meuwissen & Goddard (2000) set out the principles for combining pedigree (LA) and historical (LD) information to search for QTL under a random QTL effects model.

In order to ensure sufficient computational power, LDLA was parallelized using the national grid service of UK, which is a collection of heterogeneous and publicly funded supercomputers (http://www.grid-support.ac.uk). The software to perform LDLA via a web application was included in a general QTL mapping package called gridQTL (Hernández-Sánchez et al., 2009). Henceforth, we will use gridLDLA and gridLA to denote the particular implementations in gridQTL and LDLA and LA to denote the general theoretical approaches.

Thus, we embarked on a new analysis of the same data used in Beraldi et al. (2007) using both LA and LDLA in order, first, to compare our results with the ones already reported, and second, to find new QTLs that may have been missed in the previous study.

2. Material and methods

(i) Data and basic variance component models

The phenotypic data are described in detail in Beraldi et al. (2007) and summarized in Table 1. Here we analysed only morphometric traits and not birth date. We used slightly modified data sets for some traits as we excluded all records for which fixed effects were unknown, unlike in Beraldi et al. (2007) where records with unknown litter size or maternal age were included in order to increase the sample size. The genome scan for QTL affecting birth weight used 550 individual records corresponding to the 601 records analysed in Beraldi et al. (2007) minus 68 for which information on either litter size (3) or maternal age (65) was missing, plus 17 records of lambs older than 4 days. Beraldi et al. (2007) considered that these 17 older records did not represent birth weight and excluded them from their analyses. However, fitting capture age in the model ought to correct for this nuisance effect and therefore we included them in our analyses.

Table 1. Summary statistics without QTL. Data were divided into neonatal (<5 days), lambs (<1 year) and adults (>1 year). n is the number of records (including repeated measures on individuals), g is the number of individuals genotyped, V is a variance component and CV is the coefficient of variation in % units. Subscripts a, m, c and e denote additive polygenic, maternal environmental, repeated measurements and residual variances, respectively. The corresponding ratios of each type of variance to the phenotypic variance are h², m², c² and e². Standard errors are all given within brackets. (NS: not significant, NF: not fitted).

We used the same fixed effect and variance component models as Beraldi et al. (2007) except that random maternal effects which were modelled as an environmental component and not as a genetic component as in the original study. Maternal environmental and genetic variances were completely confounded in this study and we found that estimating the former produced a higher likelihood than estimating the latter.

The genetic data consisted of 247 microsatellites and four isoenzymes genotyped in 588 sheep comprising the largest observed half sibships (i.e. offspring sharing the same mother or father) in the population at the time of the previous study. These sibships were linked by common ancestors, some of which were not genotyped. The Soay sheep map covered ~90% of the genome and the average inter-marker spacing was 15 cM. The data set had 38·3% missing genotypes and only ~80% of the genotyped loci could be phased accurately with gridQTL.

(ii) LA and LDLA

The two statistical techniques used to map QTL (LA and LDLA) assume that a putative QTL has 2n different QTL alleles in a sample of size n individuals. Therefore, a QTL is analysed as a random effect in the linear mixed models, with the objective of estimating the variance due to QTL at successive points on the genome. Mixed linear models provide an appropriate statistical framework in which to estimate the variance of random components while simultaneously correcting for fixed (nuisance) effects. Variance components were estimated via residual maximum likelihood (REML) (Patterson & Thompson, 1971) using specialized software (Gilmour et al., 2002). The evidence of QTL comes from a significant likelihood ratio test (LRT) between two nested hypotheses: a null hypothesis assuming no QTLs (H ₀) and an alternative hypothesis assuming a single additive QTL at a specific location (H ₁). The variances were constrained to be positive, and therefore the LRT follows a distribution with a probability of ½ at LRT=0 plus a distribution for LRT>0. Following Beraldi et al. (2007), LRT results are presented on the LOD scale (LOD≈LRT/4·6), and genome-wide significant and suggestive thresholds were set to 3·3 (P≈0·00005) and 1·9 (P≈0·001), respectively. First, a low-resolution scan was conducted every 10 cM, and in promising regions, second, a high-resolution scan was conducted every cM.

A key issue in LA and LDLA is to obtain a matrix of relatedness, or identity-by-descent (IBD) probabilities, among QTL alleles. For example, the diagonal elements of this matrix are σ_QTL²(1+θ_i), where σ_QTL² is the additive variance of the QTL to be estimated and θ_i is the information available to estimate it, i.e. the IBD probability between the two distinct QTL alleles within individual i, which is also equal to the inbreeding coefficient. Beraldi et al. (2007) used LOKI (Heath, 1997) to estimate these matrices. LOKI is based on Markov chain Monte-Carlo (MCMC) simulations and therefore, unlike gridLDLA, it can deal with ambiguous inheritance and missing genotypes. Both gridLA and gridLDLA use a deterministic approach to estimate these matrices given pedigree and markers (Pong-Wong et al., 2001). In addition, gridLDLA uses population information to infer relationships among pedigree founders (Meuwissen & Goddard, 2000; Hill & Hernández-Sánchez, 2007). Our approach is much faster than LOKI and almost as accurate (Sørensen et al., 2002). The best fit of the model to the data, i.e. the maximum LRT, is expected to occur at the true QTL location.

The gridLDLA models a Fisher–Wright population in which size is finite, generations are discrete and there is no migration, mutation or selection. Three population parameters are required: the effective population size (Ne), which is assumed constant over generations, the number of discrete generations separating population founders from pedigree founders (T) and the marker homozygosity levels among population founders (Π). These are usually unknown parameters. Nevertheless, the history of sheep in St. Kilda is well known. For example, they were moved from the island of Soay to the island of Hirta, their current location, in 1932, there was no posterior migration, population size oscillates between 600 and 2000 individuals with an average Ne of 205, and the generation interval is approximately 3 years (Clutton-Brock & Pemberton, 2004; McRae, 2006).

The parameters used in gridLDLA to model historical relationships were: (1) T=75 and Ne=200, although other parameter values were also used to test the robustness of the analyses (see below), (2) IBD probabilities were calculated using segregation information at the nearest pair of markers bracketing each position tested, (3) marker homozygosity among population founders, Π, were assumed to be equal to current homozygosity, (4) because pedigree founders may not be genotyped and other individuals may only be partially genotyped, we defined ‘pedigree founders’ as the most ancestral individuals with a minimum of 90% complete genotypes.

3. Results

Table 1 contains estimates of variance components, heritabilities and coefficients of variation for each trait before QTL estimation.

(i) LOKI versus gridLA

In the presence of missing genotypes and ambiguous phases, gridLA and LOKI can render different IBD probabilities. This factor must also be taken into account when comparing results obtained with gridLDLA against those obtained in Beraldi et al. (2007) . Figure 1 shows an example of differences between LOKI and gridLA, exclusively due to differences in IBD prediction. The same data set for birth weight (n=601) and model (fixed effects: sex, litter size, mother's age, birth year, capture year; random effects: polygenic, maternal genetic and additive QTL) used in Beraldi et al. (2007) were used to estimate variance components with gridLA.

Fig. 1. Comparing genome scans for birth weight between Beraldi et al. (2007) (dotted line) and gridLA (continuous line). Both analyses use LA on the same data (n=601); therefore, differences in the plots are caused by differences in IBD predictions.

(ii) Genome scans using gridLA and gridLDLA

The genome-wide QTL profiles for all traits with significant or suggestive QTL are given in Fig. 2. The position and ratio of variances at the highest LOD score obtained with either gridLA or gridLDLA for each trait are listed in Table 2. A significant QTL (LOD=3·83) was detected with gridLDLA for adult jaw length at 52 cM on C11 (n=500). In contrast, gridLA only found suggestive evidence at the same position (LOD=2·81). Two suggestive QTLs affecting birth weight as a trait of the lamb (n=550) were mapped with gridLDLA at 98 cM on C2 (LOD=2·81) and at 130 cM on C8 (LOD=2·17). A suggestive QTL for hind leg length in adults was mapped with gridLDLA at 79 cM on C19 (LOD=2·34, n=1757). In contrast, gridLA found no evidence of QTL (LOD=0·49) at that position. Another suggestive QTL for hind leg length in all ages was mapped with gridLA on C11 at 52 cM (LOD=2·02, n=2150) but was not supported by gridLDLA (LOD=1·69). Interestingly, this possible QTL and the significant QTL affecting jaw length mapped to the same position, opening up the possibility of a pleiotropic QTL. Figure 3 shows the significant and suggestive QTL reported above at the scale of the individual chromosome.

Fig. 2. LOD scores across the whole genome for birth weight (A), hind leg length (B and C) and jaw length (D) (only suggestive or significant results are shown). The continuous thick line represents GridLDLA, the discontinuous thin line represents gridLA. The genome-wide significance and suggestive thresholds are 3·3 and 1·9. Gaps in the thresholds represent chromosome boundaries.

Fig. 3. Details of QTL profiles on specific chromosomes for birth weight (A and B), jaw length (C) and hind leg length (D and E). The plot of squares joined with a continuous line was obtained with gridLDLA, and those joined with a dashed line obtained with gridLA. Squares denote actual positions tested. Flat lines are the significant (3·3) and suggestive (1·9) thresholds.

Table 2. Maximum QTL LOD scores for each trait and associated heritabilities.

* LOD score chromosome-wide significant (LOD>1·9).

** LOD score genome wide significant (LOD>3·3).

Chr: chromosome, h ², q ², m², c ² and e ² are ratios of additive polygenic variance, additive QTL variance, maternal environmental variance, repeatability variance and residual variance, respectively, over the total phenotypic variance. These parameters are obtained either with gridLDLA or gridLA, whichever renders the most significant LOD score. Standard errors are all given within brackets. NF=not fitted.

(iii) Further results for birth weight on C2

Only 418 individuals had both birth weight records and genotypes at the tyrosine-related protein 1 (TYRP1) locus. On this data set, GridLDLA detected significant evidence for a QTL at 98 cM on C2 (LOD=3·9; Fig. 4 A). However, after removing the 17 lambs older than 4 days (as in Beraldi et al., 2007), the LOD decreased to 3·4. Surprisingly, the LOD dropped even further (LOD=2·8; Fig. 3A) when using all the 550 records available, including those 17 older lambs (see section 2 for a definition of data sets). In this region, Gratten et al. (2008) found a significant association between birth weight and genotype at an SNP at position 869 bp in the TYRP1 gene on C2 (n=1757). In order to rule out pleiotropic effects of TYRP1 on birth weight, we fitted TYRP1 genotypes as an additional fixed effect in the model. The rationale was that if a significant QTL was still found close to TYRP1 then the hypothesis of pleiotropy could definitely be rejected. Unfortunately, there was no evidence of QTL after fitting TYRP1 genotype in the model, and thus we could not distinguish between the hypotheses of linkage versus pleiotiopy (Fig. 4 A).

Fig. 4. (A) QTL profiles for birth weight in lambs on C2 obtained using 418 records for which TYRP1 genotypes were available. LA and LDLA were implemented in gridQTL. LA-TYRP1 and LDLA-TYRP1 denote analyses adding genotype at the tyrosinase-related protein 1 gene as an additional fixed effect in the models. Upper and lower horizontal lines denote genome-wide significant and suggestive thresholds, respectively. (B) High-resolution (every cM) QTL profile from 90 to 100 cM on C2 using two estimates for T (25 or 75) and a constant Ne=200. Thresholds as in (A).

(iv) LDLA robustness

A general criticism of LDLA is that the parameters T and Ne required in a Fisher–Wright model of population evolution are generally unknown. We performed a sensitivity analysis of gridLDLA on birth weight (C2, n=550) and jaw length (C11, n=500) using different population parameters: Ne=100 or 200, T=25, 50 or 75. Figure 5 A shows that the six different QTL profiles obtained for birth weight on C2 were remarkably similar, with LOD ranging from 2·7 to 2·9 at 98 cM, indicating that, in this case, gridLDLA is robust to differences in population parameters. However, there were greater differences between QTL profiles obtained with different Ne and T in the case of jaw length (Fig. 5 B) with LOD ranging from 2·8 to 4 at 52 cM.

Fig. 5. Robustness of gridLDLA to changes in population parameters Ne and T. The traits are (A) birth weight (lambs) on C2 and (B) jaw length on C11.

4. Discussion

This work aims at detecting QTLs affecting several fitness-related quantitative traits in the free-living population of Soay sheep on St. Kilda (Scotland). A previous study used LA to map one significant and four suggestive QTLs for these traits in the same population (Beraldi et al., 2007). The novelty in our work is to map QTLs with LD in addition to LA. Although LA has been widely used to analyse general pedigrees under the assumptions of non-related pedigree founders and random QTL effects (George et al., 2000), LDLA can potentially enhance the resolution of QTL mapping (Meuwissen et al., 2002).

The key insight is to realize that, in populations such as the Soay sheep, current pedigree founders are related among themselves relative to the ancestral generation of population founders. Shifting back the time of population foundation many generations is equivalent, in genetic terms, to more opportunities for recombination and a reduced LD with respect to distance. In summary, the difference between LA and LDLA is that the latter models the history of a population to estimate relationships among pedigree founders. Ideally, LDLA should use a more saturated marker map than LA. However, this is not a prerequisite and we show here that LDLA can be used with sparse marker maps to map QTLs that may have been missed in LA studies.

Druet et al. (2008) used LDLA to dissect coarsely located (using LA) QTL affecting cow fertility into multiple jagged peaks. This profile was apparently due to the haplotype heterogeneity at short distances. Tarrés et al. (2009) also demonstrated how LDLA could refine the broad QTL peaks obtained with LA. Meuwissen & Goddard (2007) showed that LDLA was more powerful for QTL detection and rendered less false positive results than single marker and haplotype association analyses. Those authors stated that the advantage of LDLA is that it models ‘the inheritance of chromosome segments without recombination from a common ancestor’. LA also does that, but the ancestors are recent and LD is conserved over many cM compared to considering more distant ancestors in LDLA. Finally, Lee & van der Werf (2005) reported that pedigree information did not change QTL profiles obtained with LD mapping using saturated marker maps. So, once short-range (historical) disequilibrium is captured with, for example, dense SNP chips, the long-range disequilibrium observed within pedigrees is redundant.

The version of LDLA used in this study, gridLDLA, assumes a Fisher–Wright model of population evolution defined by discrete generations, constant size and genetic drift. Our study population follows this evolutionary model very approximately. For example, the population of Hirta was founded in 1932 and has existed since without immigration. The average Ne was calculated as ~200 (McRae, 2006). In the absence of better information, the homozygosity levels among population founders were assumed to be equal to that in the pedigree founders. We chose T=75 because although a generation interval of 3 years from 1932 and 2007 equates to T of ~25, it is likely that population founders were already related, in which case a larger T is more appropriate. A sensitivity analysis changing Ne and T showed that gridLDLA was robust to changes in these parameters. Nevertheless, it is important to confirm the results of LDLA using several histories and to contrast those results against LA. Moreover, it is possible that a single set of global parameters does not perform well in all regions, because there may be variation across those regions regarding selection intensity and genetic drift. Finally, it is possible to extend the mathematical models to predict IBD under other evolutionary scenarios. For example, Meuwissen & Goddard (2007) developed theory to calculate IBD probabilities under a mutation-drift equilibrium model, thus eliminating the need to estimate, or guess, T. However, QTLs might not be in equilibrium or proving that they are may be difficult. GridQTL could also incorporate theory to accommodate variation in Ne and migration (Hill & Hernández-Sánchez, 2007). More complex population histories may be accommodated using simulations, e.g. overlapping generations.

The first objective of this work was to compare our results against those in Beraldi et al. (2007) . However, note that we did not analyse the exact same data and applied slightly modified models, e.g. we estimated maternal environmental rather than genetic effects. These authors used the software LOKI to obtain IBD probabilities for LA. LOKI uses MCMC to obtain the probability distribution of full haplotypes including loci with missing genotypes or ambiguous phases. In contrast, gridLA implements a deterministic procedure that increases computational speed at the expense of only using loci for which phase is completely known. Theoretically, given a perfect data set, i.e. no missing genotypes and perfectly known phases, both approaches should calculate the same IBD probabilities. However, we observed small differences in QTL profiles between LOKI and gridLA when analysing birth weight on C2 despite using exactly the same data and phenotypic models (Fig. 1). These differences were entirely due to differences in IBD probabilities. This is expected given that 38% of all genotypes were missing (some ungenotyped individuals were included to link different sibships) and that phases among known genotypes could be perfectly reconstructed in only ~80% of the cases. Nevertheless, differences were not substantial.

We found that gridLDLA was generally more powerful than gridLA at finding QTL. For example, the highest LOD scores obtained with gridLDLA for birth weight (n=550), hind leg length in adults (n=1757) and jaw length (n=500) were 2·81, 2·34 and 3·83, respectively, compared to 1·06, 0·49 and 2·81 obtained with gridLA. Nevertheless, a suggestive QTL (LOD=2·02) on C11 for hind leg length in all ages (n=2150) was only detected with gridLA. Simulation studies showed that LDLA is on an average more powerful than LA under a correct population model (Lee & van der Werf, 2004), although the contrary was true in a real data analysis (Meuwissen et al., 2002).

GridLDLA found a significant QTL affecting jaw length on C11 (n=500) and two suggestive QTL affecting birth weight on C2 and C8 (n=550). The latter QTL was also detected with gridLA. Beraldi et al. (2007) reported the QTL on C8 for birth weight, but did not detect the one on C2. Finally, a suggestive QTL affecting hind leg length in adults (n=1757) was mapped on C19. Beraldi et al. (2007) also reported a significant QTL for jaw length on C11 (LOD=3·59, n=566), but at a different position (37 cM) and suggestive QTLs for birth weight on C8 (LOD=2·54, n=601) at the same position (130 cM) and hind leg length in adults but on C15 (LOD=2·89, n=1954). Our sample sizes were smaller than those in Beraldi et al. (2007) as we excluded, unlike them, records for which fixed effects, e.g. litter size, were missing.

A very interesting result was found when analysing birth weight with gridLDLA using a subset (n=418) of the original data (n=601) that excluded records without TYRP1 genotypes or missing factors but included lambs older than 4 days with TYRP1 genotypes and no missing factors. A significant QTL was found with gridLDLA at 98 cM on C2 (LOD=3·9) but not with gridLA (LOD=1·4). Gratten et al. (2008) previously found a significant effect of a non-synonymous point mutation at TYRP1 (C2) on birth weight using association tests (n=1757), excluding the possibility of false associations with the transmission disequilibrium test (TDT; n=492) (Hernández-Sánchez et al., 2003). TYRP1 is completely linked to marker BMS0678 located at 96·8, and thus only 1·2 cM away from the peak detected with gridLDLA at 98 cM. A high-resolution scan every 1 cM between 90 and 100 cM rendered a plateau of significant LOD scores when using gridLDLA with Ne=200 and T=75, where the maximum LOD was 4·2 at 97 cM, just 0·2 cM to the right of TYRP1 (Fig. 4 B). The two closest markers in this region, OarFCB0128 at 90·9 cM and BMS0678 at 96·8 cM, showed the strongest LD on C2 (results not shown). The relatively long-range LD that still exists in this population undoubtedly contributed to detecting this QTL as the marker density is low (average distance between consecutive markers on C2 is 14·6 cM with a standard deviation of 7·9 cM). Gratten et al. (2008) hypothesized that a novel QTL affecting birth weight must be tightly linked to TYRP1 as this gene is only known to affect coat colour. However, they could not provide a better estimate of QTL position, because they only applied single-marker tests. Moreover, Beraldi et al. (2007) did not detect this QTL, probably because they only used LA. Only gridLDLA could detect and provide a good estimate of the position of this QTL. Adding TYRP1 genotype as a fixed factor explained almost all of the QTL variation, and thus we could not distinguish between pleiotropy and linkage. The additional 17 lambs older than 4 days contributed considerably to detecting this QTL. A potential explanation is that this QTL may be active perinatally rather than inside the womb. Gratten et al. (2008) included lambs between 0 and 4 months of age in their analyses.

In conclusion, LDLA can be more powerful than LA at detecting QTL as long as the assumptions underlying the historical model in LDLA are approximately fulfilled. Fortunately, the history of Soay sheep on St. Kilda is well documented and population parameters such as Ne have been accurately estimated. More importantly, the population has remained closed since its foundation. Migration can reduce the power of LDLA and is not yet accounted for in the gridLDLA models. Thus, LDLA can be an invaluable tool for geneticists interested in QTL mapping in wild as well as managed populations.