A Million-Cow Validation of a Chromosome 14 Region Interacting with All Chromosomes for Fat Percentage in U.S. Holstein Cows

A genome-wide association study (GWAS) of fat percentage (FPC) using 1,231,898 first lactation cows and 75,198 SNPs confirmed a previous result that a Chr14 region about 9.38 Mb in size (0.14–9.52 Mb) had significant inter-chromosome additive × additive (A×A) effects with all chromosomes and revealed many new such effects. This study divides this 9.38 Mb region into two sub-regions, Chr14a at 0.14–0.88 Mb (0.74 Mb in size) with 78% and Chr14b at 2.21–9.52 Mb (7.31 Mb in size) with 22% of the 2761 significant A×A effects. These two sub-regions were separated by a 1.3 Mb gap at 0.9–2.2 Mb without significant inter-chromosome A×A effects. The PPP1R16A-FOXH1-CYHR1-TONSL (PFCT) region of Chr14a (29 Kb in size) with four SNPs had the largest number of inter-chromosome A×A effects (1141 pairs) with all chromosomes, including the most significant inter-chromosome A×A effects. The SLC4A4-GC-NPFFR2 (SGN) region of Chr06, known to have highly significant additive effects for some production, fertility and health traits, specifically interacted with the PFCT region and a Chr14a region with CPSF1, ADCK5, SLC52A2, DGAT1, SMPD5 and PARP10 (CASDSP) known to have highly significant additive effects for milk production traits. The most significant effects were between an SNP in SGN and four SNPs in PFCT. The CASDSP region mostly interacted with the SGN region. In the Chr14b region, the 2.28–2.42 Mb region (138.46 Kb in size) lacking coding genes had the largest cluster of A×A effects, interacting with seventeen chromosomes. The results from this study provide high-confidence evidence towards the understanding of the genetic mechanism of FPC in Holstein cows.


Introduction
The fat percentage (FPC) in Holstein cattle has the strongest genetic effects among dairy traits, and the gene of diacylglycerol O-acyltransferase 1 (DGAT1) of chromosome 14 (Chr14) has been widely confirmed to contain the most significant effects of FPC [1][2][3][4][5][6][7], including the effect of the K232A variant [3][4][5] and the effect of the single-nucleotide polymorphism (SNP) marker rs109421300 (ARS-BFGL-NGS-4939) in DGAT1 [1,6].As an example showing how much more significant the FPC additive effects are than the additive effects of other dairy traits, the log 10 (1/p) value as a measure of statistical significance was 5150 for FPC, and it was 1320, 820, 374 and 371 for the protein percentage, milk yield, fat yield and protein yield, respectively, from a previous large-scale genome-wide association study (GWAS) using 60,671 SNPs and 294,079 Holstein cows (2019 GWAS) [1].Although the exact reasons for the highly significant effect of DGAT1 on FPC were not completely understood, the antagonism between milk and fat yields of DGAT1 was a likely genetic mechanism [1,8,9], and the antagonism of the most significant SNP (rs109421300) was extreme, with the lowest milk yield and the highest fat yield among all SNPs on the cattle genome [1].In addition, a large chromosome segment (6.79 Mb in size) containing DGAT1

Results and Discussion
The genome-wide tests of inter-chromosome A×A effects for FPC using 1,231,898 first lactation cows and 75,198 SNPs essentially confirmed the previous result that the 10 Mb Chr14 region interacted with all chromosomes for FPC, confirmed the structure of this 10 Mb region and revealed many new results.This study identified 2763 pairs of significant inter-chromosome A×A effects with log 10 (1/p) > 32 for FPC (Figure 1, Table S1).Of these, 2761 pairs each involved the 0.14-9.52Mb Chr14 region (9.38 Mb in size) with 107 SNPs and one of the remaining 29 chromosomes with 1050 SNPs, whereas only two pairs did not involve the Chr14 region, one pair between Chr03 and Chr05 and one pair between Chr05 and Chr20 (Figure 1a).Among the 2763 inter-chromosome A×A effects, the four most significant effects were between Chr06 and Chr14 (Figure 1b).The total number of SNPs in the 2763 SNP pairs was 1160.The circular plots of A×A effects between Chr14 and other chromosomes, as well as the Manhattan plots of the statistical significance of the A×A effects of all non-Chr14 chromosomes, are provided in Figure S1.Detailed test results, including statistical significance and estimated A×A effects and values for each of the 2763 SNP pairs, are provided in Table S1.

