### Description of animals and genotyping data

The Kinsella beef composite population has been described in our recent studies on QTL mapping, candidate gene identification and genomic selection [e.g., 14-17]. Here we recapitulate the essential details of this population. It was produced by crossing between Angus, Charolais, or University of Alberta hybrid bulls and a hybrid dam line. The hybrid dam line was obtained by crossing among three composite cattle lines, namely beef synthetic 1 (SY1), beef synthetic 2 (SY2) and dairy × beef synthetic (SD) for more than 10 years after 30 years (1960-1990) single-sire crossbreeding. SY1 was composed of approximately 33% each of Angus and Charolais, 20% Galloway, 5% Brown Swiss, and small amounts of other breeds. SY2 was composed of approximately 60% Hereford and 40% other beef breeds mainly including Augus, Charolais and Galloway. SD was composed of approximately 60% dairy cattle (Holstein, Brown Swiss, or Simmental) and approximately 40% of other breeds, mainly including Angus and Charolais [22]. The blood samples of 1023 beef steers were collected and genotyped using the Illumina Infinium genotyping system with the BovineSNP50 Beadchip. All steers were produced from multi-sire breeding group natural service on pasture. The sire genotype of each calf was determined in a parentage test by using the BovineSNP50 Beadchip, but the parentage of about 100 animals remained unknown because these animals were either sires at initial crossing or sires without progeny. There were 116 sire families with varying family sizes ranging from one to 54 progeny per family. It is estimated that there have been about 4-5 generations since initial crossing.

A total of 51,828 SNP markers were originally obtained in the genotyped data. These markers were distributed across 29 autosomes and one sex chromosome in the entire bovine genome. For our analyses, we only used 43,124 SNPs after removing those markers (i) with monomorphism, (ii) with unknown genomic position and (iii) on the sex chromosome, (iv) with minor allele frequency (MAF) of ≤ 2% [1], and (v) with a Chi-square value >600 for the HWD test.

### Components of zygotic linkage disequilibrium

For two loci, each with two alleles, *A* and *a* at locus *A* and *B* and *b* at locus *B*, there are nine possible genotypes (ten if the coupling and repulsion double heterozygotes are distinguishable). Following Yang [23], we wrote frequencies of these genotypes as, {P}_{\mathit{vz}}^{\mathit{uy}}={P}_{\mathit{uy}}^{\mathit{vz}}, which result from union of gametes *uy* and *vz* with *u* *v* = *A* or *a*, and *y* *z* = *B* or *b*. The genotypic frequencies at individual loci are the marginal totals of the appropriate two-locus genotypic frequencies. For example, the frequency of genotype *AA* is,

{P}_{A\cdot}^{A\cdot}={P}_{\mathit{AB}}^{\mathit{AB}}+{P}_{\mathit{Ab}}^{\mathit{AB}}+{P}_{\mathit{Ab}}^{\mathit{Ab}}\text{.}

(1)

With the genotypic frequencies at locus *A*, the frequency of allele *A* is, {p}_{A}={P}_{A\cdot}^{A\cdot}+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.{P}_{a\cdot}^{A\cdot} and that of allele *a* is *p*_{
a
} = 1- *p*_{
A
}.

Departures from HWE at locus *A* are, {D}_{A}={P}_{A\cdot}^{A\cdot}-{p}_{A}^{2}={P}_{a\cdot}^{a\cdot}-{p}_{a}^{2}=-\left(\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.{P}_{a\cdot}^{A\cdot}-{p}_{A}{p}_{a}\right) and those at locus *B* are, {D}_{B}={P}_{\cdot B}^{\cdot B}-{p}_{B}^{2}={P}_{\cdot b}^{\cdot b}-{p}_{b}^{2}=-\left(\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.{P}_{\cdot b}^{\cdot B}-{p}_{B}{p}_{b}\right).

In a random mating population, HWD disappears (i.e., *D*_{
A
} = *D*_{
B
} = 0). In a non-random mating population, nonzero HWD is measured by the fixation index which can be either positive when there is inbreeding or negative when inbreeding is avoided. For example, the HWD at locus *A* can be written as *D*_{
A
} = *f*_{
A
}*p*_{
A
}*p*_{
a
}, where {f}_{A}=\frac{{P}_{A\cdot}^{A\cdot}-{p}_{A}^{2}}{{p}_{A}{p}_{a}}=\frac{{P}_{a\cdot}^{a\cdot}-{p}_{a}^{2}}{{p}_{A}{p}_{a}}=1-\frac{{P}_{a\cdot}^{A\cdot}}{2{p}_{A}{p}_{a}}, is the fixation index at locus *A*, with -1 ≤ *f*_{
A
} ≤ +1.

