Two diallelic loci A and B, with alleles A₁ and A₂ at the A locus and alleles B₁ and B₂ at the B locus, are considered. To incorporate selection and MAF threshold, we consider the A locus as the causal variant under selection and B locus as the molecular marker under the MAF filtering. Four possible haplotypes (i.e., A₁B₁, A₁B₂, A₂B₁, and A₂B₂) are possible at these two loci. In a population of size N_e, the frequency counts of these four haplotypes can take on k=(2Ne+3)!3!(2Ne)! possible values, where for each of these possibilities, the sum of the four counts is 2N_e. For example, when N_e is equal to two, all possible combinations of haplotype frequency counts are given in Figure S1. In general, let X be a k × 4 matrix with each row representing a possible combination of haplotype frequency counts for some value of N_e. Thus, the rows of X represent the collection of k lines with different allele frequencies at two loci. Let P_t denote a k × 1 vector with element i indicating the probability of the frequency counts in row i of X at generation t. Thus, P_t gives the distribution of allele frequencies at generation t. We show below that the distribution in generation t+1 can be written in general as

where A is a k × k transition matrix that is derived below in various circumstances of recombination, mutation, and selection.

First, consider a line with haplotype frequency counts xiT from ith row of X at generation t, e.g., xiT=[f11,f12,f21,f22], where the frequency counts for the four haplotypes, A₁B₁, A₁B₂, A₂B₁, and A₂B₂, are denoted by f₁₁, f₁₂, f₂₁ and f₂₂. Ignoring the effects of recombination, mutation, or selection, sampling of 2N_e gametes from this line can be modeled by a multinomial process with sample size n = 2N_e and probability vector θi=xiT2Ne for the four haplotype probabilities. Thus, the distribution of the frequency count in the next generation is given by the k × 1 vector m_i, where the element j of m_i is the probability of getting xjT from the Multinomial(n, θ_i) distribution. Given that element i of P_t is the probability of xiT in generation t, the distribution of allele frequencies in the next generation is given by:

where M is the k × k matrix with columns m_i, as described above, for i = 1, …, k. The above formula shows how the distribution of allele frequencies change due to genetic drift, ignoring recombination, mutation, and selection. Now we accommodate recombination in computing P_t+1. In gametes produced from generation t, the probability of a non-recombinant A₁B₁ is (1-c)×f112Ne, where c is the recombination rate between locus A and locus B. On the other hand, a recombinant A₁B₁ can be produced in one of four ways. They and their associated probabilities are:

Thus, accounting for recombination, the probability of a A₁B₁ haplotype changes from f112Ne to:

Similarly, due to recombination, the probabilities of the other three haplotypes become:

and

Now, to see how mutation alters these probabilities, we let θi* be a vector of the four probabilities from Equations (8) through (11) computed using the four haplotype frequencies in xiT. In modeling mutation, we assume that an A₁ or B₁ allele mutates to an A₂ or B₂ allele with probability u and an A₂ or B₂ allele mutates to an A₁ or B₁ allele with probability v. Then, haplotype probabilities following mutation can be computed as:

where

Furthermore, to incorporate selection in the model, a selection coefficient s that reduces the allele frequency of A₁ at locus A, which is a causal variant, is considered. Conditional on the haplotype probabilities in θi**, i.e., pA1B1**, pA1B2**, pA2B1**, and pA2B2**, the haplotype probabilities following selection can be computed as

Finally, to compute the distribution of allele frequencies in generation t+1 using Equation (6), where the forces of the genetic drift, recombination, mutation and selection are simultaneously considered, the matrix A is defined such that element j of column i contains the probability of getting xjT from a Multinomial (2Ne,θi***T) distribution, for i, j = 1, …, k.

The value of r² can be computed for a sub-population or line with haplotype frequency counts in any row of X. For example, consider the frequency counts in row 11 of the X matrix given in Figure S1, where f₁₁ = 0, f₁₂ = 2, f₂₁ = 1 and f₂₂ = 1. From these frequency counts, Pr(A₁) = 1/2, Pr(B₁) = 1/4, and Pr(A₁B₁) = 0 are obtained, and r² calculated using Equations (1) and (2) is 1/3 for a line with haplotype frequency counts in x11T from Figure S1. The distribution of r² in generation t is given by values of r² corresponding to the frequency counts in each row of X together with the probabilities for haplotype frequency counts given in P_t. Note that haplotype frequency counts in rows of X that have indeterminate values of r² (i.e., when one or other allele is extinct) are not used to compute the distribution of r². Similarly, the distribution of r² with MAF cutoff is given by considering only the values of r² and their probabilities corresponding to the rows of X with MAF at the B locus ≥5%. MAF of 0.05 is used as a threshold in the present study.