Structure of the 9.38 Mb Chr14 Region Interacting with All Chromosomes for FPC
The significant inter-chromosome A×A effects in the 9.38 Mb Chr14 region (Figure 2a) were in two sub-regions that we named Chr14a and Chr14b, where Chr14a was the 0.14-0.88Mb region (0.74 Mb in size) with 19 SNPs [11], and Chr14b was the 2.21-9.52Mb region (7.31Mb in size) with 88 SNPs [12].These two sub-regions were separated by a 1.3 Mb gap at 0.9-2.2Mb, without significant inter-chromosome A×A effects with log10(1/p) > 32 (Figure 2b) but with over 40 coding genes [13].The Chr14a region had the largest number of significant inter-chromosome A×A effects, 2148 out of the 2761 pairs or 78%, interacting with all the remaining chromosomes, and the Chr14b region had 613 of the 2761 pairs or 22%, interacting with all the remaining chromosomes except Chr24 (Table S1).Details of Chr14a and Chr14b are described below.The significant inter-chromosome A×A effects in the 9.38 Mb Chr14 region (Figure 2a) were in two sub-regions that we named Chr14a and Chr14b, where Chr14a was the 0.14-0.88Mb region (0.74 Mb in size) with 19 SNPs [11], and Chr14b was the 2.21-9.52Mb region (7.31Mb in size) with 88 SNPs [12].These two sub-regions were separated by a 1.3 Mb gap at 0.9-2.2Mb, without significant inter-chromosome A×A effects with log 10 (1/p) > 32 (Figure 2b) but with over 40 coding genes [13].The Chr14a region had the largest number of significant inter-chromosome A×A effects, 2148 out of the 2761 pairs or 78%, interacting with all the remaining chromosomes, and the Chr14b region had 613 of the 2761 pairs or 22%, interacting with all the remaining chromosomes except Chr24 (Table S1).Details of Chr14a and Chr14b are described below.

Inter-Chromosome A×A Effects of Chr14a
The Chr14a region was a gene-dense area with at least 46 coding genes [11], including the 15 genes shown in Figure 2c, which shows multiple locations with large clusters of inter-chromosome A×A effects.An SNP upstream of LOC789384 (rs136939758) was the upstream boundary of Chr14a with 139 inter-chromosome A×A effects.An SNP in PLEC with one of the lowest log 10 (1/p) values (32.40) was the downstream boundary of Chr14a.
Two regions within Chr14a were particularly interesting: the PPP1R16A-FOXH1-CYHR1-TONSL (PFCT) region and the CPSF1-ADCK5-SLC52A2-DGAT1-SMPD5-PARP10 (CASDSP) region.The PFCT region was only 29 Kb in size with four SNPs, but had the largest number of inter-chromosome A×A effects (1141 pairs) with all chromosomes, including the most significant inter-chromosome A×A effects (Figures 1b and 2c, Table 1).The CASDSP region was 280 Kb in size, with six SNPs.Unlike PFCT, which interacted with all chromosomes, CASDPS mainly (42 out of 98 pairs) interacted with the SLC4A4-GC-NPFFR2 (SGN) region of Chr06 (86.75-87.40Mb).The epistasis effects between the SGN region of Chr06 and the PFCT and CASDSP regions of Chr14 are an important new discovery in this study.

Inter-Chromosome A×A Effects of Chr14a
The Chr14a region was a gene-dense area with at least 46 coding genes [11], including the 15 genes shown in Figure 2c, which shows multiple locations with large clusters of inter-chromosome A×A effects.An SNP upstream of LOC789384 (rs136939758) was the upstream boundary of Chr14a with 139 inter-chromosome A×A effects.An SNP in PLEC with one of the lowest log10(1/p) values (32.40) was the downstream boundary of Chr14a.

