3: Population Structure and Correlations Among Loci

Other approaches to population structure

There is a broad spectrum of methods to describe patterns of population structure in population genetic datasets. We’ll briefly discuss two broad-classes of methods that appear often in the literature: assignment methods and principal components analysis.

Assignment Methods

Here we’ll describe a simple probabilistic assignment to find the probability that an individual of unknown population comes from one of \(K\) predefined populations. For example, there are three broad populations of common chimpanzee (Pan troglodytes) in Africa: western, central, and eastern. Imagine that we have a chimpanzee whose population of origin is unknown (e.g. it’s from an illegal private collection). If we have genotyped a set of unlinked markers from a panel of individuals representative of these populations, we can calculate the probability that our chimp comes from each of these populations.

We’ll then briefly explain how to extend this idea to cluster a set of individuals into \(K\) initially unknown populations. This method is a simplified version of what population genetics clustering algorithms such as STRUCTURE and ADMIXTURE do.

Principal components analysis

Principal component analysis (PCA) is a common statistical approach to visualize high dimensional data, and used by many fields. The idea of PCA is to give a location to each individual data-point on each of a small number principal component axes. These PC axes are chosen to reflect major axes of variation in the data, with the first PC being that which explains largest variance, the second the second most, and so on. The use of PCA in population genetics was pioneered by Cavalli-Sforza and colleagues and now with large genotyping datasets, PCA has made a comeback.

Consider a dataset consisting of N individuals at \(S\) biallelic SNPs. The \(i^{th}\) individual’s genotype data at locus \(\ell\) takes a value \(g_{i,\ell}=0,1,\; \text{or} \; 2\) (corresponding to the number of copies of allele \(A_1\) an individual carries at this SNP). We can think of this as a \(N \times S\) matrix (where usually \(N \ll S\)).

Denoting the sample mean allele frequency at SNP \(\ell\) by \(p_{\ell}\), it’s common to standardize the genotype in the following way

\[\frac{g_{i,\ell} - 2 p_{\ell}}{\sqrt{2 p_{\ell}(1-p_{\ell})}} \label{eqn:std_allele_freq}\]

i.e. at each SNP we center the genotypes by subtracting the mean genotype (\(2p_{\ell}\)) and divide through by the square root of the expected variance assuming that alleles are sampled binomially from the mean frequency (\(\sqrt{2 p_{\ell} (1-p_{\ell})}\)). Doing this to all of our genotypes, we form a data matrix (of dimension \(N \times S\)). We can then perform principal component analysis of this data matrix to uncover the major axes of genotype variance in our sample. Figure \(\PageIndex{6}\) shows a PCA from using the same chimpanzee data as in Figure \(\PageIndex{5}\).

It is worth taking a moment to delve further into what we are doing here. There’s a number of equivalent ways of thinking about what PCA is doing. One of these ways is to think that when we do PCA we are building the individual by individual covariance matrix and performing an eigenvalue decomposition of this matrix (with the eigenvectors being the PCs). This individual by individual covariance matrix has entries the \([i,~j]\) given by

\[\frac{1 }{S-1} \sum_{\ell=1}^S \frac{(g_{i,\ell} - 2p_{\ell})(g_{j,\ell} - 2p_{\ell})}{2 p_{\ell}(1-p_{\ell})} \label{eqn:kinship_mat}\]

Note that this is the sample covariance of our standardized allele frequencies (\ref{eqn:std_allele_freq}), and is very similar to those we encountered in discussing \(F\)-statistics as correlations ( \ref{eqn:Fascorr}), except now we are asking about the covariance between two individuals above that expected if they were both drawn from the total sample at random (rather than the covariance of alleles within a single individual). So by performing PCA on the data we are learning about the major (orthogonal) axes of the kinship matrix.

As an example of the application of PCA, let’s consider the case of the putative ring species in the greenish warbler (Phylloscopus trochiloides) species complex. This set of subspecies exists in a ring around the edge of the Himalayan plateau. collected \(95\) greenish warbler samples from \(22\) sites around the ring, and the sampling locations are shown in Figure \(\PageIndex{7}\).