Starting with an allele frequency of 0.5 at each locus and linkage equilibrium between the two loci, the expected value of r² was computed over generations given some values of the mutation rate, recombination rate, selection coefficient and effective population size. Mutation rate of u = v = 0 is used to represent the absence of mutation, while u = v = 1 × 10⁻⁹ is used to represent the existence of mutation. Similarly, a selection coefficient of s = 0 is used to represent the absence of selection, while s = 0.1 or s = 0.01 is used to represent the existence of selection.

In our analysis, different population parameters are considered. Populations with N_e ranging from 5 to 50 in intervals of 5 are compared under different recombination rates, mutation rates (i.e., u = 0 or u = 1e⁻⁹) and selection coefficients for the causal variant (i.e., s = 0.1 or s = 0.01 for locus A), in the absence or presence of filtering by MAF (MAF≥0.05 for locus B). The recombination rates ranged from 0.01 to 0.1 in intervals of 0.01 and from 0.1 to 0.5 in intervals of 0.1.

The recombination rate for adjacent markers in a bovine 777 k SNP panel was considered to mimic a realistic scenario. The recombination rate for adjacent markers was estimated to be 6.25 × 10⁻⁵ using the Kosambi map function given that the mean distance between adjacent markers is around 5 kb (Espigolan et al., 2013), and the average genetic distance per unit of physical distance in bovine genome is 1.25 cM/Mb (Arias et al., 2009). A period of more than 3,000 generations was used to ensure that the haplotype frequencies had reached their equilibrium values. At equilibrium, the distribution of LD stays constant over generations (i.e., P_t+1 = P_t). That distribution was used to describe, quantify, and compare rE2 at different equilibria in the absence or presence of mutation, selection or MAF filtering. Results from a population similar to the international black and white Holstein dairy cattle are presented in these cases, for which the effective population size is estimated to be about 50 (Kim and Kirkpatrick, 2009; Wray et al., 2019).

Furthermore, the exact distributions of rE2 derived under different scenarios were used to validate the approximate deterministic formulas from Sved and Feldman (1973) and Hill (1975). To generalize our results to larger effective population size, a non-linear regression model of the form 1a+b×Ne×c, following Sved's formula, was considered. The parameters in that model were estimated by non-linear least squares using the data generated from our transition-matrix approach.

The authors state that all data necessary for confirming the conclusions presented in the article are represented fully within the article.

In this section, Sved's and Hill's approximations were compared to the exact distribution of rE2 derived from our transition-matrix approach. The relationship between equilibrium value of r² (rE2), recombination rate and mutation rate is shown in Figure 1 for these three approaches. In the absence of mutation and selection, Sved's formula showed consistency with the exact values from our transition-matrix approach (Figure 1A). On the other hand, Hill's formula had a better fit to the exact values of rE2 (Figure 1B) in the presence of mutation. However, neither Sved's nor Hill's deterministic formulae are accurate enough to describe rE2 in the presence of mutation. Therefore, we proposed an adjusted non-linear regression model to correct Sved's approximation to predict rE2 for larger effective population sizes in the presence of mutation. The mean square of errors is used to evaluate the fit of these three models and shows that our calibrated non-linear regression model was significantly better than Sved's or Hill's approximations.

The trajectory of evolving LD over generations computed from our transition-matrix approach is shown in this section, in the absence or presence of mutation and selection, where selection was applied with either a low selection coefficient of 0.01 or a high selection coefficient of 0.1. Starting with an allele frequency of 0.5 at each locus and linkage equilibrium between the two loci, expected values of LD [i.e., E(r²)] over generations are displayed in Figure 2. Note that the trajectory may be quite different, depending on the initial conditions used.

