{"id":193938,"date":"2025-10-01T04:46:07","date_gmt":"2025-10-01T04:46:07","guid":{"rendered":"https:\/\/www.newsbeep.com\/us\/193938\/"},"modified":"2025-10-01T04:46:07","modified_gmt":"2025-10-01T04:46:07","slug":"mapping-genes-for-resilient-dairy-cows-by-means-of-across-breed-genome-wide-association-analysis-bmc-genomics","status":"publish","type":"post","link":"https:\/\/www.newsbeep.com\/us\/193938\/","title":{"rendered":"Mapping genes for resilient dairy cows by means of across-breed genome-wide association analysis | BMC Genomics"},"content":{"rendered":"

Data processing and evaluation as well as statistical analyses were performed in software R 4.2.1 [16<\/a>]. Genotype preparation, phasing, imputation and population analysis were done by using PLINK 1.9 and PLINK 2.0 [17<\/a>] and Beagle 5.4 [18<\/a>, 19<\/a>]. ASReml-R 4.1.0 [20<\/a>] was used for univariate estimation of breeding values and GCTA 1.94.1 [21<\/a>] for GWAS with de-regressed GEBVs as phenotypes.<\/p>\n

Materials<\/p>\n

The dataset consists of 5,703 dairy cows of the breeds German Holstein (HF, \\(n=\\) 2,300), German Fleckvieh (FV, \\(n=\\) 2,330) and German Brown Swiss (BS, \\(n=\\) 1,073), originating from 48 farms. The farms are members of the breeding organization Rinderunion Baden-W\u00fcrttemberg e.V. (Herbertingen, Germany) and work with automatic milking systems. The milking data, consisting of yield per milking and start time of milking, from the automatic milking systems were longitudinally recorded between October 2017 and April 2024 and were provided by the Landeskontrollverband Baden-W\u00fcrttemberg e.V. (Stuttgart, Germany).<\/p>\n

Animals in our study were genotyped with 50\u00a0K SNP chips. For HF, the Eurogenomics-MD-Chip was used in several versions. The merging of the chips and the imputation of undetermined SNPs were done by Vereinigte Informationssysteme Tierhaltung w.V. (Verden, Germany). Genotyping of FV and BS was done by Landeskontrollverband Bayern (Munich, Germany) and Bayerische Landesanstalt f\u00fcr Landwirtschaft (Grub, Germany) using Illumina BovineSNP50v2 BeadChip or DAC-BS50 Illumina. For both breeds, we obtained a standardized data set over the different chips, whereby missing SNPs were called with the population mean. We excluded SNPs on sex chromosomes, on undefined position, with minor allele frequency (MAF)\u2009<\u20090.03 or missing call rate\u2009>\u20090.01. The MAF threshold was stricter than in comparable studies [13<\/a>, 22<\/a>] in order to cause fewer false-positive results in the GWAS and was based on a study with a similar data set to our BS sample [23<\/a>]. SNPs that are not included in the high-density reference panels were excluded, as they would be automatically removed during imputation. All genotypes were updated to the reference genome assembly ARS-UCD1.2 by manually replacing the existing coordinates with the updated ones [24<\/a>]. SNPs, which are not listed in the ARS-UCD1.2, were excluded. This left 31,416 SNPs for HF, 34,238 SNPs for FV and 31,527 SNPs for BS.<\/p>\n

A detailed overview of the size of the raw genotype and phenotype data and remaining data sets after filtering per breed can be found in the Supplementary (see Additional File 1).<\/p>\n

MethodsCalculation of resilience indicator traits<\/p>\n

Milk yields per milking were added up to daily milk yields, whereby the first milking was divided proportionally between the current and the previous day. Daily milk yields were filtered as follows. Days where the number of animals milked on the farm deviated by three standard deviations from the mean number of animals normally milked on the farm as well as the day before and after were removed for all animals on this farm. Furthermore, the day before and after an individual gap of an animal, meaning one or more days without available data, was removed. Lactations were considered from day 10 to day 305, and those with less than 50% coverage respectively 148 data points were excluded. For each lactation of each cow, a lactation curve was predicted using spline interpolation with pspline in R [25<\/a>] with degrees of freedom set to 5 and data points weighted by the strength of assumed disturbances. Details are described in Ke\u00dfler et al. [6<\/a>].<\/p>\n

We considered three resilience indicator traits. The variance of the deviations between observed and predicted absolute daily milk yield (\\({v}_{d}\\)) was calculated as:<\/p>\n

$${v}_{d}\\left({\\mathbf{y}}_{il}\\right)=\\frac{1}{n-1}*\\sum_{t=10}^{n}{\\left(\\left({{\\varvec{y}}}_{{\\varvec{i}}{\\varvec{l}}{\\varvec{t}}}-{\\widehat{{\\varvec{y}}}}_{{\\varvec{i}}{\\varvec{l}}{\\varvec{t}}}\\right)-\\overline{\\left({{\\varvec{y}}}_{{\\varvec{i}}{\\varvec{l}}}-{\\widehat{{\\varvec{y}}}}_{{\\varvec{i}}{\\varvec{l}}}\\right)}\\right)}^{2}$$<\/p>\n

with \\({\\varvec{y}}\\) a \\(n\\)-vector with daily milk yields of a cow \\(i\\) in a certain lactation \\(l\\), \\(\\widehat{\\mathbf{y}}\\) the vector with interpolated daily milk yields and \\(t\\) the lactation day. The variance of relative daily milk yields (\\({v}_{r}\\)) in considered lactation period was calculated as:<\/p>\n

$${v}_{r}\\left({\\mathbf{y}}_{il}\\right)=\\frac{1}{n-1}*\\sum_{t=10}^{n}{\\left(\\left({\\alpha }_{ilt} {\\mathbf{y}}_{ilt}\\right)-\\overline{\\left({\\alpha }_{il} {\\mathbf{y}}_{il}\\right)}\\right)}^{2}$$<\/p>\n

with \\({\\varvec{y}}\\) a \\(n\\)-vector with daily milk yields of a cow \\(i\\) in a certain lactation \\(l\\) and the scaling factor \\({\\alpha }_{il}=\\frac{100}{\\sum_{t=10}^{n}{y}_{ilt}}\\) with \\(t\\) the lactation day. The two parameters were logarithmized to conform to the assumption of normal distribution according to previous analyses in Ke\u00dfler et al. [6<\/a>].<\/p>\n

The autocorrelation of the deviations between observed and predicted absolute daily milk yield (\\({r}_{Auto}\\)) was calculated as:<\/p>\n

$${r}_{Auto}\\left({\\mathbf{y}}_{il}\\right)=\\frac{cov\\left({\\mathbf{y}}_{il}-{\\widehat{\\mathbf{y}}}_{il}\\right)}{var\\left({\\mathbf{y}}_{il}-{\\widehat{\\mathbf{y}}}_{il}\\right)}$$<\/p>\n

and in detail:<\/p>\n

$${r}_{Auto}\\left({\\mathbf{y}}_{il}\\right)=\\frac{\\frac{1}{n-2}*\\sum_{t=2}^{n}\\left(\\left(\\left({{\\varvec{y}}}_{{\\varvec{i}}{\\varvec{l}}{\\varvec{t}}}-{\\widehat{{\\varvec{y}}}}_{{\\varvec{i}}{\\varvec{l}}{\\varvec{t}}}\\right)-\\stackrel{-}{\\left({{\\varvec{y}}}_{{\\varvec{i}}{\\varvec{l}}}-{\\widehat{{\\varvec{y}}}}_{{\\varvec{i}}{\\varvec{l}}}\\right)}\\right)*\\left(\\left({{\\varvec{y}}}_{{\\varvec{i}}{\\varvec{l}},{\\varvec{t}}-1}-{\\widehat{{\\varvec{y}}}}_{{\\varvec{i}}{\\varvec{l}}{\\varvec{t}}}\\right)-\\stackrel{-}{\\left({{\\varvec{y}}}_{{\\varvec{i}}{\\varvec{l}}}-{\\widehat{{\\varvec{y}}}}_{{\\varvec{i}}{\\varvec{l}}}\\right)}\\right)\\right)}{\\frac{1}{n-1}*\\sum_{t=1}^{n}{\\left(\\left({{\\varvec{y}}}_{{\\varvec{i}}{\\varvec{l}}{\\varvec{t}}}-{\\widehat{{\\varvec{y}}}}_{{\\varvec{i}}{\\varvec{l}}{\\varvec{t}}}\\right)-\\stackrel{-}{\\left({{\\varvec{y}}}_{{\\varvec{i}}{\\varvec{l}}}-{\\widehat{{\\varvec{y}}}}_{{\\varvec{i}}{\\varvec{l}}}\\right)}\\right)}^{2}}$$<\/p>\n

Genetic analysis<\/p>\n

The resilience indicators described above were genetically analyzed separately for each breed and two different lactation groups. The first one, labelled \\(P\\) for primiparous, included phenotypes from the first lactation. For the second one, named \\(M\\) for multiple lactations, we added information also from higher lactations of the cows in \\(P\\) as repeated measurements. Thus, a univariate analysis was performed for each of six groups: first lactation HF, all lactation HF, first lactation FV, all lactation FV, first lactation BS, all lactation BS. The animal model for the data sets with information about cows in the first lactation was:<\/p>\n

with \\(y\\) vector of observations, \\(u\\) vector of additive genetic effects, \\(e\\) vector of residuals, and incidence matrices \\(X\\) and \\(Z\\). The animal model for all lactations of the considered cows was:<\/p>\n

where \\(pe\\) is the vector of permanent environmental effects and \\(W\\) is the incidence matrix. The additive genetic effects \\(u\\), the permanent environmental effects \\(pe\\) and the residual effects \\(e\\) were normally distributed with \\(u \\sim N(0,G{\\sigma }_{u}^{2})\\), \\(pe \\sim N\\left(0,{I}_{pe}{\\sigma }_{pe}^{2}\\right)\\) and \\(e \\sim N\\left(0,{I}_{e}{\\sigma }_{e}^{2}\\right)\\), whereby \\(I\\) are identity matrices. The genomic relationship matrix \\(G\\) was built according to VanRaden et al. [26<\/a>]. The vector of fixed effects \\(b\\) comprises in each model age at first calving in month, herd-year-season effect and completeness of lactation data. In the second model the effect of the lactation separated in first and higher lactation was added. Season in herd-year-season was divided into spring from March to May, summer from June to August, autumn from September to November, and winter from December to February and completeness of lactation data was graded in\u2009\u2265\u200990%, 80%\u2014<\u200990%, 70%\u2014<\u200980%, 60%\u2014<\u200970%, 50%\u2014<\u200960%. In all analyses, effect levels with fewer than six animals are excluded. This left 1,370 HF cows for \\(P\\) and 1,538 HF cows for \\(M\\), for FV 1,135 cows in \\(P\\) and 1,551 in \\(M\\) and for BS \\(P\\) includes data from 310 and \\(M\\) from 488 cows. For further analysis, we de-regressed the cow\u2019s additive genetic effects, i.e. their GEBV, by dividing the effects by their reliabilities and standardized them to a mean of 0 and a standard deviation of 1 within each group. We also calculated the heritability as the proportion of additive genetic variance of the total variance.<\/p>\n

Additionally, a selection index (\\(SI\\)) was computed that combines the resilience indicators into a single trait. The selection index is a linear combination of the GEBVs of individual resilience indicators, with weights chosen to optimally describe overall resilience. Weighting factors were determined in Ke\u00dfler et al. [8<\/a>], by maximizing the selection response in health traits when selecting on resilience indicators based on the performance trait daily milk yield. The weighting factors were adopted from this study, which were estimated in the same data set [8<\/a>]. The selection index was calculated as<\/p>\n

$$SI\\left({\\mathbf{y}}_{i}\\right)=0.65*{\\widetilde{v}}_{d}\\left({\\mathbf{y}}_{i}\\right)+ 0.35* {\\widetilde{v}}_{r}\\left({\\mathbf{y}}_{i}\\right),$$<\/p>\n

where \\({\\widetilde{v}}_{d}\\left({\\mathbf{y}}_{i}\\right)\\) and \\({\\widetilde{v}}_{r}\\left({\\mathbf{y}}_{i}\\right)\\) were GEBVs of the logarithmized, negated \\({v}_{d}\\left({\\mathbf{y}}_{i}\\right)\\) and \\({v}_{r}\\left({\\mathbf{y}}_{i}\\right)\\), for cow \\(i\\). The selection index was calculated in both lactation groups, \\(P\\) and \\(M\\), and all three breeds.<\/p>\n

Imputation<\/p>\n

For further analyses we imputed the 50\u00a0K genotypes to high-density genotypes. The following reference panels were used. For HF, we used data from 1,278 Australian Holsteins with 585,517 SNPs that were genotyped with the BovineHD Genotyping BeadChip [27<\/a>]. High-density data from 3,439 FV cows and bulls including 629,028 SNPs genotyped with Illumina BovineHD Bead chip [28<\/a>] and 192 BS individuals with 613,140 SNPs genotyped with Illumina Bovine-HD [23<\/a>] were available. Similar to the 50\u00a0K genotype data, SNPs whose MAF was\u2009<\u20090.03, whose missing call rate was\u2009>\u20090.01 or which were found on sex chromosomes were removed. The positions of considered SNPs were updated to the reference genome ARS-UCD1.2 [24<\/a>]. 522,726 SNPs in HF, 580,873 SNPs in FV and 503,971 SNPs in BS remained.<\/p>\n

Reference panel data were phased using beagle 5.4 within each dataset. 50\u00a0K and high-density data of each breed were concorded with conform-gt and then imputed with beagle 5.4. According to [29<\/a>], the effective population size Ne was set to 100. Apart from that, the default settings were used. SNPs with a dosage R-Squared\u2009<\u20090.75 are removed to maintain a sufficiently high imputation quality [19<\/a>]. The final imputed data set contains 2,300 HF cows with 503,263 SNPs, 2,330 FV cows with 569,056 SNPs and 1,073 BS cows with 493,637 SNPs. For the across-breed population analysis, the imputed genotypes were merged with PLINK 1.9 [30<\/a>], leaving 401,491 SNPs shared across all breeds.<\/p>\n

Population analysis<\/p>\n

Population analyses were done in PLINK 1.9 for MDS and PLINK 2.0 for LD analysis [17<\/a>]. Genomic relationship between the breeds was studied with a MDS based on imputed genotypes, considering only SNPs that remained in all three breeds after all filtering steps described above. The LD was analysed across breeds and breed-specific, using the maximum number of available SNPs in the imputed genotype datasets of each breed, respectively the merged imputed dataset after all filtering steps. The LD was calculated as the \\({r}^{2}\\) [31<\/a>]:<\/p>\n

$${r}^{2}= \\frac{{\\left({p}_{{A}_{1}{B}_{1}}-{p}_{{A}_{1}}{p}_{{B}_{1}}\\right)}^{2}}{{p}_{{A}_{1}}{p}_{{A}_{2}}{p}_{{B}_{1}}{p}_{{B}_{2}}}$$<\/p>\n

with \\(p\\) the allele frequency, \\({A}_{1}\\) and \\({A}_{2}\\) alleles on gene locus \\(A\\) and \\({B}_{1}\\) and \\({B}_{2}\\) alleles on gene locus \\(B\\). The LD between pairs of SNPs that were a maximum of 1,000 SNPs or 3,000 kilobases apart was calculated.<\/p>\n

Genome-wide association study<\/p>\n

For each of the six groups defined above, a GWAS was performed for the resilience indicator traits and the selection index. We optimized the GWAS model by minimizing the difference of the genomic inflation factor lambda \\({\\lambda }_{GC}\\) to 1 in each breed. \\({\\lambda }_{GC}\\) is calculated as the ratio between the median of the observed chi-squared test statistics and the expected median of 0.455 [32<\/a>]. The value should be close to 1 or slightly above [22<\/a>, 33<\/a>]. We therefore decided for a mixed linear model based association analysis (\\(MLMA\\)) for HF, a \\(MLMA\\) excluding the GRM of the chromosome where the SNP is located (\\(MLM-LOCO\\)) for BS and a \\(MLM-LOCO\\) corrected for a fixed effect of the principal components calculated based on the whole genome (\\({\\text{MLM}-\\text{LOCO}}_{PC Genome}\\)) to account for population stratification for FV [21<\/a>] (\\({\\lambda }_{GC}\\) of the selected GWAS are provided in the results, all \\({\\lambda }_{GC}\\) are provided in Additional File 2).<\/p>\n

GWAS were performed using the following model:<\/p>\n

$$y = \\mu + Xb+ Zu + e,$$<\/p>\n

with \\(y\\) vector of de-regressed, standardized additive genetic effects of each animal derived from univariate analysis for each resilience indicator trait and the selection index, the population mean \\(\\upmu\\), fixed effects \\(b\\) with the incidence matrix \\(X\\), including the reference allele coded with 0, 1 or 2, the vector of additive genetic effects \\(u\\) with \\(u \\sim N(0,G{\\sigma }_{u}^{2})\\) and \\(G\\) the genomic relationship matrix created according to Yang et al. [34<\/a>], and the residual term \\(e\\) with \\(e \\sim N(0,{{I}_{e}\\sigma }_{e}^{2})\\). In \\(MLM-LOCO\\) and \\({\\text{MLM}-\\text{LOCO}}_{PC Genome}\\), \\(G\\) is calculated for each chromosome separately and leaves out the chromosome, where the SNP to be tested for association is located. The vector \\(b\\) includes the additive effect of the candidate SNP to be tested for association and for \\({\\text{MLM}-\\text{LOCO}}_{PC Genome}\\) in FV a covariate of 20 principal components, which were calculated based on the whole genome with software GCTA [35<\/a>]. For the calculation of GWAS in GCTA, the reference allele was specified in accordance with ARS-UCD1.2.<\/p>\n

For across-breed evaluation, we pooled the p-values in data set \\(P\\) and data set \\(M\\) per resilience indicator using the package poolr in R [36<\/a>]. We used the Fisher method, which is based on the chi2-statistic:<\/p>\n

$${X}^{2}=-2\\sum_{i=1}^{k}\\text{ln}({p}_{i})$$<\/p>\n

with K number of tests (three) and pi the p-values of the GWAS.<\/p>\n

The results of the single GWAS and pooled GWAS are shown as Manhattan plots, whereby two significance thresholds were defined. Genome-wide significance was determined using Bonferroni correction with, with \\(N\\) number of SNPs, and a less stringent, suggestive level of significance was set to \\(p\\le \\frac{0.05}{N}\\). For SNPs showing significant association at the stringent significance level, the Cattle QTL Database release 53 [27<\/a>] was used to search for QTL in the range from 25 kB up- and downstream.<\/p>\n","protected":false},"excerpt":{"rendered":"Data processing and evaluation as well as statistical analyses were performed in software R 4.2.1 [16]. Genotype preparation,…\n","protected":false},"author":2,"featured_media":193939,"comment_status":"","ping_status":"","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[50],"tags":[13114,112631,112634,257,200,37224,18666,112632,112633,5081,93117,19647,48141,6057,79,112630],"class_list":["post-193938","post","type-post","status-publish","format-standard","has-post-thumbnail","category-genetics","tag-animal-genetics-and-genomics","tag-autocorrelation","tag-dairy-cattle","tag-general","tag-genetics","tag-gwas","tag-life-sciences","tag-linkage-disequilibrium","tag-mds","tag-microarrays","tag-microbial-genetics-and-genomics","tag-plant-genetics-and-genomics","tag-proteomics","tag-resilience","tag-science","tag-variance"],"_links":{"self":[{"href":"https:\/\/www.newsbeep.com\/us\/wp-json\/wp\/v2\/posts\/193938","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/www.newsbeep.com\/us\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/www.newsbeep.com\/us\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/www.newsbeep.com\/us\/wp-json\/wp\/v2\/users\/2"}],"replies":[{"embeddable":true,"href":"https:\/\/www.newsbeep.com\/us\/wp-json\/wp\/v2\/comments?post=193938"}],"version-history":[{"count":0,"href":"https:\/\/www.newsbeep.com\/us\/wp-json\/wp\/v2\/posts\/193938\/revisions"}],"wp:featuredmedia":[{"embeddable":true,"href":"https:\/\/www.newsbeep.com\/us\/wp-json\/wp\/v2\/media\/193939"}],"wp:attachment":[{"href":"https:\/\/www.newsbeep.com\/us\/wp-json\/wp\/v2\/media?parent=193938"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/www.newsbeep.com\/us\/wp-json\/wp\/v2\/categories?post=193938"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/www.newsbeep.com\/us\/wp-json\/wp\/v2\/tags?post=193938"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}