It is thought that these warblers spread from the south, northward in two different directions around the inhospitable Himalayan plateau, establishing populations along the western edge (green and blue populations) and the eastern edge (yellow and red populations). When they came into secondary contact in Siberia, they were reproductively isolated from one another, having evolved different songs and accumulated other reproductive barriers from each other as they spread independently north around the plateau, such that P. t. viridanus (blue) and P. t. plumbeitarsus (red) populations presently form a stable hybrid zone.

Alcaide et al. (2014) obtained sequence data for their samples at 2,334 snps. In Figure \(\PageIndex{9}\) you can see the matrix of kinship coefficients, using \ref{eqn:kinship_mat}, between all pairs of samples. You can already see a lot about population structure in this matrix. Note how the red and yellow samples, thought to be derived from the Eastern route around the Himalayas, have higher kinship with each other, and blue and the (majority) of the green samples, from the Western route, form a similarly close group in terms of their higher kinship.

We can then perform PCA on this kinship matrix to identify the major axes of variation in the dataset. Figure \(\PageIndex{10}\) shows the samples plotted on the first two PCs. The two major routes of expansion clearly occupy different parts of PC space. The first principal component distinguishes populations running North to South along the western route of expansion, while the second principal component distinguishes among populations running North to South along the Eastern route of expansion. Thus genetic data supports the hypoth- esis that the greenish warblers speciated as they moved around the Himalayan plateau. However, as noted by Alcaide et al. (2014), it also suggests additional complications to the traditional view of these warblers as an unbroken ring species, a case of speciation by continu- ous geographic isolation. The Ludlowi subspecies shows a significant genetic break, with the southern most MN samples clustering with the Trochiloides subspecies, in both the PCA and kinship matrix (Figure \(\PageIndex{10}\) and Figure \(\PageIndex{9}\)), despite being much more geographically close to the other Ludlowi samples. This suggests that genetic isolation is not just a result of geographic distance, and other biogeographic barriers must be considered in the case of this broken ring species.

Finally, while PCA is a wonderful tool for visualizing genetic data, care must be taken in its interpretation. The U-like shape in the case of the greenish warbler PC might be consistent with some low level of gene flow between the red and the blue populations, pulling them genetically closer together and helping to form a genetic ring as well as a geographic ring. However, U-like shapes are expected to appear in PCAs even if our populations are just arrayed along a line, and more complex geometric arrangements of populations in PC space can result under simple geographic models . Inferring the geographical and population-genetic history of species requires the application of a range of tools; see and for more discussion of the greenish warblers.

Correlations between loci, linkage disequilibrium, and recombination

Up to now we have been interested in correlations between alleles at the same locus, e.g. correlations within individuals (inbreeding) or between individuals (relatedness). We have seen how relatedness between parents affects the extent to which their offspring is inbred. We now turn to correlations between alleles at different loci.

Recombination

To understand correlations between loci we need to understand recombination a bit more carefully. Let us consider a heterozygous individual, containing \(AB\) and \(ab\) haplotypes. If no recombination occurs between our two loci in this individual, then these two haplotypes will be transmitted intact to the next generation. While if a recombination (i.e. an odd number of crossing over events) occurs between the two parental haplotypes, then \(\frac{1}{2}\) the time the child receives an \(Ab\) haplotype and \(\frac{1}{2}\) the time the child receives an \(aB\) haplotype. See Figure \(\PageIndex{11}\). Effectively, recombination breaks up the association between loci. For linked markers we’ll define the recombination fraction (\(x\)) to be the probability of an odd number of crossing over events between our loci in a single meiosis. The recombination fraction between a pair of loci can range from \(0\) to \(\frac{1}{2}\), with \(c=\frac{1}{2}\) corresponding markers far enough apart on a chromosome that many recombination events occur between them (loci on different automosomes also have a \(c=\frac{1}{2}\)). In practice we’ll often be interested in relatively short regions such that recombination is relatively rare, and so we might think that \(c=c_{BP}L \ll \frac{1}{2}\), where \(c_{BP}\) is the average recombination rate (in Morgans) per base pair (typically \(\sim 10^{-8}\) ) and L is the number of base pairs separating our two loci.

Linkage disequilibrium