It was established [11, 12] that the total zygotic LD between loci *A* and *B* could be defined in terms of zygotic LDs for individual genotypes with each zygotic LD being a complex function of digenic, trigenic and quadrigenic disequilibria. For example, the zygotic LD for double homozygote *AABB* ({\omega}_{\mathit{AB}}^{\mathit{AB}}) would simply be the deviation of the frequency of double homozygote from the product of the corresponding homozygotes at loci *A* and *B*

{\omega}_{\mathit{AB}}^{\mathit{AB}}={P}_{\mathit{AB}}^{\mathit{AB}}-{P}_{A\cdot}^{A\cdot}{P}_{\cdot B}^{\cdot B}=2{p}_{A}{D}_{.B}^{\mathit{AB}}+2{p}_{B}{D}_{A.}^{\mathit{AB}}+2{p}_{A}{p}_{B}{D}_{..}^{\mathit{AB}}+2{p}_{A}{p}_{B}{D}_{.B}^{A.}+{({D}_{..}^{\mathit{AB}})}^{2}+{({D}_{.B}^{A.})}^{2}+{D}_{\mathit{AB}}^{\mathit{AB}}

(2)

where each genic disequilibrium (*D*) is the deviation of a frequency from that based on random association of genes and accounting for any lower order disequilibria. The usual gametic LD ({D}_{..}^{\mathit{AB}}) would be the deviation of frequency of gamete *AB* from the product of frequencies of allele *A* at locus *A* and allele *B* at locus *B*{D}_{..}^{\mathit{AB}}={P}_{..}^{\mathit{AB}}-{p}_{A}{p}_{B} with {P}_{..}^{\mathit{AB}}={P}_{\mathit{AB}}^{\mathit{AB}}+{P}_{\mathit{Ab}}^{\mathit{AB}}+{P}_{\mathit{aB}}^{\mathit{AB}}+{P}_{\mathit{ab}}^{\mathit{AB}}.

When zygotes arise from random union of gametes as often assumed in most LD studies, all non-gametic disequilibria including HWD would disappear (*e.g.*, {D}_{A.}^{A.}={D}_{.B}^{A.}={D}_{.B}^{\mathit{AB}}={D}_{\mathit{AB}}^{\mathit{AB}}=0). In this case, the zygotic LD for genotype *AABB* ({\omega}_{\mathit{AB}}^{\mathit{AB}}) would reduce to, {\omega}_{\mathit{AB}}^{\mathit{AB}}=2{p}_{A}{p}_{B}{D}_{..}^{\mathit{AB}}+{({D}_{..}^{\mathit{AB}})}^{2}.

This formula is the basis for possible use of double homozygosity to measure gametic LD in a random mating population [24, 25].

Since the two types of double heterozygote (*AB*/*ab* and *Ab*/*aB*) in our unphased SNP data could not be distinguished, we used the composite LD (Δ_{
AB
}) and a composite quadrigenic component (Δ_{
AABB
}) in place of gametic and quadrigenic disequilibria. Thus, the zygotic LD for genotype *AABB* ({\omega}_{\mathit{AB}}^{\mathit{AB}}) in equation (1) was rewritten as

{\omega}_{\mathit{AB}}^{\mathit{AB}}={P}_{\mathit{AB}}^{\mathit{AB}}-{P}_{A\cdot}^{A\cdot}{P}_{\cdot B}^{\cdot B}=2{p}_{A}{D}_{\mathit{ABB}}+2{p}_{B}{D}_{\mathit{AAB}}+2{p}_{A}{p}_{B}{\Delta}_{\mathit{AB}}+{\Delta}_{\mathit{AB}}^{2}+{\Delta}_{\mathit{AABB}}

(3)

where

{\Delta}_{\mathit{AB}}={P}_{\cdot \cdot}^{\mathit{AB}}+{P}_{\cdot B}^{A\cdot}-2{p}_{A}{p}_{B}={D}_{\cdot \cdot}^{\mathit{AB}}+{D}_{\cdot B}^{A\cdot}

(4)

and

{\Delta}_{\mathit{AABB}}={D}_{\mathit{AB}}^{\mathit{AB}}-2{D}_{\cdot \cdot}^{\mathit{AB}}{D}_{\cdot B}^{A\cdot}

(5)

It should be noted from equations (1) and (2) that the two trigenic disequilibria in (2) were rewritten without superscripts for notational simplicity.

### Maximum likelihood estimation