Inter-Chromosome A×A Effects between Chr06 and Chr14a
The significant inter-chromosome A×A effects of Chr06 (314 effects) were distributed in the 6.11-114.38Mb region (Figure 3a,b; Table S1).However, Chr06 nearly exclusively interacted with Chr14a, with 312 of the 314 significant inter-chromosome A×A effects between Chr06 and Chr14a.Only two A×A effects involved Chr14b, between an SNP in the ABCG2 gene of Chr06 and two SNPs in Chr14b (Figure 3a,b; Table S1), noting that ABCG2 had significant effect on milk yield [1].Although the 312 pairs between Chr06 and Chr14a were distributed in the 6.11-114.38Mb region, the A×A effects between Chr14a and the SGN region of chr06 were most interesting because the Chr14a and SGN regions were two of the most important regions affecting milk production, and the SGN region also had highly significant effects on somatic cell score, daughter pregnancy rate and cow conception rate [1,2,[14][15][16][17].
The PFCT region had four SNPs each interacting with the same eleven SNPs in SGN, and the CPSF1-ADCK5-SLC52A2-DGAT1 (CASD) region had three SNPs each interacting with the same nine SNPs in SGN (Figure 3c, Table S1).The SLC4A4 gene had ten significant SNPs and four of these ten SNPs in the 10.65 Kb tail region (86751807-86762457 bp), i.e., 2.5% of the gene interacted with both the PFCT and CASD regions of Chr14a, noting that the size of SLC4A4 was 427.295Kb.The GC gene had no significant inter-chromosome A×A effects for FPC.The GC-NPFFR2 region had seven significant SNPs and two of these seven SNPs interacted with the PFCT region (Table 1).The top four most significant interchromosome A×A effects were between an SNP in GC-NPFFR2 (rs42766480) of Chr06 and four SNPs in PFCT, rs110984572 in PPP1R16A-FOXH1 (#1), rs137472016 in CYHR1-TONSL (#2), rs137727465 in CYHR1 (#3) and rs10914637 in PPP1R16A (#4) (Table 1; Figure 3c,d).These results showed that rs42766480 likely had a major role in the interactions with the PFCT region.However, this SNP was not ranked high (highest ranking #742) among the A×A effects between the SGN and CASD regions.The most significant inter-chromosome A×A effect between CASD and SGN was between rs211309638 in ADCK5-SLC52A2 and rs110352004 in GC-NPFFR2, with a ranking of #37 among all effects.Other than this SNP, the effect rankings of SNPs in CPSF1, DGAT1 and SMPD5 were in the range of #153-#393 (Table 2), less significant than those of the PFCT region (ranking #1-#66, Table 1).
inter-chromosome A×A effects was that the SGN region only interacted with Chr14a, including the PCFT and CASD regions.These results of Chr14a and the SGN region of Chr06 were particularly interesting because Chr14a and the SGN region of Chr06 were two of the most significant regions for milk production traits and the SGN region also was highly significant for somatic cell score and two fertility traits (daughter pregnancy rate and cow conception rate) [1,14].The most important feature of PFCT was that this small 29 Kb region of Chr14a interacted with all chromosomes.The most important feature of CASD was that this region mainly interacted with the SGN region of Chr06.The most important feature of the Chr06 inter-chromosome A×A effects was that the SGN region only interacted with Chr14a, including the PCFT and CASD regions.These results of Chr14a and the SGN region of Chr06 were particularly interesting because Chr14a and the SGN region of Chr06 were two of the most significant regions for milk production traits and the SGN region also was highly significant for somatic cell score and two fertility traits (daughter pregnancy rate and cow conception rate) [1,14].

Other Inter-Chromosome A×A Effects of Chr14a
The Chr14a region had other highly significant inter-chromosome A×A effects, in addition to those with the SGN region of Chr06, including those with Chr02, Chr05, Chr17, Chr29 and ChrX (Table 3, Figure 4).An SNP in LOC789384 (rs109208977), an SNP between ZNF250 and ZNF16 (rs110508680) and an SNP in ARHGAP39 each had a large cluster of inter-chromosome A×A effects, with 260, 195 and 231 inter-chromosome A×A effects, respectively (Figure 2c, Table S1).Of the 20 highly significant effects not involving the SGN region of Chr06 (Table 3), 18 involved SNPs in the PCFT region, further showing a major role of the PCFT region in the inter-chromosome A×A effects for FPC.

Inter-Chromosome A×A Effects of Chr14b
The Chr14b region about 7.31 Mb in size (Figure 2b,d) [12] was nearly ten times as large as Chr14a (0.74 Mb in size) and had multiple locations with large clusters of interchromosome A×A effects.This region was divided into two sub-regions for convenience in describing the results: the 2.28-2.42Mb region (138.46Kb in size) as 'Chr14b1' and the remaining 7.17 Mb region as 'Chr14b2', with Chr14b1 accounting for 2% and Chr14b2 for 98% of Chr14b.
The Chr14b1 region had six significant SNPs in noncoding regions with two uncharacterized coding genes (LOC112449593, LOC112449592) and TRNAC-GCA [18].These six SNPs had the largest number of inter-chromosome A×A effects (119 pairs) in Chr14b, involving seventeen chromosomes (chromosomes 2, 3, 4, 8, 12, 13, 16, 17, 18, 20, 21, 23, 25, 26, 28, 29, X) (Table S1).The most significant inter-chromosome A×A effect of Chr14b (log 10 (1/p) = 64.37,#5 ranking) was between rs134537992 in Chr14b1 and rs42368654 about 96.6 Kb downstream of the LMX1A gene of Chr03 (Figure 5a, Table 4).This SNP of Chr03 interacted with 14 SNPs in Chr14b, two SNPs immediately downstream of this SNP interacted with an SNP in KCNK9 in Chr14b2, and two other SNPs interacted with an SNP downstream of the FAM135B gene in Chr14b2.All the other Chr03 SNPs (96 total) interacted with Chr14a (Figure 5a, Table S1).Although nearly the entire Chr03 interacted with Chr14a, the small region interacting with Chr14b1 had the most significant effect of Chr03.Other than the Chr03 region near LMX1A, the most significant effect of Chr14b1 was with Chr21 and Chr23 (Figure 5b,c; Table 4), noting that Chr23 almost interacted with Chr14a only, except for the inter-chromosome A×A effects with Chr14b1 (Figure 5c).effects of Chr14b2, except one, were in or near PTK2, TRAPPC9, KCNK9 and GPR20 (Table S1).Based on the limited gene information in Chr14b1, two alternative hypotheses for the inter-chromosome A×A effects of the six SNPs in the Chr14b1 region could be made: (1) the noncoding sequences in the Chr14b1 region had biological functions in the form of interactions with other chromosomes for FPC, and (2) any or all LOC112449593, LOC112449592 and TRNAC-GCA genes were responsible for the interactions between the six SNPs in the Chr14b1 region and the seventeen chromosomes due to LD with the significant SNPs.Hypothesis (1) should be the most likely reason for the interactions involving the six SNPs and implies major biological functions of the noncoding regions in the Chr14b1 region in the form of inter-chromosome A×A effects for FPC.Hypothesis (2) implies linked effects of the six SNPs through LD with any or all LOC112449593, LOC112449592 and TRNAC-GCA genes.

Inter-Chromosome A×A Effects of Chr20 and Chr05 Interacting with Chr14
Chr20 and Chr05, along with Chr14 and Chr06, also had highly significant additive effects for milk production traits.Therefore, it was of interest to determine whether the Chr20 and Chr05 regions affecting milk production traits also interacted with the Chr14 region for FPC.
The inter-chromosome A×A effects of Chr20 covered a large distance of 63.52 Mb (6.58-70.10Mb).Chr14a interacted with the 6.58-28.8Mb region (mostly the 20-28 Mb region) and Chr14b with the 30.61-42.14Mb region, whereas both Chr14a and Chr14b interacted with the remaining regions of Chr20 (Figure 6a).The most significant interchromosome A×A effect of Chr20 (log 10 (1/p) = 57.12)was that between rs136653182 about 332 Kb downstream of the ITGA1 gene of Chr20 and rs109208977 in LOC789384 of Chr14a (Figure 6a, Table 6).The 20-28 Mb region of Chr20 interacting with Chr14a was near the location with the most significant effects of Chr20 at 31-33 Mb.In contrast, the 30.61-42.14Mb region of Chr20 interacting with Chr14b had the most significant effects for milk yield among Chr20 SNPs.In particular, the NNT gene had highly significant effects for milk yield and had an SNP interacting with two SNPs in the Chr14b1 region that had the largest cluster of inter-chromosome A×A effects of Chr14b (Figure 4b, Table 3).The most significant A×A effect in the 30.61-42.14Mb region of Chr20 was that between rs133536911 about 10.59 Kb downstream of the FGF10 gene of Chr20 and rs134537992 in Chr14b1, and rs133536911 also interacted with four other SNPs of Chr14b (Table S1).These results showed that the inter-chromosome A×A effects between Chr20 and the Chr14 region involved the Chr20 region with highly significant effects for milk yield.about 701 Kb downstream of the MGST1-SLC15A5 region had inter-chromosome A×A effects with an SNP in Chr03 and an SNP in Chr20, and these two inter-chromosome A×A effects were the only ones not involving Chr14a or Chr14b (Table S1), noting that EPS8 had highly significant additive effects for FPC.The 23-44 Mb region had significant SNP additive effects for milk and fat yields, and this large region interacted with both Chr14a and Chr14b for FPC (Table S1).The inter-chromosome A×A effects of Chr05 were mostly in the 0.5-10.8Mb region interacting with Chr14a (Figure 4c, Table 3).The most significant SNP of chr05 was an intergenic SNP (rs109208465) about 71.13 Kb upstream of the BBS10 gene that interacted with six SNPs in PFCT (log 10 (1/p) = 59.25-62.89)and four SNPs in CASD and SMPD5 (log 10 (1/p) = 32.40-34.56) of Chr14a.However, the 0.5-10.8Mb Chr05 region did not have highly significant effects for milk and fat yields or FPC [1].The MGST1-SLC15A5 region (93.51-93.63Mb) of Chr05 had highly significant additive effects on fat yield and FPC but this region did not interact with Chr14a or Chr14b for FPC.The SNP closest to the MGST1-SLC15A5 region was rs134855280 at 92.59 Mb, which had a significant inter-chromosome A×A effect with an SNP in LOC789384 of Chr14a.It was interesting that an SNP in EPS8 about 701 Kb downstream of the MGST1-SLC15A5 region had inter-chromosome A×A effects with an SNP in Chr03 and an SNP in Chr20, and these two inter-chromosome A×A effects were the only ones not involving Chr14a or Chr14b (Table S1), noting that EPS8 had highly significant additive effects for FPC.The 23-44 Mb region had significant SNP additive effects for milk and fat yields, and this large region interacted with both Chr14a and Chr14b for FPC (Table S1).

Patterns of A×A Epistasis Effects
The A×A values of the four allelic combinations (AB, Ab, aB, ab) of each pair of loci typically had large absolute values for the most positive and negative allelic combinations.Let AC1, AC2, AC3 and AC4 represent the four allelic combinations from the most positive combination to most negative combination and let aa1-aa4 represent the A×A values of AC1-AC4, where 'AC' stands for 'allelic combination'.Then, AC1 and AC4 had the largest absolute A×A values, whereas AC2 and AC3 had considerably smaller absolute A×A values than those of AC1 and AC4 (Table S1, Tables 7 and 8).Consequently, the size of the A×A effect of two loci as a contrast of aa1-aa4 (Equations ( 1) and ( 2) in Materials and Methods) was mostly determined by the A×A values of AC1 and AC4.Therefore, the discussion of A×A patterns focused on the A×A values of AC1 and AC4, which had two patterns: (1) the two A×A values involved the same chr14 allele and two non-Chr14 alleles such as the 1_1 allelic combination for AC1 and 1_2 allelic combination for AC4, and (2) the A×A values involved two chr14 alleles and the same non-Chr14 allele, such as 1_1 for AC1 and 2_1 for AC4.Most Chr14a A×A values (1554 out of 2148, or 72%) had pattern (1), whereas most Chr14b A×A values (346 out of 613, or 56%) had pattern (2).It was interesting that no AC1 and AC4 of any SNP pair involved completely different alleles, such as 1_1 for AC1 and 2_2 for AC4, or 1_2 for AC1 and 2_1 for AC4.Table 7 shows examples of the Chr14a A×A values where the AC1 and AC4 of each SNP pair had the same Chr14a allele and two different non-Chr14 alleles.SNP rs109421300 of DGAT1 should be a highly recognizable SNP because this SNP had the most significant effects for all five production traits: milk, fat and protein yields and fat and protein percentages.Allele 1 of rs109421300 had an extreme antagonism between fat yield and milk and protein yields, with the most positive effect for fat yield and FPC and most negative effects for milk and protein yields [1].In this study, allele 1 of rs109421300 was the common Chr14a allele of AC1 and AC4 for all nine pairs of A×A values, and each of the nine SNPs in the SGN region of Chr06 had both alleles in AC1 and AC4.The AC1 and AC4 for seven of the nine SNP pairs had similar absolute values, indicating that these AC1 and AC4 were approximately symmetric.The combination of allele 1 of rs109421300 with a Chr06 allele was positive (aa1, Table 7), whereas the combination of allele 1 of rs109421300 with the alternative Chr06 allele of each Chr06 SNP was negative (aa4, Table 7).The other A×A values of Chr14a had similar patterns.The results of Table 7 indicate that one allele of a Chr14a SNP interacted with both alleles of a non-Chr14 SNP for most of the significant SNP pairs involving Chr14a.Table 8 shows examples of the Chr14b A×A values where the AC1 and AC4 of each SNP pair had the same non-Chr14 allele and two different Chr14b alleles.Allele 1 of SNP rs42368654 downstream of LMX1A of Chr03 was the common allele of AC1 and AC4 with five SNPs of Chr14b1.The size of AC1 was 2-5 times as large as that of AC4 for the five A×A values, indicating that the interaction between allele 1 of rs42368654 of Chr03 and one Chr14b1 was the main contributor of the five A×A values.Of the twenty SNP pairs in Table 8, the size of AC1 was larger than that of AC4 for eighteen pairs, AC4 was larger than AC1 for two pairs, and AC1 and AC2 had approximately the same sizes for two pairs.

Holstein Population and SNP Data
The Holstein population in this study had 1,231,898 first lactation cows with phenotypic observations of fat percentage (FPC) and genotypes of 78,964 original and imputed SNPs.The SNP genotypes were from 32 SNP chips with various densities and were imputed to 78,964 SNPs via the FindHap algorithm [19] as a routine procedure for genomic evaluation by the Council on Dairy Cattle Breeding (CDCB) [20].The phenotypic values used in the GWAS analysis were the phenotypic residuals after removing fixed non-genetic effects available from the December 2022 U.S. Holstein genomic evaluation by the CDCB.Basic statistics of the cows and phenotypic data of FPC are given in Table S2.With the requirement of a 0.05 minor allele frequency, the number of SNPs for the GWAS analysis was 75,198 SNPs.A strict criterion of log 10 (1/p) > 32 was used to declare the statistical significance of any inter-chromosome epistasis effect.This requirement ensured that any significant effect had better statistical significance than that of the highest statistical significance of log 10 (1/p) = 32 shown in Figure 2a,b.The log 10 (1/p) > 32 requirement was stricter than the requirement of log 10 (1/p) > 12 for the Bonferroni correction with 0.05 genomewide false positives.The SNP and gene positions were those from the ARS-UCD1.3cattle genome assembly [21].Genes containing or in the proximity of highly significant effects were identified as candidate genes affecting FPC.

GWAS Analysis
The A×A value of each of the four allelic combinations (AB, Ab, aB, ab) of two loci was calculated as the deviation of the mean of the allelic combination from the population mean and the additive values of the two alleles in the allelic combination [22,23]: where The A×A effect of two loci was calculated as a contrast of the four A×A values and this contrast was further expressed as the A×A contrast of the nine genotypic values for epistasis testing [24]: where αα = A×A effect of the two loci as a contrast of the four A×A values of the four allelic combinations of the two loci; ( aa) AB , (aa) Ab , (aa) aB and (aa) ab are the four A×A values of the four allelic combinations of AB, Ab, aB and ab, respectively, defined by Equation (1); µ AB = the mean of genotypic values with the AB allelic combination, µ Ab = the mean of genotypic values with the Ab allelic combination, µ aB = the mean of genotypic values with the aB allelic combination, µ ab = the mean of genotypic values with the ab allelic combination; g = column vector of the nine SNP genotypic values of the two loci: g AABB , g AABb , g AAbb , g AaBB , g AaBb , g Aabb , g aaBB , g aaBb , g aabb ; s a×a = row vector of the A×A contrast coefficients of nine SNP genotypic values; and L a×a = A×A effect of the two loci as a contrast of the nine SNP genotypic values.In the absence of allelic interactions between the two loci, the A×A effect of Equation ( 2) is expected to be null because each A×A value is expected to be null.In the presence of an allele × allele interaction between the two loci, the [( aa) AB − (aa) Ab ] − [(aa) aB − (aa) ab ] definition of the A×A effect indicates that the allelic difference of locus 2 changes in the presence of the two different alleles of locus 1, whereas the [(aa) AB − (aa) aB ] − [(aa) Ab − (aa) ab ] definition of A×A effect indicates that the allelic difference of locus 1 changes in the presence of the two different alleles of locus 2. Therefore, a significant A×A effect expressed as L a×a = s a×a g indicates the presence of an allele × allele interaction between the two loci due to the equivalence between this expression and any of the other four expressions in Equation ( 2).The GWAS analysis of A×A effects used an approximate generalized least squares (AGLS) method.The AGLS method combines the least squares (LS) tests implemented by EPISNP1mpi [25,26] with the estimated breeding values from routine genetic evaluation using the entire U.S. Holstein population.The statistical model was y = µI + X g g + Za + e = Xb + Za + e where y = column vector of phenotypic deviation after removing fixed nongenetic effects such as heard-year-season (termed as 'yield deviation' for any trait) using a standard procedure for the CDCB/USDA genetic and genomic evaluation; µ = common mean; I = identity matrix; g = column vector of genotypic values; X g = model matrix of g; b = (µ, g′)′, X = (I, X g ); a = column vector of additive polygenic values; Z = model matrix of a; and e = column vector of random residuals.The first and second moments of Equation ( 3) are E(y) = Xb and var(y) = V = ZGZ′ + R = σ 2 a ZAZ′ + σ 2 e I, where σ 2 a = additive variance, A = additive relationship matrix and σ 2 e = residual variance.The problem of estimating the b vector that includes SNP genotypic values in Equation ( 3) is the requirement of inverting the V if the generalized least squares (GLS) method is used, or inverting the A matrix and the coefficient matrix of the mixed model equations (MME) if the MME method is used [27].However, both V and MME could not be inverted for our sample size.To avoid inverting these large matrices, the GWAS used the method of approximate GLS (AGLS), which replaces the polygenic additive values (a) with the best linear unbiased prediction based on pedigree relationships [1].The AGLS method is based on the following results: where y * = y − Z â and â is the best linear unbiased prediction (BLUP) of a. Equation ( 4) is the GLS solution, and Equation ( 5) is the MME solution of b.These two equations yield identical results, and b from either equation is termed the best linear unbiased estimator (BLUE) [27].If â is known, the LS version of BLUE given by Equation ( 5) is computationally efficient relative to the GLS of Equation ( 4 where y * = y − Z ã, and ã is the column vector of 2(PTA), with PTA being the predicted transmission ability from the routine genetic evaluation.Equation ( 6) achieves the benefit of sample stratification correction from mixed models using pedigree relationships without the computing difficulty of inverting V or A. The second approximation of the AGLS approach is the t-test using the LS rather than the GLS formula of the t-statistic, to avoid using the V inverse in the GLS formula.The significance tests for A×A SNP effects used the t-tests of the A×A contrast of the estimated two-locus SNP genotypic values [24,25]:

Conclusions
This GWAS using over 1.2 million Holstein cows confirmed that a Chr14 region about 9.38 Mb region in size had significant inter-chromosome additive × additive (A×A) effects

Figure 1 .
Figure 1.Significant inter-chromosome A×A effects of FPC (log10(1/p) > 32).(a) The 0.14-9.52Mb region of Chr14 (marked by the red arrow) had significant inter-chromosome A×A effects with all chromosomes.(b) Significant inter-chromosome A×A effects of each chromosome.

Figure 1 .
Figure 1.Significant inter-chromosome A×A effects of FPC (log 10 (1/p) > 32).(a) The 0.14-9.52Mb region of Chr14 (marked by the red arrow) had significant inter-chromosome A×A effects with all chromosomes.(b) Significant inter-chromosome A×A effects of each chromosome.

2. 1 .
Structure of the 9.38 Mb Chr14 Region Interacting with All Chromosomes for FPC

Figure 2 .
Figure 2. Inter-chromosome A×A effects of Chr14.(a) Significant inter-chromosome A×A effects of Chr14 among the top 50,000 pairs of inter-chromosome A×A effects.Effects with log10(1/p) > 32 were considered statistically significant.(b) Inter-chromosome A×A effects of the 0.14-9.52Mb region with two sub-regions of Chr14a and Chr14b that were separated by a 1.3 Mb gap region.Effects with log10(1/p) > 32 were considered statistically significant.(c) Inter-chromosome A×A effects of Chr14a with log10(1/p) > 32.(d) Inter-chromosome A×A effects of Chr14b with log10(1/p) > 32.

Figure 2 .
Figure 2. Inter-chromosome A×A effects of Chr14.(a) Significant inter-chromosome A×A effects of Chr14 among the top 50,000 pairs of inter-chromosome A×A effects.Effects with log 10 (1/p) > 32 were considered statistically significant.(b) Inter-chromosome A×A effects of the 0.14-9.52Mb region with two sub-regions of Chr14a and Chr14b that were separated by a 1.3 Mb gap region.Effects with log 10 (1/p) > 32 were considered statistically significant.(c) Inter-chromosome A×A effects of Chr14a with log 10 (1/p) > 32.(d) Inter-chromosome A×A effects of Chr14b with log 10 (1/p) > 32.

Figure 3 .
Figure 3. Significant inter-chromosome A×A effects of FPC between Chr14a and Chr06 with log 10 (1/p) > 32.(a) Inter-chromosome A×A effects of Chr06.(b) Inter-chromosome A×A effects between Chr14a and the SGN region Chr06.(c) Manhattan plot of statistical significance of interchromosome A×A effects of Chr06.(d) Manhattan plot of statistical significance of inter-chromosome A×A effects between Chr14a and the SGN region Chr06.

Figure 4 .
Figure 4. Examples of significant inter-chromosome A×A effects of Chr14 for FPC with log10(1/p) > 32.(a) Inter-chromosome A×A effects of Chr02.(b) Inter-chromosome A×A effects of Chr29.(c) Interchromosome A×A effects of Chr17.(d) Inter-chromosome A×A effects of ChrX.

Figure 4 .
Figure 4. Examples of significant inter-chromosome A×A effects of Chr14 for FPC with log 10 (1/p) > 32.(a) Inter-chromosome A×A effects of Chr02.(b) Inter-chromosome A×A effects of Chr29.(c) Inter-chromosome A×A effects of Chr17.(d) Inter-chromosome A×A effects of ChrX.

Figure 5 .
Figure 5. Inter-chromosome A×A effects of Chr03 and Chr10 for FPC.(a) Inter-chromosome A×A effects between Chr03 and Chr14.(b) Inter-chromosome A×A effects between Chr21 and Chr14.(c) Inter-chromosome A×A effects between Chr23 and Chr14.(d) Inter-chromosome A×A effects between Chr10 and Chr14.

Figure 5 .
Figure 5. Inter-chromosome A×A effects of Chr03 and Chr10 for FPC.(a) Inter-chromosome A×A effects between Chr03 and Chr14.(b) Inter-chromosome A×A effects between Chr21 and Chr14.(c) Inter-chromosome A×A effects between Chr23 and Chr14.(d) Inter-chromosome A×A effects between Chr10 and Chr14.

Figure 6 .
Figure 6.Inter-chromosome A×A effects of Chr20 and Chr05 for FPC.(a) Inter-chromosome A×A effects between Chr20 and Chr14.(b) Manhattan plot of statistical significance of inter-chromosome A×A effects of Chr20.(c) Inter-chromosome A×A effects between Chr05 and Chr14.(d) Manhattan plot of statistical significance of inter-chromosome A×A effects of Chr05.
aa) ik = A×A value of allelic combination of the i th allele of locus 1 (A or a allele) and the k th allele of locus 2 (B or b allele), µ ik = the mean of the genotypic values with allelic combination of the i th allele of locus 1 (A or a allele) and the k th allele of locus 2 (B or b allele), µ = the population mean of genotypic values, a i = µ i − µ (i = A or a) = additive value of i th allele of locus 1 (A or a allele), a k = µ k − µ (k = B or b) = additive value of the k th allele of locus 2 (B or b allele), µ i = the mean of genotypic values with the i th allele of locus 1 (A or a allele), and µ k = the mean of genotypic values with the k th allele of locus 2 (B or b allele).
b = (X′V −1 X) − X′V −1 y (4) b = (X′R −1 X) − (X′R −1 y − X′R −1 Z â) = (X′X) − X′(y − Z â) = (X′X) − X′y * ), requiring the V inverse, or the joint MME solutions of b and â, requiring the inverse of the coefficient matrix of the MME.The AGLS method uses two approximations.The first approximation is to use ã from routine genetic evaluation as an approximation of â in Equation (5): b = (X′X) − X′(y − Z ã) = (X′X) − X′y * = A×A contrast of the nine genotypic values defined by Equation (1); var(L a×a = standard deviation of L a×a ; s a×a = row vector of A×A contrast coefficients; v 2 = (y − X b) ′ (y − X b)/(n − k) = estimated residual variance; ĝ = column vector of the AGLS estimates of the nine SNP genotypic values of the two loci; and (X ′ X) − gg = submatrix of (X ′ X) − corresponding to ĝ.

Table 1 .
Inter-chromosome A×A effects between four SNPs in the PFCT region of Chr14a and six SNPs in SLC4A4-NPFFR2 of Chr06.

Table 2 .
Inter-chromosome A×A effects between four SNPs in CASDS region of Chr14a and the SGN region of Chr06.

Table 3 .
Top inter-chromosome A×A effects of the Chr14 excluding those with the SGN region of Chr06.
Pos is the chromosome position of the SNP.

Table 3 .
Top inter-chromosome A×A effects of the Chr14 excluding those with the SGN region of Chr06.

Table 4 .
Top 20 inter-chromosome A×A effects of the Chr14b1 region.

Table 4 .
Top 20inter-chromosome A×A effects of the Chr14b1 region.
Pos is the chromosome position of the SNP.Chr31 is the nonrecombining region of ChrX.

Table 5 .
Top 20inter-chromosome A×A effects of the Chr14b2 region.

Table 6 .
Top inter-chromosome A×A effects of the Chr05 and Chr20.

Table 7 .
Patterns of A×A epistasis effects between Chr14a and the SGN region of Chr06.

Table 8 .
Patterns of A×A epistasis effects of Chr14b.AC4 are the four allelic combinations of the two loci, and aa1-aa4 are the A×A epistasis values of AC1-AC4.