In the absence of mutation and selection, the expected value of r² over generations increases to its maximum and reaches its “equilibrium” status. Value of rE2 increases as c decreases due to the decreasing breakdown of LD with lower recombination rates. Typically, when the recombination rate is relatively small (e.g., c = 6.25 × 10⁻⁵), the expected value of LD at “equilibrium” is almost equal to 1. However, at this “equilibrium” stage, frequencies of lines (P_t) keep changing though the distribution of allele frequencies remain unchanged. This is because forces of recombination and drift are balanced in lines with determinate r², and frequencies of lines with determinate r² almost proportionally decrease toward 0.

In the absence of selection but with mutation, the expected value of r² increases and stays at an apparent equilibrium for several generations and then decreases to its true equilibrium value, where the mutation-drift equilibrium is reached. During the apparent equilibrium that is initially observed, the expected value of r² is identical to the equilibrium value in the absence of mutation and selection, where the forces of drift and recombination are balanced. When this stage is first reached, the effect of mutation is negligible because the mutation rate is low relative to the frequency of lines with segregating loci. The apparent equilibrium ends when the frequency is high for lines where one or both loci are fixed (r² is indeterminate) and the frequency of lines with segregating loci becomes close to the rate of mutation. Then, alleles introduced due to mutation into these lines have a noticeable effect on the distribution of allele frequencies and expected LD changes. At the end of the apparent equilibrium, the most frequent lines have only one haplotype. In the lines with two haplotypes, which have low frequencies, r² is either indeterminate or has value 1.0. When mutation introduces a third haplotype in a line where r² is 1.0, it drops in value. Further, when mutation introduces a third haplotype in a line where r² is indeterminate, the value of r² will be close to zero. Thus, as mutation becomes noticeable, the expected value of LD decreases. When the true equilibrium is reached, the loss of variation by fixation is balanced by its replenishment by mutation. In other words, the frequency of lines where loci are segregating will remain non-zero. Therefore, when mutation is present, the probabilities of allele frequencies (P_t) stay constant over generations when the true equilibrium is reached.

In the presence of both mutation and selection, lines with segregating loci become low in frequency at an earlier stage due to selection, particularly when the selection coefficient is large (e.g., s = 0.1), and thus, mutation starts to have a noticeable effect on allele frequency at an earlier stage in the presence of selection relative to when selection is absent. When selection coefficient is relatively large (e.g., s = 0.1), the expected value of r² reaches its maximum value early and then decreases sharply to two plateaus. When the recombination rate is low, e.g., c ≤ 0.01, the difference between these two plateaus is small.

In the absence of selection, allele frequencies at two loci tend to be similar. However, when locus A is under selection, allele frequencies at these two loci tend to be different such that LD is lower than when selection is absent. This effect, however, is negligible when the selection coefficient is small (e.g., s = 0.01). Thus, a selection coefficient of 0.1 is used to present results when mutation and selection are present in the next section.

The expected values of r² over generations in the presence of filtering by MAF (i.e., MAF ≥0.05 for locus B) are shown in Figure S2, where a similar pattern as for E(r²) over generations are observed.

In a population where mutation is present and N_e = 50, the equilibrium value of r² (rE2) is shown in Figure 3, when selection is absent or present, for different recombination rates. As explained in last section, the equilibrium values of LD in the absence of selection are always higher than those in the presence of selection. The equilibrium value of LD (i.e., rE2) with filtering by MAF (i.e., MAF ≥5% at locus B) is higher than its corresponding value in the absence of filtering.

To understand why filtering by MAF increases rE2, the k lines with different haplotype frequencies were divided into two groups: lines where MAF at the locus B is <5% and lines where MAF ≥5%. Next, the transition-matrix approach was used to compute the frequency of each of these lines at equilibrium. In lines where r² is defined, the equilibrium frequency of each line is plotted against its r² value in Figure 4; the color blue is used for the first group of lines where MAF at the locus B is <5% and red is used for the second group of lines where MAF is ≥5%; frequencies in the absence of selection are given in Figure 4A and those in the presence of selection are given in Figure 4B.

In the absence of selection (Figure 4A), a proportion of about 0.41 of the lines were in the first group and a proportion of about 0.59 were in the second group. It can be seen from this figure that most of the r² values in the first group were small, where around 45% of the lines had an r² value < 0.001. In contrast, the second group had many lines with large r² values; around 15% of the lines had r² = 1. This difference in the distribution of r² in the two groups shows why filtering by MAF at the B locus (removing lines from the first group, which had an abundance of low r² values) would increase the expected value of r².