The (horrible) phrase linkage disequilibrium (LD) refers to the statistical non-independence (i.e. a correlation) of alleles in a population at different loci . It’s a fantastically useful concept; LD is key to our understanding of diverse topics, from sexual selection and speciation to the limits of genome-wide association studies.

Our two biallelic loci, which segregate alleles \(A/a\) and \(B/b\), have allele frequencies of \(p_A\) and \(p_B\) respectively. The frequency of the two locus haplotype AB is \(p_{AB}\), and likewise for our other three combinations. If our loci were statistically independent then \(p_{AB} = p_Ap_B\), otherwise \(p_{AB} \neq p_Ap_B\) We can define a covariance between the \(A\) and \(B\) alleles at our two loci as

\[D_{AB} = p_{AB} - p_Ap_B \label{eqn:LD_def}\]

and likewise for our other combinations at our two loci (\(D_{Ab},~D_{aB},~D_{ab}\)). Gametes with two similar case alleles (e.g. A and B, or a and b) are known as coupling gametes, and those with different case alleles are known as repulsion gametes (e.g. a and B, or A and b). Then, we can think of \(D\) as measuring the excess of coupling to repulsion gametes. These \(D\) statistics are all closely related to each other as \(D_{AB} = - D_{Ab}\) and so on. Thus we only need to specify one \(D_{AB}\) to know them all, so we’ll drop the subscript and just refer to \(D\). Also a handy result is that we can rewrite our haplotype frequency \(p_{AB}\) as

\[p_{AB} = p_Ap_B+D. \label{eqn:ABviaD}\]

If \(D=0\) we’ll say the two loci are in linkage equilibrium, while if \(D>0\) or \(D<0\) we’ll say that the loci are in linkage disequilibrium (we’ll perhaps want to test whether \(D\) is statistically different from \(0\) before making this choice). Linkage disequilibrium is a horrible phrase, as it risks muddling the concepts of genetic linkage and linkage disequilibrium. Genetic linkage refers to the linkage of multiple loci due to the fact that they are transmitted through meiosis together (most often because the loci are on the same chromosome). Linkage disequilibrium merely refers to the covariance between the alleles at different loci; this may in part be due to the genetic linkage of these loci but does not necessarily imply this (e.g. genetically unlinked loci can be in LD due to population structure).

Our linkage disequilibrium statistic \(D\) depends strongly on the allele frequencies of the two loci involved. One common way to partially remove this dependence, and make it more comparable across loci, is to divide \(D\) through by its maximum possible value given the frequency of the loci. This normalized statistic is called \(D^{\prime}\) and varies between \(+1\) and \(-1\). In Figure \(\PageIndex{12}\) there’s an example of LD across the TAP2 region in human and chimp. Notice how physically close SNPs, i.e. those close to the diagonal, have higher absolute values of \(D^{\prime}\) as closely linked alleles are separated by recombination less often allowing high levels of LD to accumulate. Over large physical distances, away from the diagonal, there is lower \(D^{\prime}\). This is especially notable in humans as there is an intense, human-specific recombination hotspot in this region, which is breaking down LD between opposite sides of this region.

Another common statistic for summarizing LD is \(r^2\) which we write as

\[r^2 = \frac{D^2}{p_A(1-p_A) p_B(1-p_B) }\]

As \(D\) is a covariance, and \(p_A(1-p_A)\) is the variance of an allele drawn at random from locus \(A\), \(r^2\) is the squared correlation coefficient.

fraction.

Figure \(\PageIndex{14}\) shows \(r^2\) for pairs of SNPs at various physical distances in two population samples of Mus musculus domesticus. Again LD is highest between physically close markers as LD is being generated faster than it can decay via recombination; more distant markers have much lower LD as here recombination is winning out. Note the decay of LD is much slower in the advanced-generation cross population than in the natural wild-caught population. This persistence of LD across megabases is due to the limited number of generations for recombination since the cross was created.

Figure \(\PageIndex{16}\): The decay of LD from an initial value of D0 = 0.25 due
to recombination over t generations, plotted across possible recombination fractions (c) between our pair of loci. Code here.

The decay of LD due to recombination

