Almost all of the dairy cattle in Korea are Holstein, which have been produced by selected domestic seed bull or imported semen. Therefore, the Korean dairy cattle population in this manuscript is the Holsteins population in Korea. The country of origin for imported semen is shown in Table S4. The 96 Holstein samples were collected from four field dairy farms located in Gyeonggi and Gangwon provinces. The samples were collected by the Animal and Plant Quarantine Agency of Korea (a governmental agency) by random sampling for a survey on population disease resistance in Korea. The collected samples for this study represent the Korean Holstein cattle population despite the limitations of the sample size and is treated as a single population sample representing the dairy cattle of Korea. There are few regional differences in the dairy cattle in Korea because the supply and management of dairy cattle semen or seed bull takes place at the national level.

The Illumina BovineHD Genotyping BeadChip, which contains 777 kb SNP located across all autosomes, sex chromosomes, and mitochondria was used. These informative SNPs were selected from the Bovine Genome Database. Genotyping data were analyzed using Plink for quality control (GENO>0.05, MAF<0.1, HWE test p≤ 0.0001), which removed 226,156 SNP from the analysis because of poor genotyping quality (Purcell et al., 2007). As the Illumina BovineHD Genotyping BeadChip contains numerous SNP at a high density, the distance between adjacent SNPs is very short making the computing time long. To reduce the computational time, we used 70,000 randomly selected SNPs including only autosomal SNPs after quality control. Interbull (Uppsala, Sweden; www.interbull.org) uses BovineSNP50 Genotyping BeadChips, which contains 54,609 SNP to estimate genetic parameters and hence, we considered 70,000 SNPs sufficient for our purposes.

We estimated linkage disequilibrium using the “LDcorSV” package implemented in R. This package provides a set of functions to measure the existing r² and the r² corrected by the sample structure (Mangin et al., 2012). Estimating the r² of all pairwise per chromosome using high-density SNP chips is time-consuming. To reduce computing time, the input for each chromosome was split into files containing 100 loci, and r² was calculated for all syntenic marker pairs in each file as has been done in previous studies (Flury et al., 2010).

The standard measures of existing r² are respectively equivalent to the covariance and the correlation between alleles at two different loci (Hill and Robertson, 1968), computed as: where P_A, P_a, P_B, and P_b are the respective frequencies of alleles A, a, B and b, and D is P_AB−P_AP_B. Sample structure information is required for the corrected r². We used the genomic relationship matrix based on SNP data instead of pedigree data. In this study, the genomic relationship value of individual pairs was the number of common SNP between two individuals/number of total allele sites. In this way, we proposed a simple genomic relationship matrix (96×96). A process of calculating r² corrected by the genomic relationship structure based on SNP was identical to the existing r² calculation.