Following Weir and Cockerham [19] and Weir [20], we used the procedure of statistical inference based on the assumption of multinomial sampling of individual diploids from a population. The observed frequencies and disequilibria with tildes (~) were maximum likelihood (ML) estimates of corresponding parametric values. Since the additive models described earlier allowed for defining the same number of parameters as there would be degrees of freedom, the ML estimates were simply replacing all parametric values of frequencies and disequilibria with corresponding observed values. For example, the ML estimates of composite LD were simply given by, \begin{array}{lll}{\tilde{\Delta}}_{\mathit{AB}}\hfill & =\hfill & {\tilde{P}}_{\cdot \cdot}^{\mathit{AB}}+{\tilde{P}}_{\cdot B}^{A\cdot}-2{\tilde{p}}_{A}{\tilde{p}}_{B}\hfill \\ =\hfill & {\tilde{D}}_{\cdot \cdot}^{\mathit{AB}}+{\tilde{D}}_{\cdot B}^{A\cdot}\hfill \end{array}

However, the ML estimates might be biased because they would involve quadratic terms of multinomial variables. For example, the expectation of the squared gene frequency of allele *A* over replicate samples of size *n* would be, E({\tilde{p}}_{A}^{2})={p}_{A}^{2}+[{p}_{A}(1-{p}_{A})+{D}_{A}]/2n

where *D*_{
A
} is the HWD measure at locus *A*[20]. With the sufficiently large sample (*n* = 1023 animals) in our data set, we invoked large-sample theory for statistical inference about genic disequilibria. Thus, we ignored the possible biases of order 1/*n*.

### Hypothesis testing and power

With a ML estimate (\tilde{D} or \tilde{\Delta}) of a given genic disequilibrium *D* or *Δ*, along with its sampling variance, [Var(\tilde{D}) or Var(\tilde{\Delta})] being given in Appendix in the Additional file 1 section, we constructed a test statistic, {X}^{2}={\tilde{D}}^{2}/Var(\tilde{D}) or {X}^{2}={\tilde{\Delta}}^{2}/Var(\tilde{\Delta})

to test the hypothesis of zero disequilibrium (i.e., H_{0}: *D* = 0 or H_{0}: *Δ* = 0). Assuming the asymptotic normality of the ML estimate, *X*^{2} under the hypothesis of zero disequilibrium would be distributed as chi-square with one degree of freedom.

As usual, each chi-square test would commit two kinds of error: a true hypothesis may be rejected (Type I error) or a false hypothesis may not be rejected (Type II error). The probability of Type I error is measured by the significance level whereas the probability of Type II error is often related to the power of the test. Generally, as the power is the probability of rejecting a false hypothesis, it equals to one minus the probability of Type II error. In the present study, however, we adopted a different use of the power as proposed by Weir and Cockerham ([19], p. 100) and Weir ([20], p. 110): we calculated the power when the hypothesis being tested is true. In this particular case, a power value equals to the significance level.

### Chi-square statistic and correlation

In the past, the squared correlation (*r*^{2}) has been routinely used as a measure of gametic LD ({r}_{\mathit{GLD}}^{2}), composite LD ({r}_{\mathit{CLD}}^{2}), or zygotic LD ({r}_{\mathit{ZLD}}^{2}). We used a chi-square statistic ({X}_{i}^{2}=n{r}_{i}^{2}*i* = *GLD* *CLD* or *ZLG*) to test for the significance of the LD estimate. It is known from the literature [21] that the relationship of{X}_{i}^{2}=n{r}_{i}^{2} would hold exactly only for a 2 × 2 contingency table. This was the case for GLD and ZLD, but not for CLD. When dropping out three- and four-gene disequilibria in testing for zero composite LD ({\tilde{\Delta}}_{\mathit{AB}}=0), we would obtain an approximate chi-square statistic, {X}_{\mathit{AB}}^{2}\approx n{\tilde{\Delta}}_{\mathbf{A}B}^{2}/[({\tilde{\pi}}_{A}+{\tilde{D}}_{A})({\tilde{\pi}}_{B}+{\tilde{D}}_{B})].

Which would be equal to n{r}_{\mathit{CLD}}^{2} as given in Weir [26]. Similar approximations or restrictions would be needed if the relationship of{X}_{i}^{2}=n{r}_{i}^{2} were desired for three- and four-gene disequilibria. Thus, to avoid such approximations or restrictions, we used a generalized measure of square correlation {\phi}^{2}={X}^{2}/n[21] in place of *r*^{2} as a standardized measure of genic disequilibria. As pointed out above, the relationship of {\phi}^{2}={r}^{2} would hold only for a 2 × 2 contingency table.

### Data analysis

All data analysis and required computation were done using SAS 9.3 [27]. The calculation of gametic and composite LD was carried out using PROC ALLELE of SAS/Genetics 9.3. In this calculation, the SNP marker data was read in as columns of genotypes using the GENOCOL and DELIMITER= options in the PROC ALLELE statement; gametic LD was calculated if the HAPLO= EST option in the PROC ALLELE statement was invoked, whereas composite LD was calculated if the HAPLO= NONEHWD option was specified. Zygotic LD and its components as well as hypothesis testing were calculated using SAS Macro language and SAS/IML procedure.