We will now examine what happens to LD over the generations if, in a very large population (i.e. no genetic drift and frequencies of our loci thus follow their expectations), we only allow recombination to occur. To do so, consider the frequency of our \(AB\) haplotype in the next generation, \(p_{AB}^{\prime}\). We lose a fraction \(c\) of our \(AB\) haplotypes to recombination ripping our alleles apart but gain a fraction \(cp_A p_B\) per generation from other haplotypes recombining together to form \(AB\) haplotypes. Thus in the next generation

\[p_{AB}^{\prime} = (1-c)p_{AB} + cp_Ap_B \label{new_hap_freq}\]

The last term above, in \ref{new_hap_freq}, is \(c(p_{AB}+p_{Ab})(p_{AB}+p_{aB})\) simplified, which is the probability of recombination in the different diploid genotypes that could generate a \(p_{AB}\) haplotype.

We can then write the change in the frequency of the \(p_{AB}\) haplotype as

\[\Delta p_{AB} = p_{AB}^{\prime} -p_{AB} = -c p_{AB} + cp_Ap_B = - c D\]

[fig:LD_time]

[fig:LD_recom]

So recombination will cause a decrease in the frequency of \(p_{AB}\) if there is an excess of \(AB\) haplotypes within the population (\(D>0\)), and an increase if there is a deficit of \(AB\) haplotypes within the population (\(D<0\)). Our LD in the next generation is

\[\begin{aligned} D^{\prime} & = p_{AB}^{\prime} - p'_{A}p'_{B} \nonumber\\ & = (p_{AB} + \Delta p_{AB}) - (p_{A} + \Delta p_{A})(p_{B} + \Delta p_{B}) \nonumber\\ & = p_{AB} + \Delta p_{AB} - p_{A}p_{B} \nonumber\\ & = (1-c) D\end{aligned}\]

where we can cancel out \(\Delta p_{A}\) and \(\Delta p_{B}\) above because recombination only changes haplotype, not allele, frequencies. So if the level of LD in generation \(0\) is \(D_0\), the level \(t\) generations later (\(D_t\)) is

\[D_t= (1-c)^t D_0\]

Recombination is acting to decrease LD, and it does so geometrically at a rate given by \((1-c)\) . If \(c \ll 1\) then we can approximate this by an exponential and say that

\[D_t \approx D_0 e^{-ct} \label{eqn_LD_decay}\]

which follows from a Taylor series expansion, see Appendix \ref{eqn:Taylor_geo}.

Exercise \(\PageIndex{5}\)

You find a hybrid population between the two mouse subspecies described in the question above, which appears to be comprised of equal proportions (\(50/50\)) of ancestry from the two subspecies. You estimate LD between the two markers to be \(D=0.0723\). On the basis of previous work you estimate that the two loci are separated by a recombination fraction of 0.1. Assuming that this hybrid population is large and was formed by a single mixture event, can you estimate how long ago this population formed?

A particularly striking example of the decay of LD generated by the mixing of populations is offered by the LD created by the interbreeding between humans and Neanderthals . Neanderthals and modern humans diverged from each other likely over half a million years ago, allowing time for allele frequency differences to accumulate between the Neanderthal and modern human populations. The two populations spread back into secondary contact when humans moved out of Africa over the past hundred thousand years or so. One of the most exciting findings from the sequencing of the Neanderthal genome was that modern-day people with Eurasian ancestry carry a few percent of their genome derived from the Neanderthal genome, via interbreeding during this secondary contact . To date the timing of this interbreeding, looked at the LD in modern humans between pairs of alleles found to be derived from the Neanderthal genome (and nearly absent from African populations). In Figure \ref{fig:LD_Neanderthal} we show the average LD between these loci as a function of the genetic distance (\(c\)) between them, from the work of .

Assuming a recombination rate \(r\), we can fit the exponential decay of LD predicted by \ref{eqn_LD_decay} to the data points in this figure; the fit is shown as a red line. Doing this we estimate \(t=1200\) generations, or about 35 thousand years (using a human generation time of 29 years). Thus the LD in modern Eurasians, between alleles derived from the interbreeding with Neanderthals, represents over thirty thousand years of recombination slowly breaking down these old associations.