The reason for the lower value of r² in the first group is that, due to filtering by MAF at the B locus, this group has a large proportion of lines with recent mutations at the B locus. Consider a line where alleles are segregating at the A locus but fixed at the B locus. Such a line would belong to neither of the groups because r² is not defined in a line where one locus is fixed. However, the introduction of a new allele at the B locus due to mutation will result in r² becoming defined, but, typically, at a very low level because it results from the association in a single haplotype. The MAF in this line with the new mutation will be 12Ne, and it will belong to the first group provided that 12Ne<0.05.

Figure 4B gives the distribution of r² for the two groups in the presence of selection, and it can be seen that still there is a greater abundance of low r² values in the first group. The difference between the two groups, however, is smaller than in the absence of selection, and this explains the smaller effect of filtering on equilibrium value of r² when selection is present (Figure 3).

The two existing deterministic formulas (i.e., Sved's and Hill's formula) for rE2 are not accurate as shown in Figure 1. Thus, a non-linear regression formula by recombination rate was calibrated using the exact values of rE2 calculated from our transition-matrix approach. rE2 is computed as the expectation of r² at equilibrium (i.e., mean r² weighted by the frequency of each of the k possible lines at equilibrium). In the following, this formula is referred to as the calibrated non-linear regression formula. The form of the non-linear regression formula follows that of Sved's formula as presented below:

To study the extrapolation of rE2 from small population sizes (i.e., N_e ≤ 50) to predict those for a larger population size, we split the values of rE2 calculated from our transition-matrix approach into training and validation sets. Exact values of rE2 with N_e ≤ 40 are included in the training set, and the remaining with N_e >40 are used for validation. The prediction accuracy is assessed using mean square error (MSE). As shown in Table S1, prediction accuracy from the calibrated non-linear regression formula is substantially higher than those from Sved's or Hill's formulas. Finally, parameters in Equation (13) are estimated using all available values of rE2, and parameter estimates are shown in Table 1. Typically, the estimates of β₀ and β₁ are inversely related to the values of recombination rates.

To further validate the accuracy of the extrapolation approach, the expected value of LD at equilibrium between adjacent markers in a 777k SNPs chip (i.e.,c = 6.25 × 10⁻⁵) was predicted for a Nellore cattle population, for which the effective population size is about 100 (Brito et al., 2013), using Equation (13) with β₀ = 2.26 and β₁ = 1557.11. In a real genotyped Nellore cattle population (Espigolan et al., 2013), the estimated mean and standard deviation for rE2 was 0.17 and 0.20, respectively. The predicted rE2 is 0.08, 0.449, and 0.976 using the calibrated non-linear regression formula, Hill's formula, and Sved's formula, respectively. Only the predicted rE2 from the calibrated non-linear regression formula falls in the range 0.17±0.2, while the other two are out of the measured range from this real population.

We further studied the accuracy of the extrapolation approach for very large effective population size, e.g., Ne = 10, 000 in human. However, rE2 did not align with the biological expectation, which indicates that our non-linear formulas for extrapolation of rE2 to extremely large population size does not work because the non-linear regression formulas are derived from data with N_e ≤ 50.

The maximum effective population size (N_e) presented in this paper is 50 (Kim and Kirkpatrick, 2009; Wray et al., 2019), which is the estimated N_e for the international black and white Holstein dairy cattle population. Larger N_e was not studied in this paper due to computational limitations. When N_e = 50, the transition matrix A is square of order k=4Ne2+8Ne=176,851, which requires more than 250 gigabytes memory to store. In addition, it takes around 5 h to complete the computational process for the analysis of 3,500 generations. Note that the memory complexity to store A is O(Ne6). When N_e = 100, the estimated processing memory to store the transition matrix A is more than 14 terabytes. This computational problem may be addressed by parallel computing. This idea, however, needs further investigation and is out of the scope of this paper. Thus, to generalize the transition-matrix approach a non-linear regression model was calibrated using exact values of rE2. The proposed calibrations may provide a better description of the relationship between LD, effective population size, recombination rate, and mutation.