Details on the physical position of the markers can be found in the Illumina product literature (http://support.illumina.com/array/array_kits/bovinehd_dna_analysis_kit/downloads.ilmn). To determine linkage disequilibrium in relation to the physical distance between markers, marker pairs were divided into distance bins. We established two kinds of classes (0 to 0.5Mb, 0 to 5Mb) and subsequently, applicable marker pairs to each class were put into 50 distance bins with bin ranges dependent on the class (see Table S2). The two types of mean r² for each of the distance bins were then plotted against the median of the distance bin range (Mb). This estimation of r² was performed on a chromosome by chromosome basis; the pooled results are presented in Figure 2.

Also r² was calculated for a random selection of non-syntenic SNP. Across the autosome, 700 SNP were randomly selected and r² values were calculated for only non-syntenic marker pairs. This resulted in a total of 250,096 pairwise comparisons that were not corrected by the genomic relationship structure.

Under the assumption of an isolated population with random mating, Sved (1971) derived an approximate expression for the expectation of r² (Sved, 1971): where N is the effective population size, and c is the recombination frequency. In this study, as in previous studies, c was replaced by linkage distance in Morgans (Hayes et al., 2003; Tenesa et al., 2007; Thevenon et al., 2007; De Roos et al., 2008; Villa-Angulo et al., 2009; Corbin et al., 2010; Flury et al., 2010; Qanbari et al., 2010). This was justified by the approximation of the more precise equation for E(r²) given by Sved (1971). Based on this formula, a non-linear least-squares approach to statistically model the observed r² was implemented within R (nls function) using this model: where y_i is the value of r² for marker pair i, at linkage distance d_i in Morgans. Parameters a and b were estimated iteratively using the least-squares method. Chromosome-specific megabase-to-centimorgan conversion rates were calculated based on total physical chromosome length as stated on the UCSC Web site, and total chromosome genetic length derived from the bovine linkage map (Arias et al., 2009) (see Table S2). We used marker pairs for which r² values were higher than the mean of r² for non-syntenic marker pairs. This model was applied to each chromosome in turn and the parameters were estimated. Similar to Corbin (2010), estimated parameters were combined by meta-analysis in R using an inverse variance method for pooling and random effects method based on the DerSimonian-Laird method (DerSimonian and Laird, 1986; Corbin et al., 2010).

Equation (2) can predict the effective population size at a given point in time, expressed as generations in the past (Hayes et al., 2003; De Roos et al., 2008). where N_T is the effective population size t generations ago, c is the distance between markers in Morgans, r_c² is the mean value of r² for marker pairs c Morgans apart, and c = (2t)⁻¹ when assuming linear growth (Hayes et al., 2003). As previously mentioned, marker pairs with linkage distances less than the mean of r² for non-syntenic marker pairs (<0.014248) were excluded from this analysis. To estimate N_T, the number of prior generations was selected and a suitable range for c was calculated. The binning process was designed to ensure sufficient marker pairs within each bin and obtain a representative r² mean (see Table S3). This process was performed for markers pooled across autosomes.

For the 96 genotyped individuals, we were able to establish genomic relationship structure per chromosome using SNP information of each chromosome. We proposed r² corrected by the genomic relationship structure per chromosome using “LDcorSV” package implemented in R (Mangin et al., 2012). A process of estimating N_e based on genomic relationship per chromosome and ancestral N_e was identical to that described above.

Of the 734,862 genotyped autosomal SNPs, 508,707 (69.2%) remained after quality control processes and 70,000 were selected for analysis. The number of SNPs per autosome remaining after filtering and selection ranged from 1,300 to 4,280 and was closely related to chromosome length and total number of SNP as shown in Figure 1. The minor allele frequency of remaining SNP followed a uniform distribution and averaged (±SD) to be 0.31±0.11. The average distance between marker pairs (±standard deviation) for this analysis was 1,202.77±932 kb, with the distance between markers ranging from 0.134 kb to 8,398 kb.

Linkage disequilibrium declined with increasing distances between marker pairs, as shown in Figure 2a and b. The most rapid decrease was seen over the first four bins, with the mean r² decreasing by about half. The mean existing r² decreased more slowly with increasing distance and was constant after 0.2 Mb of distance. The mean r² corrected by the genomic relationship structure based on SNP with distance between marker pairs was slightly less than that of the existing r². However, change patterns in the r² means with distance for both methods were similar. The decline in linkage disequilibrium showed slight differences with log-transformed distance (Figures S1 and S2). According to the existing r², 23,797 of the 3,385,800 marker pairs were in complete linkage disequilibrium; 12,225 of these were adjacent pairs. For r² corrected by the genomic relationship structure based on SNP, 23,794 of the 3,385,800 marker pairs were in complete linkage disequilibrium; 12,223 of these were adjacent pairs. The mean (±standard deviation) r² between random non-syntenic marker pairs was 0.014248±0.0197. 835,324 using the existing r² measure and 861,676 using r² corrected by the genomic relationship structure based on SNP are less than the mean of non-syntenic marker pairs.

The non-linear regression model of the declining linkage disequilibrium with distance resulted in both parameters, a and b, being significantly different from zero. For parameter estimation using the existing r², the mean estimate and 95% confidence interval by meta-analysis across autosomes for parameters a and b were 2.95 (2.84; 3.07) and 106.32 (92.95; 119.70), respectively. The parameters a and b for r² corrected by the genomic relationship structure based on SNP were 2.87 (2.75; 2.99) and 122.28 (107.21; 137.34), respectively. The predicted r² from the non-linear regression equation was similar to the mean observed r² of regions with massed bins with discrepancies occurring in other regions for both two types of r² (Figures S1 and S2). Parameter b showed greater variability between chromosomes than parameter a. No such relationship was observed between the estimated parameters (a, b) and chromosome length (cM) as shown in Figure 3 and 4. This reason for the lack of relationship and interpretation of parameter b in this non-linear regression model as an estimated N_e is demonstrated in the discussion.

We observed rapid increase in N_e over the past 10 generations with the values increasing fivefold (close to 500) by 10 generations ago (Figure 5). From the past 10 generations to the past 100 generations, N_e increased slowly. Our results suggest that a continuous increase in N_e occurred over the last 100,000 generations (Figure S3). Overall, the tendencies for the two different r² were similar, but N_e based on r² corrected by the genomic relationship structure based on SNP was slightly higher than N_e based on the existing r². In Figure 5, although the recent difference between the two estimated measures was small, the difference in value increased with time.

Linkage disequilibrium corrected by the genomic relationship structure per chromosome declined with increasing distance between marker pairs, as shown in Figure S4a and b. Over the first four bins, the decline was greater than two previous results. The mean r² corrected by the genomic relationship structure per chromosome with distance between marker pairs was less than the two previous mean but the change pattern was similar. (Figure S5) In total, 23,783 of marker pairs were in complete linkage disequilibrium; 12,220 of these were adjacent pairs. For parameter estimation using r² corrected by the genomic relationship structure per chromosome, the mean estimate and 95% confidence interval by meta-analysis across autosomes for parameters a and b were 2.11 (2.04; 2.18) and 361.73 (335; 339.47), respectively. Parameter b was approximately three times greater than the two previous values for parameter b. The decline in linkage disequilibrium was almost linear with log-transformed distance and the decline in linkage disequilibrium with log-transformed distance was better fitted with predicted values than that of the other two. No such relationship was observed between parameters (a, b) and chromosome length (cM) in r² corrected by the genomic relationship structure per chromosome as shown in Figure S6. Similar to previous patterns, we observed a rapid increase in N_e over the past 10 generations and then a slow increase up to 100 generations ago. However, the values increased fivefold (close to 1,500) at 10 generation than two kinds of r² (Figure S7). From 10 generation ago, our results suggest that continuous increase in N_e has occurs over the 100,000 generations. (Figure S8) Overall, these tendencies were similar with that of the previous two but N_e based on r² corrected by the genomic relationship structure per chromosome was greater than the two previous kinds of r² measures. Although recent difference with two estimated measures was big, their difference in value decreased with generations ago.