- Email: [email protected]
Genome-wide association study of milk components in Chinese Holstein cows using single nucleotide polymorphism
Journal Pre-proof
Genome-wide association study of milk components in Chinese Holstein cows using single nucleotide polymorphism Tianzhen Wang , Jiao Li , Xue Gao , Wenqin Song , Chengbin Chen , Dawei Yao , Jing Ma , Lingyang Xu , Yi Ma PII: DOI: Reference:
S1871-1413(19)31077-7 https://doi.org/10.1016/j.livsci.2020.103951 LIVSCI 103951
To appear in:
Livestock Science
Received date: Revised date: Accepted date:
2 August 2019 19 January 2020 28 January 2020
Please cite this article as: Tianzhen Wang , Jiao Li , Xue Gao , Wenqin Song , Chengbin Chen , Dawei Yao , Jing Ma , Lingyang Xu , Yi Ma , Genome-wide association study of milk components in Chinese Holstein cows using single nucleotide polymorphism, Livestock Science (2020), doi: https://doi.org/10.1016/j.livsci.2020.103951
This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2020 Published by Elsevier B.V.
Highlights
A total of 18 SNPs located in BTA14 were detected which were associated with at least one milk traits using mixed linear model after Bonferroni correction.
One 2.5 Mb region (range from 1.5 to 4.0 Mb) with high linkage disequilibrium (LD) on BTA14 were identified for fat percentage trait.
Five novel SNPs as candidate variants for fat percentage, protein percentage and cheese merit traits, also our haplotype-based analysis confirmed the region of potential association for these traits.
Genome-wide association study of milk components in Chinese Holstein cows using single nucleotide polymorphism Tianzhen Wang1,2, Jiao Li3, Xue Gao4, Wenqin Song2, Chengbin Chen2, Dawei Yao1, Jing Ma1,2, Lingyang Xu4,*, Yi Ma1,*
1. Tianjin Institute of Animal Science and Veterinary, Tianjin, 300381, China 2. College of Life Science, Nankai University, Tianjin, 300071, China 3. National Animal Husbandry Station, 100125, China 4 Innovation Team of Cattle Genetic Breeding, Institute of Animal Sciences, Chinese Academy of Agricultural Sciences, Beijing 100193, China
Email addresses: Tianzhen Wang: [email protected] Jiao Li: [email protected] Xue Gao: [email protected] Wenqin Song: [email protected] Chengbin Chen: [email protected] Dawei Yao: [email protected] Jing Ma: [email protected] Lingyang Xu: [email protected] Yi Ma: [email protected] *Corresponding Authors: Lingyang Xu: [email protected] & Yi Ma: [email protected] Running title: Genome-wide association for milk traits Address: 17 km, Jinjing Road, Xiqing District, Tianjin, 300381, P.R.China Tel: +86-022-83726967
Fax: +86-022-83726967
Abstract：Genome-wide association study (GWAS) has been widely utilized to identify candidate genes and causal loci for important traits in farm animals. In this study, we performed a GWAS for five milk- related traits using the GeneSeek Genomic Profiler Bovine 50K SNP chip in 300 Chinese Holstein cows. Using a mixed linear model, we detected 28 candidate SNPs surpassing the threshold level (P-value <1.2-6). A total of 18 SNPs located in Bos Taurus autosome (BTA)14 were significantly associated with at least one milk traits after Bonferroni correction. Interestingly, we observed a 2.5 Mb region (range from 1.5 to 4.0 Mb) with high linkage disequilibrium (LD) on BTA14 for fat percentage trait. Among them, four SNPs were identified for five traits including fat percentage, fat yield, protein percentage, milk yield and cheese merit, one SNP was associated with four traits fat percentage, fat yield, protein percentage and cheese merit, and nine SNPs were associated with fat percentage and fat yield traits. Moreover, our haplotype-based analysis confirmed the region of potential association for these traits. Notably, despite most identified SNPs have been reported as previous studies, we also found five novel SNPs as candidate variants for fat percentage, protein percentage and cheese merit traits.
Keyword: Milk traits, GWAS, haplotype-based analysis, mixed linear model, Chinese Holstein cows
1. Introduction Milk is an important source of nutrients as it plays a significant role in human nutrition all over the world. It is popular for its components that are also an important index for milk quality. Researchers have paid much attention to improve milk performance-related traits including fat yield (FY), fat percentage (FP), milk yield (MY), protein percentage (PP), and cheese merit (CM). CM is a composite index, which reflects the economic value of components like milk fat and protein. Genome-wide association study (GWAS) has been widely used to detect candidate genes for complex disease and quantitative traits in humans (Cordell, 2009) and farm animals (Goddard and Hayes, 2009). This approach has enabled us to understand genetic mechanisms underlying complex traits, and further promoted the rapid genetic improvement for these traits (Meuwissen et al., 2001; Zimin et al., 2009). With the development of high throughput technology, SNP array provides more opportunities for GWAS in farm animals. Several GWAS have detected loci and genes for milk
production traits in dairy cows. For instance, a previous study (Maxa et al., 2012) detected two SNPs for milk yield on BTA4, two SNPs for fat yield on BTA14 and BTA23, and one SNP was associated with fat percent on BTA1. Minozzi et al. (2013) identified 21 SNPs between 1.2 Mb and 2.8 Mb on BTA14 associated with fat yield and milk yield. Several studies identified a list of candidate SNPs or regions that were associated with multiple milk traits including milk production (Jiang et al., 2010), milk protein composition (Schopen et al., 2011), fat production, protein production and fat and protein deviation (Nayeri et al., 2016). Cole et al. (2011) identified a number of candidate genes, including GNAS for milk, fat and protein yields, DGAT-NIBP for fat percentage, FBP2 for protein yields and percentage, MGMT and PDGFRA for protein percentage. Moreover, many studies have suggested DGAT1 and TRPS1 located on BTA14 were associated with milk production, fat (Minozzi et al., 2013; Fang et al., 2017; Pausch et al., 2017; Pegolo et al., 2017; Ning et al., 2018) and lactation persistency (Do et al., 2017). For Chinese dairy cows, several studies in whole-genome level have investigated the genetic basis for economically important traits including milk fatty acid (Li et al., 2015), milk production (Ning et al., 2018), mastitis resistance (Cai et al., 2018), milk composition (Gao et al., 2017), pigmentation (Fan et al., 2014), female fertility (Liu et al., 2017) and body size (Zhang et al., 2017). Several candidate regions on BTA5, BTA6, BTA14 and BTA20 for milk production traits, fat percentage, milk yield, protein percentage have been identified (Ning et al., 2018). Previous findings also suggested several quantitative trait loci (QTLs) have pleiotropic effects for milk-related traits (Jiang et al., 2010; Pryce et al., 2010; Buitenhuis et al., 2014). However, the investigation of genetic mechanisms for new traits such as like cheese merit, milk composition and identification of the causal variants with pleiotropic effects have not been explored in Chinese dairy cows. The objectives of this study were to 1) detect the novel associated SNPs for five milk traits using the Mix Linear Model in Chinese Holstein cows, 2) evaluate the contribution of candidate variants with pleiotropic effects, and 3) identify candidate haplotypes for milk- related traits in Chinese Holstein cows.
2. Materials and Methods 2.1. Ethics statement All procedures strictly followed the guidelines developed by the China Council on Animal Care, and all protocols were approved by the Science Research Department of the Institute of Animal Science,
Tianjin Academy of Agricultural Science. We collected all the samples during regular quarantine inspections on the farms while all farm owners agreed to the use of these animals. 2.2. Sample genotyping and quality control Individuals (300 samples) from two farms in Tianjin were genotyped using the GeneSeek Genomic Profiler Bovine 50k SNP chip (Neogen Corporation, http://www.neogenchina.com.cn/). Genomic DNA was extracted from blood samples. The GeneSeek Genomic Profiler Bovine 50K SNP chip with 47,843 SNPs was utilized for individual genotyping. These SNPs were uniformly distributed on the genome with a mean inter-marker space of 59 kb. For quality control, individuals and their autosomal SNP were processed using PLINK v1.07 software (Purcell et al., 2007). First, all individuals with more than 10% missing SNP genotypes were excluded. Then, SNPs with either less than 95% call rate, more than 1% missing genotype, minor allele frequencies (MAF) less than 5%, and Hardy-Weinberg equilibrium (HWE) test P-value <10E-7 were also removed. 2.3. Predicted Transmitting Ability The association analysis was conducted by Predicted Transmitting Ability (PTA) which is the parent average and correction with the phenotype of offspring. It represents the expected performance of the daughters of a particular bull in relation to the average genetic herds (Cole et al., 2011; de Oliveira et al., 2014). PTA for five traits fat yield (FY), fat percentage (FP), protein percentage (PP), milk yield (MY) and cheese merit (CM) were calculated according to the standard method of CDCB (https://www.uscdcb.com/) from the U.S. Department of Agriculture (Beltsville, MD). 2.4. SNP-based GWAS using EMMAX Population stratification was assessed based on the inflation factor. Principal component analysis (PCA) was used to adjust the population stratification in GWAS for milk-related traits (Klei et al., 2008; Zhang et al., 2018). We performed GWAS for five traits (FP, FY, PP, MY and CM) using 40,501 high quality SNPs, the threshold for P-value was set to 0.05/36,532= 1.37E-06 (-log10(P) = 5.86) after Bonferroni correction. SNPs surpassing this threshold were considered as candidate variants. In this study, the mixed linear model implemented in EMMAX software was used to detect significant SNPs (Kang et al., 2008), and the equitation was as follow: y=Xβ+Zu+e
(1)
where, y is an n × 1 vector of PTA adjusted by non-genetic effects, and X is an n × q matrix of fixed
effect including SNPs and confounding variables. β is a q × 1 vector representing coefficients of the fixed effects. Z is an n × t incidence matrix relating phenotypes to the corresponding random polygenic effect. u is the random effect of the mixed model following a N (0, Iσ2g K), where K is the t × t kinship matrix. e is an n × n matrix of residual effect following a N (0, Iσ2g ). This approach can avoid repetitive variance component estimation procedure for association mapping using a mixed model (Kang et al., 2010). PCA was estimated within the population using PLINK (v1.9) (Chang et al., 2015). The statistical model for milk-related traits included the first two PCA components as covariates. Bonferroni correction for multiple tests was adopted for the value of the P-value at the genome significant of 0.05. Manhattan plot and quantile-quantile (Q-Q) plot were generated using qqman package version 0.1.4. All the statistical analyses were carried out based on R version 3.5.1 (The R Foundation, Vienna, Austria) 2.5. Candidate genes The SNPs coordinates were determined according to UMD3.1 (Zimin et al., 2009). We searched genes within 100 kb window size at upstream and downstream of associated SNPs (Table S1). In addition, we queried the significant SNPs for studied traits based on the QTLdb database (www.animalgenome.org/cgi-bin/QTLdb). Moreover, to evaluate the linkage disequilibrium (LD) pattern around the candidate regions, we estimated the parameter r2 using Haploview v4.2 software (Barrett et al., 2005). 2.6 Haplotype-based analysis The haplotype-based analysis was conducted for the target region in studied population. In this study, first, haplotype block was estimated using a standard expectation-maximum (EM) algorithm as described elsewhere (Barrett et al., 2005). Then, haplotype-based association analysis for five traits (FP, FY, PP, MY and CM) was performed using the command (--hap-assoc option) implemented in PLINK v1.07 software (Purcell et al., 2007). 3. Results 3.1. Quality control After quality control, one individual was removed due to more than 10% missing genotypes. A total of 4,523, 2,780 and 39 SNPs were deleted with missing genotyping rate > 0.01, MAF < 0.05 and HWE< 10E-7 respectively and 3969 SNPs in BTA0 and BTA30 were removed. In addition, 25 samples without PTA were removed due to incompletely pedigree information. Finally, 274 individuals with
36,532 high-quality SNPs were left for the subsequent GWAS analysis. 3.2. Population stratification assessment PCA based on the variance-standardized relationship matrix showed that the first two principal components (PC1 and PC2) contributed more than 93% of the variance. All individuals were clustered into three groups, as displayed in Fig S1, therefore, the population stratification based on PCA analysis result was considered and incorporated into the mixed linear model. 3.3. GWAS analysis In this study, 18 SNPs were significantly associated with five traits. The summary statistics of GWAS results are presented in Table 1 and the Manhattan plot for five traits is shown in Fig 1. We found the several most significant SNPs including Chr14_1757935 for FP (-log10(P) = 34.69), FY (-log10(P) = 17.57), PP (-log10(P) = 11.63), CM (-log10(P) = 7.84) and Chr14_1765835(-log10(P) = 8.54) for MY respectively. The region plot of P-values of 18 SNPs for each trait is shown in Fig 2a. For FP, we found 18 SNPs were significant associated with FP on BTA14.
Of those,
the
ARS-BFGL-NGS-94706, Chr14_1765835, BovineHD1400000246, Hapmap-24715-BTC-001-971, ARS-BFGL-NGS-26520, UA-IFASA-9288 located in the genes of VSP28, SLC52A2, MROH1, MAPK15, VAMP2, PTK2 (Table 1). For FY, nine SNPs were significantly associated with FY on BTA14.
These
SNPs
including
ARS-BFGL-NGS-
94706,
Chr14_1765835,
and
Hapmap24715-BTC-001973 were located in three genes (VPS28, SLC52A2, and MAPK15). For PP, MY and CM, we found four significant SNPs associations on BTA14. Of those, five SNPs including ARS-BFGL-NGS-57820, Chr14_1653693, Chr14_1757935, Chr14_1765835 and Chr14_2022745 were embedded or neighbored with FOXH1, ADCK5, SLC52A2, GRINA genes respectively. All the 18 significant SNPs were found embedded with discovered QTLs as reported in previous studies (Table S2). (Jiang et al., 2010; Pryce et al., 2010; Cole et al., 2011; Marques et al., 2011; Minozzi et al., 2013; Wang et al., 2013; Buitenhuis et al., 2014; Jiang et al., 2014; Capomaccio et al., 2015; Lehnert et al., 2015; Buitenhuis et al., 2016; Iso-Touru et al., 2016). These identified SNPs associated with multiple milk traits (FP, FY, MY, and PP), which may have pleiotropic effects related to milk traits in dairy cows. We found several novels QTLs for milk traits, for instance, the candidate region with 0.37 Mb (range 1.65 to 2.02 Mb) on BTA14 associated with three traits including FP, PP, CM, and ARS-BFGL-NGS-57820 was detected as a novel candidate for CM in the current study. 3.4. LD analysis of candidate region
To further determine correction between significant SNPs, we performed LD analysis in 2.5 Mb region (contains 31 SNPs) by computing pairwise r2 from Hapmap30383-BTC-005848 to UA-IFASA-9288 (1.5~4.0 Mb) on BTA14. As a result, we found two haploblocks (Fig 2b) with ~198 kb size
in
this region.
ARS-BFGL-NGS-57820,
Among them,
Chr14_1653693,
one
block had four
SNPs
BTA-34956-no-rs,
ARS-BFGL-NGS-94706,
while
the
five
SNPs
(Chr14_1757-935, Chr14_1765835, BovineHD1400000246, ARS-BFGL- NGS-71749, Chr14_202 2745) were found in other block. Moreover, we observed nine significant SNPs in candidate region with 35 candidate genes C14H8orf33, ZNF34, RPL8, COMMD5, ARHGAP39, MIR2306, C14H8orf82, LRRC24, LRRC14, RECQL4, MFSD3, GPT, PPIR16A, FOXH1, KIFC2, CYHR1, TONSL, VPS28, SLC39A4, CPSF1, ADCK5, SLC52A2, SCRT1, DGAT1, HSF1, MROH1, MIR1839, MAF1, SHARPIN, CYC1, GPAA1, EXOSC4, OPLAH, GRINA and MIR2309 respectively (Table S1). For CM and MY, regional plots of candidate region (1.33 to 2.60 Mb) on BTA14 were presented in Fig 3. The most two significant SNPs are Chr14_1757935 and Chr14_1765835 (-log10(P) = 7.84, 8.54), we also found nine SNPs with strong LD around them (CM: ARS-BFGL-NGS-57820, Chr14_1653693, Chr14_1765835, Chr14_2022745. MY: ARS-BFGL-NGS-57820, Chr14_1653693, Chr14_1757935). 3.5 Haplotype-based analysis The target region located at BTA14 (range from 1.5 to 4.0 Mb) were selected for haplotype –based analysis. Two blocks were identified with eight haplotypes across the studied population (Table S3), and these blocks were comprised of seven SNPs. Of those, the length of the large block was 190 kb, while 7.9 kb for the small block. In this study, we found three significant haplotypes (shown in Table 2) for each trait at the threshold level (P-value < 10-5). Among them, significant haplotypes H1-1 (AAGCT), H2-1 (TG), H2-3 (CC) were identified for FP, FY, MY, PP, while the most significant haplotype was found for FP trait (-log10(P) =21.88, 26.09,28.19). Notably, three haplotypes including H2-1 (-log10(P) = 4.29) and H2-3 (-log10(P) = 4.25), and H1-5 (AGGTC, -log10(P) = 4.25) showed significant association for CM. 4. Discussion In this study, we performed a GWAS for five important milk-related traits (fat percentage, fat yield, protein percentage, milk yield and cheese merit) using GeneSeek Genomic Profiler Bovine 50K SNP chip in Chinese Holstein cows. Some of these traits have been explored in previous studies; however, elucidation of the molecular mechanism underpinning milk traits in other populations can offer
additional insights for understanding the genetic basis of these traits in dairy cows. In the present study, we first estimated the population stratification in our population. Using PCA analysis implemented in PLINK v1.9 (Chang et al., 2015), we observed our data can be clustered into three groups, and the inflation factor (λ) suggested obvious population stratification for the studied population. Therefore, population stratification and SNP effects were considered as fixed effects, and the polygenic effects were considered as random factors in mix linear model (Banerjee et al., 2008; Ferreira and Purcell, 2009; Kim and Xing, 2009; O'Reilly et al., 2012). These methods not only can increase power for detecting pleiotropic genetic variants but also avoid false-positive result in association tests (Stephens, 2013). Previous studies found several regions and genes on BTA14 were widely associated with many important traits based on multiple statistical association tests. For instance, DGAT1 located at ~1.8Mb on BTA14 is generally accepted as a major gene affecting milk production traits (milk yield, protein content, fat content, and fatty acid composition) (Grisart et al., 2002; Schennink et al., 2007; Schennink et al., 2008; Conte et al., 2010). CYP11B1 located at ~1.33 Mb has significant effects on milk yield, protein yield, fat percentage, protein percentage (Kaupe et al., 2007). In this study, we detected 18 significant SNPs for milk-related traits on BTA14. Based these SNPs, we obtained 29 candidate genes including C14H8orf33, FOXH1, CYHR1, VPS28, SCRT1, DGAT1, OPLAH, MROH1, MAF1, GRINA, LY6D, LY6H, ZNF34, RPL8, GPT, PPP1R16A, CPSF1, HSF1, GPAA1, PUF60, MAPK15, EEF1D, LYNX1, PTK2, ARHGAP39, TONSL, LRRC14, ADCK5, which have been reported in previous studies. Among them, 22 genes (C14H8orf33, FOXH1, CYHR1, VPS28, DGAT1, OPLAH, MAF1, GRINA, LY6D, LY6H, ZNF34, RPL8, GPT, PPP1R16A, CPSF1, HSF1, GPAA1, PUF60, MAPK15, EEF1D, LYNX1, PTK2) have been detected for multiple milk traits in Chinese dairy cows (Jiang et al., 2010; Jiang et al., 2014). ARHGAP39 was identified as a candidate gene for mastitis susceptibility in Holsteins (Wang et al., 2015). TONSL, LRRC14, ADCK5 were detected as potential candidate genes for milk traits or mammary gland functions (Ibeagha-Awemu et al., 2016). MROH1 and SCRT1 were affected by FP (Capomaccio et al., 2015). In our study, we found several new candidate genes for milk-related traits. For instance, CPSF1, DGAT1, GPT, PPP1R16A, and CYHR1 were associated with CM. HSF1 was associated with FY, PP, MY, CM. ZNF34 was associated with FY, ARHGAP39, TSTA3 for FP and GRINA for FY and CM (Table S1). We also found several new candidate genes including SLC52A2, VAMP2 embedded with
three SNPs (Chr14_1765835 and ARS-BFGL-NGS-26520). Of these, SLC52A2 was associated with FP, FY, PP, MY, CM and VAMP2 was associated with FP (Table 1). In addition, we observed many significant SNPs (BTA-34956-no-rs, ARS-BFGL-NGS-57820, Chr14_1653693, Chr14_1765835, ARS-BFGL-NGS-71749,
Hapmap24715-BTC-001973,
ARS-BFGL-NGS-26520,
ARS-BFGL-NGS-22866) associated with milk traits nearside several candidate genes (Table S1). Our study identified 18 significant SNPs region located at ~2.5 Mb (from 1.65 to 2.02 Mb) on BTA14 associated with multi-traits. Of these, we found several new variants (ARS-BFGL-NGS-57820, Chr14_1653693, Chr14_1757935, Chr14_1765835, Chr14_2022745) for FP, PP and CM located within a 0.37Mb (between 1.65 to 2.02 Mb) on BTA14. Previous studies suggested many QTLs with pleiotropic effects related to milk-related traits in dairy cows (Olsen et al., 2011; Bolormaa et al., 2014; Saatchi et al., 2014). Our study identified 18 significant SNPs region located at ~2.5 Mb (from 1.65 to 2.02 Mb) on BTA14 associated with multi-traits. Of these, we found several new variants (ARS-BFGL-NGS-57820, Chr14_1653693, Chr14_1757935, Chr14_1765835, Chr14_2022745) for FP, PP and CM located within a 0.37Mb (between 1.65 to 2.02 Mb) on BTA14. Moreover, Haplotype-base analysis further confirmed the results of our SNP-based GWAS, we found three significant haplotypes for multi traits, and these haplotypes contain the identified SNPs, which have been detected by SNP-based GWAS. Haplotype H2-1, H2-3 were associated with the five traits (FP, FY, PP, MY, CM) and H1-1 was associated with four traits (FP, FY, PP, MY). The identification of candidate SNPs and genes can help improve our knowledge of milk traits and provide useful information for guarding the design of breeding programs. However, more precise genomic analysis and validation of functional variations with larger sample size are required in further study. 5. Conclusions Our study identified eighteen candidate SNPs for five traits in Chinese Holstein cows, and five of them were novel discovered in the current study. These findings suggested that nine candidate SNPs may have pleiotropic effects for milk-related traits. Our results may further help to promote the design of marker-assisted selection for a breeding program. But this candidate genes must be a further
investigated for the improvement of milk traits.
Author Contributions: Conceptualization, Wenqin Song, Lingyang Xu and Yi Ma; Data curation, Tianzhen Wang,
Chengbin Chen, Dawei Yao and Jing Ma; Investigation, Tianzhen Wang and Jiao Li; Methodology, Xue Gao and Lingyang Xu; Project administration, Yi Ma; Resources, Yi Ma; Writing – original draft, Tianzhen Wang; Writing – review & editing, Tianzhen Wang and Lingyang Xu.
Funding: This research was funded by the Innovation Team of Tianjin Cattle Research System (ITTCRS2017005) and Tianjin Seed Industry Special Project (16ZXZYNC00220).
Acknowledgments: The authors would like to thank the breeders and breeding associations who provide us samples. The SNPs data and PTA data were provided by Neogen Corporation ([email protected])
Declarations Ethics and consent to participate All animals used in the study were treated following the guidelines for the experimental animals established by the Council of China Animal Welfare. Animal experiments were approved by the Science Research Department of the Institute of Animal Science, Tianjin Academy of Agricultural Science. All farm owners agreed on the use of the animals.
Consent to publish Not applicable
Competing interests The authors declare that they have no competing interests. Availability of data The data and compute programs used in this manuscript are available from the corresponding author on request.
Figure legends
Figure 1. a) Manhattan plot of -log10(P) for FY in Chinese Holsteins dairy cows. Chromosomes are color-coded separately. The horizontal reference line indicated the genome-wide significance
levels (-log10(P) = 5.98). b) Q-Q plots for FY. Deviations from the slope lines correspond to loci that deviate from the null hypothesis. c) Manhattan plot of -log10(P) for FP. d) Q-Q plots for FP. e) Manhattan plot of -log10(P) for MY. f) Q-Q plots for five MY. g) Manhattan plot of -log10(P) for PP. h) Q-Q plots for PP. i) Manhattan plot of -log10(P) for CM. j) Q-Q plots for CM.
Figure 2. a) Association test with -log10(P) for 31 SNPs for five traits including FP, FY, PP, MY and CM. The horizontal reference line indicated the genome-wide significance levels (-log10(P) = 5.98), 18 of 31 were significant. b) Haplotype analysis for the region around BTA14: 1489496– 3956956 bp. Two haplotypes with 182 and 264 kb were estimated using Haploview software.
Figure 3. a) Regional plots of candidate region at 1.33 to 2.60 Mb on BTA14 for CM. b) Regional plots of candidate region at 1.33 to 2.60 Mb on BTA14 for MY.
All authors agree to publish.
Conflict of interest The authors declare that they have no competing interests.
References Archer-Nicholls, S., Lowe, D., Utembe, S., Allan, J., Zaveri, R.A., Fast, J.D., Hodnebrog, O., van der Gon, H.D., McFiggans, G., 2014. Gaseous chemistry and aerosol mechanism developments for version 3.5.1 of the online regional model, WRF-Chem. Geoscientific Model Development 7, 2557-2579. Banerjee, S., Yandell, B.S., Yi, N., 2008. Bayesian quantitative trait loci mapping for multiple traits. Genetics 179, 2275-2289. Barrett, J.C., Fry, B., Maller, J., Daly, M.J., 2005. Haploview: analysis and visualization of LD and haplotype maps. Bioinformatics 21, 263-265. Bolormaa, S., Pryce, J.E., Reverter, A., Zhang, Y., Barendse, W., Kemper, K., Tier, B., Savin, K., Hayes, B.J., Goddard, M.E., 2014. A multi-trait, meta-analysis for detecting pleiotropic polymorphisms for stature, fatness and reproduction in beef cattle. PLoS Genetics 10, 3. Buitenhuis, B., Janss, L.L., Poulsen, N.A., Larsen, L.B., Larsen, M.K., Sorensen, P., 2014. Genome-wide association and biological pathway analysis for milk-fat composition in Danish Holstein and Danish Jersey cattle. BMC Genomics 15, 1112. Buitenhuis, B., Poulsen, N.A., Gebreyesus, G., Larsen, L.B., 2016. Estimation of genetic parameters and detection of chromosomal regions affecting the major milk proteins and their post translational modifications in Danish Holstein and Danish Jersey cattle. BMC Genetics 17, 114. Cai, Z., Guldbrandtsen, B., Lund, M.S., Sahana, G., 2018. Prioritizing candidate genes post-GWAS using multiple sources of data for mastitis resistance in dairy cattle. BMC Genomics 19, 656. Capomaccio, S., Milanesi, M., Bomba, L., Cappelli, K., Nicolazzi, E.L., Williams, J.L., Ajmone-Marsan, P., Stefanon, B., 2015. Searching new signals for production traits through gene-based association analysis in three Italian cattle breeds. Animal Genetics 46, 361-370. Chang, C.C., Chow, C.C., Tellier, L.C., Vattikuti, S., Purcell, S.M., Lee, J.J., 2015. Second-generation PLINK: rising to the challenge of larger and richer datasets. Gigascience 4, 7. Cole, J.B., Wiggans, G.R., Ma, L., Sonstegard, T.S., Lawlor, T.J., Jr., Crooker, B.A., Van Tassell, C.P., Yang, J., Wang, S., Matukumalli, L.K., Da, Y., 2011. Genome-wide association analysis of thirty one production, health, reproduction and body conformation traits in contemporary U.S. Holstein cows. BMC Genomics 12, 408. Conte, G., Mele, M., Chessa, S., Castiglioni, B., Serra, A., Pagnacco, G., Secchiari, P., 2010. Diacylglycerol acyltransferase 1, stearoyl-CoA desaturase 1, and sterol regulatory element binding protein 1 gene polymorphisms and milk fatty acid composition in Italian Brown cattle. Journal of Dairy
Science 93, 753-763. Cordell, H.J., 2009. Detecting gene-gene interactions that underlie human diseases. Nature Reviews. Genetics 10, 392-404. de Oliveira, F.C., Borges, C.C., Almeida, F.N., e Silva, F.F., da Silva Verneque, R., da Silva, M.V., Arbex, W., 2014. SNPs selection using support vector regression and genetic algorithms in GWAS. BMC Genomics 15(7), 1-15. Do, D.N., Bissonnette, N., Lacasse, P., Miglior, F., Sargolzaei, M., Zhao, X., Ibeagha-Awemu, E.M., 2017. Genome-wide association analysis and pathways enrichment for lactation persistency in Canadian Holstein cattle. Journal of Dairy Science 100, 1955-1970. Fan, Y., Wang, P., Fu, W., Dong, T., Qi, C., Liu, L., Guo, G., Li, C., Cui, X., Zhang, S., Zhang, Q., Zhang, Y., Sun, D., 2014. Genome-wide association study for pigmentation traits in Chinese Holstein population. Animal Genetics 45, 740-744. Fang, L., Sahana, G., Su, G., Yu, Y., Zhang, S., Lund, M.S., Sorensen, P., 2017. Integrating Sequence-based GWAS and RNA-Seq Provides Novel Insights into the Genetic Basis of Mastitis and Milk Production in Dairy Cattle. Scientific Reports 7, 45560. Ferreira, M.A., Purcell, S.M., 2009. A multivariate test of association. Bioinformatics 25, 132-133. Gao, Y., Jiang, J., Yang, S., Hou, Y., Liu, G.E., Zhang, S., Zhang, Q., Sun, D., 2017. CNV discovery for milk composition traits in dairy cattle using whole genome resequencing. BMC Genomics 18, 265. Goddard, M.E., Hayes, B.J., 2009. Mapping genes for complex traits in domestic animals and their use in breeding programmes. Nature reviews. Genetics 10, 381-391. Grisart, B., Coppieters, W., Farnir, F., Karim, L., Ford, C., Berzi, P., Cambisano, N., Mni, M., Reid, S., Simon, P., Spelman, R., Georges, M., Snell, R., 2002. Positional candidate cloning of a QTL in dairy cattle: identification of a missense mutation in the bovine DGAT1 gene with major effect on milk yield and composition. Genome Research 12, 222-231. Ibeagha-Awemu, E.M., Peters, S.O., Akwanji, K.A., Imumorin, I.G., Zhao, X., 2016. High density genome- wide genotyping-by-sequencing and association identifies common and low-frequency SNPs, and novel candidate genes influencing cow milk traits. Scientific Reports 6. Iso-Touru, T., Sahana, G., Guldbrandtsen, B., Lund, M.S., Vilkki, J., 2016. Genome-wide association analysis of milk yield traits in Nordic Red Cattle using imputed whole genome sequence variants. BMC Genetics 17, 55. Jiang, L., Liu, J., Sun, D., Ma, P., Ding, X., Yu, Y., Zhang, Q., 2010. Genome wide association studies for milk production traits in Chinese Holstein population. PloS One 5, 10. Jiang, L., Liu, X., Yang, J., Wang, H., Jiang, J., Liu, L., He, S., Ding, X., Liu, J., Zhang, Q., 2014. Targeted resequencing of GWAS loci reveals novel genetic variants for milk production traits. BMC Genomics 15, 1105. Kang, H.M., Sul, J.H., Service, S.K., Zaitlen, N.A., Kong, S.Y., Freimer, N.B., Sabatti, C., Eskin, E., 2010. Variance component model to account for sample structure in genome-wide association studies. Nature Genetics 42, 348-354. Kang, H.M., Zaitlen, N.A., Wade, C.M., Kirby, A., Heckerman, D., Daly, M.J., Eskin, E., 2008. Efficient control of population structure in model organism association mapping. Genetics 178, 1709-1723. Kaupe, B., Brandt, H., Prinzenberg, E.M., Erhardt, G., 2007. Joint analysis of the influence of CYP11B1 and DGAT1 genetic variation on milk production, somatic cell score, conformation, reproduction, and productive lifespan in German Holstein cattle. Journal of Animal Science 85, 11-21. Kim, S., Xing, E.P., 2009. Statistical estimation of correlated genome associations to a quantitative trait
network. PLoS Genetics 5, 8. Klei, L., Luca, D., Devlin, B., Roeder, K., 2008. Pleiotropy and principal components of heritability combine to increase power for association analysis. Genetic Epidemiology 32, 9-19. Lehnert, K., Ward, H., Berry, S.D., Ankersmit-Udy, A., Burrett, A., Beattie, E.M., Thomas, N.L., Harris, B., Ford, C.A., Browning, S.R., Rawson, P., Verkerk, G.A., van der Does, Y., Adams, L.F., Davis, S.R., Jordan, T.W., MacGibbon, A.K., Spelman, R.J., Snell, R.G., 2015. Phenotypic population screen identifies a new mutation in bovine DGAT1 responsible for unsaturated milk fat. Scientific Reports 5, 8484. Li, X., Buitenhuis, A.J., Lund, M.S., Li, C., Sun, D., Zhang, Q., Poulsen, N.A., Su, G., 2015. Joint genome-wide association study for milk fatty acid traits in Chinese and Danish Holstein populations. Journal of Dairy Science 98, 8152-8163. Liu, A., Wang, Y., Sahana, G., Zhang, Q., Liu, L., Lund, M.S., Su, G., 2017. Genome-wide Association Studies for Female Fertility Traits in Chinese and Nordic Holsteins. Scientific Reports 7, 8487. Marques, E., Grant, J.R., Wang, Z., Kolbehdari, D., Stothard, P., Plastow, G., Moore, S.S., 2011. Identification of candidate markers on bovine chromosome 14 (BTA14) under milk production trait quantitative trait loci in Holstein. Journal of animal breeding and genetics = Zeitschrift fur Tierzuchtung und Zuchtungsbiologie 128, 305-313. Maxa, J., Neuditschko, M., Russ, I., Forster, M., Medugorac, I., 2012. Genome-wide association mapping of milk production traits in Braunvieh cattle. Journal of Dairy Science 95, 5357-5364. Meuwissen, T.H., Hayes, B.J., Goddard, M.E., 2001. Prediction of total genetic value using genome-wide dense marker maps. Genetics 157, 1819-1829. Minozzi, G., Nicolazzi, E.L., Stella, A., Biffani, S., Negrini, R., Lazzari, B., Ajmone-Marsan, P., Williams, J.L., 2013. Genome wide analysis of fertility and production traits in Italian Holstein cattle. PloS One 8, 1-10. Nayeri, S., Sargolzaei, M., Abo-Ismail, M.K., May, N., Miller, S.P., Schenkel, F., Moore, S.S., Stothard, P., 2016. Genome-wide association for milk production and female fertility traits in Canadian dairy Holstein cattle. BMC Genetics 17, 75. Ning, C., Wang, D., Zheng, X., Zhang, Q., Zhang, S., Mrode, R., Liu, J.F., 2018. Eigen decomposition expedites longitudinal genome-wide association studies for milk production traits in Chinese Holstein. Genetics, Selection, Evolution : GSE 50, 12. O'Reilly, P.F., Hoggart, C.J., Pomyen, Y., Calboli, F.C., Elliott, P., Jarvelin, M.R., Coin, L.J., 2012. MultiPhen: joint model of multiple phenotypes can increase discovery in GWAS. PloS One 7, e34861. Olsen, H.G., Hayes, B.J., Kent, M.P., Nome, T., Svendsen, M., Larsgard, A.G., Lien, S., 2011. Genome-wide association mapping in Norwegian Red cattle identifies quantitative trait loci for fertility and milk production on BTA12. Animal Genetics 42, 466-474. Pausch, H., Emmerling, R., Gredler-Grandl, B., Fries, R., Daetwyler, H.D., Goddard, M.E., 2017. Meta-analysis of sequence-based association studies across three cattle breeds reveals 25 QTL for fat and protein percentages in milk at nucleotide resolution. BMC Genomics 18, 853. Pegolo, S., Dadousis, C., Mach, N., Ramayo-Caldas, Y., Mele, M., Conte, G., Schiavon, S., Bittante, G., Cecchinato, A., 2017. SNP co-association and network analyses identify E2F3, KDM5A and BACH2 as key regulators of the bovine milk fatty acid profile. Scientific Reports 7, 17317. Pryce, J.E., Bolormaa, S., Chamberlain, A.J., Bowman, P.J., Savin, K., Goddard, M.E., Hayes, B.J., 2010. A validated genome-wide association study in 2 dairy cattle breeds for milk production and fertility traits using variable length haplotypes. Journal of Dairy Science 93, 3331-3345. Purcell, S., Neale, B., Todd-Brown, K., Thomas, L., Ferreira, M.A., Bender, D., Maller, J., Sklar, P., de
Bakker, P.I., Daly, M.J., Sham, P.C., 2007. PLINK: a tool set for whole-genome association and population-based linkage analyses. American Journal of Human Genetics 81, 559-575. Saatchi, M., Schnabel, R.D., Taylor, J.F., Garrick, D.J., 2014. Large-effect pleiotropic or closely linked QTL segregate within and across ten US cattle breeds. BMC Genomics 15, 442. Schennink, A., Heck, J.M., Bovenhuis, H., Visker, M.H., van Valenberg, H.J., van Arendonk, J.A., 2008. Milk fatty acid unsaturation: genetic parameters and effects of stearoyl-CoA desaturase (SCD1) and acyl CoA: diacylglycerol acyltransferase 1 (DGAT1). Journal of Dairy Science 91, 2135-2143. Schennink, A., Stoop, W.M., Visker, M.H., Heck, J.M., Bovenhuis, H., van der Poel, J.J., van Valenberg, H.J., van Arendonk, J.A., 2007. DGAT1 underlies large genetic variation in milk-fat composition of dairy cows. Animal Genetics 38, 467-473. Schopen, G.C., Visker, M.H., Koks, P.D., Mullaart, E., van Arendonk, J.A., Bovenhuis, H., 2011. Whole-genome association study for milk protein composition in dairy cattle. Journal of Dairy Science 94, 3148-3158. Stephens, M., 2013. A unified framework for association analysis with multiple related phenotypes. PloS One 8, e65245. Wang, H., Jiang, L., Liu, X., Yang, J., Wei, J., Xu, J., Zhang, Q., Liu, J.F., 2013. A post-GWAS replication study confirming the PTK2 gene associated with milk production traits in Chinese Holstein. PloS One 8, e83625. Wang, X., Ma, P.P., Liu, J.F., Zhang, Q., Zhang, Y., Ding, X.D., Jiang, L., Wang, Y.C., Zhang, Y., Sun, D.X., Zhang, S.L., Su, G.S., Yu, Y., 2015. Genome-wide association study in Chinese Holstein cows reveal two candidate genes for somatic cell score as an indicator for mastitis susceptibility. BMC Genetics 16. Zhang, W., Gao, X., Shi, X., Zhu, B., Wang, Z., Gao, H., Xu, L., Zhang, L., Li, J., Chen, Y., 2018. PCA-based multiple-trait GWAS analysis: A powerful model for exploring pleiotropy. Animal 8, 12. Zhang, X., Chu, Q., Guo, G., Dong, G., Li, X., Zhang, Q., Zhang, S., Zhang, Z., Wang, Y., 2017. Genome-wide association studies identified multiple genetic loci for body size at four growth stages in Chinese Holstein cattle. PloS One 12, e0175971. Zimin, A.V., Delcher, A.L., Florea, L., Kelley, D.R., Schatz, M.C., Puiu, D., Hanrahan, F., Pertea, G., Van Tassell, C.P., Sonstegard, T.S., Marcais, G., Roberts, M., Subramanian, P., Yorke, J.A., Salzberg, S.L., 2009. A whole-genome assembly of the domestic cow, Bos taurus. Genome Biology 10(4), 1-10.
Table 1. Candidate SNPs list associated with milk traits fat percentage (FP), fat yield (FY), protein percentage (PP), milk yield (MY), cheese merit (CM)) on BTA14. SNP
Position
Gene
Distance
MAF
-log10(P)
Traits
Hapmap30383-B
1489496
C14H8orf33
78
0.3571
14.55
FP
BTA-34956-no-rs
1514056
RPL8
-6430
0.4891
6.75
FP
ARS-BFGL-NGS
1651311
FOXH1
3390
0.2719
27.13
FP
FY
PP
MY
CM
Chr14_1653693
1653693
FOXH1
1008
0.2646
27.76
FP
FY
PP
MY
CM
ARS-BFGL-NGS
1696470
VPS28
in
0.4088
13.77
FP
FY
Chr14_1757935
1757935
ADCK5
-1634
0.2527
33.49
FP
FY
PP
MY
CM
Chr14_1765835
1765835
SLC52A2
in
0.2792
32.55
FP
FY
PP
MY
CM
BovineHD140000
1880378
MROH1
in
0.4487
8.61
FP
1954317
OPLAH
3427
0.1752
6.94
FP
FY
PP
FY
TC-005848
-57820
-94706
0246 ARS-BFGL-NGS -71749 Chr14_2022745
2022745
GRINA
-1036
0.2719
28.36
FP
Hapmap25384-B
2217163
MAPK15
17870
0.4744
9.02
FP
2239085
MAPK15
in
0.4707
9.52
FP
BTA-35941-no-rs
2276443
MIR193A-2
4374
0.4635
10.23
FP
ARS-BFGL-NGS
2386688
VAMP2
in
0.2482
8.92
FP
2524432
LY6H
38728
0.4234
11.16
FP
2553525
LY6H
9635
0.3704
9.85
FP
2826632
LYNX1
-9946
0.3485
6.18
FP
3956956
PTK2
in
0.293
7.68
FP
TC-001997 Hapmap24715-B
FY
TC-001973
-26520 Hapmap30086-B TC-002066 Hapmap30646-B TC-002054 ARS-BFGL-NGS -22866 UA-IFASA-9288
MAF: minor allele frequencies
FY
CM
Table 2. Significant haplotypes from the haplotype-based analysis. H1-1, H1-5 is the first block with Hapmap30381-BTC-005750,
Hapmap30383-BTC-005848,
BTA-34956-no-rs,
ARS-BFGL-NGS-57820,
Chr14_1653693, and H2-1, H2-3 is the second block with Chr14_1757935, Chr14_1765835 Trait
PP
MY
FY
FP
CM
Hapname
Haplotype
Start
End
Length/bp.
P-value
H1-1
AAGCT
1463676
1653693
190017
1.18E-08
H2-1
TG
1757935
1765835
7900
1.41E-09
H2-3
CC
1757935
1765835
7900
2.56E-10
H1-1
AAGCT
1463676
1653693
190017
4.34E-06
H2-1
TG
1757935
1765835
7900
2.25E-06
H2-3
CC
1757935
1765835
7900
1.11E-07
H1-1
AAGCT
1463676
1653693
190017
6.73E-08
H2-1
TG
1757935
1765835
7900
8.13E-10
H2-3
CC
1757935
1765835
7900
3.21E-09
H1-1
AAGCT
1463676
1653693
190017
1.31E-22
H2-1
TG
1757935
1765835
7900
8.06E-27
H2-3
CC
1757935
1765835
7900
6.47E-29
H1-5
AGGTC
1463676
1653693
190017
5.66E-05
H2-1
TG
1757935
1765835
7900
5.16E-05
H2-3
CC
1757935
1765835
7900
5.64E-05
Genome-wide association study of milk components in Chinese Holstein cows using single nucleotide polymorphism Tianzhen Wang , Jiao Li , Xue Gao , Wenqin Song , Chengbin Chen , Dawei Yao , Jing Ma , Lingyang Xu , Yi Ma PII: DOI: Reference:
S1871-1413(19)31077-7 https://doi.org/10.1016/j.livsci.2020.103951 LIVSCI 103951
To appear in:
Livestock Science
Received date: Revised date: Accepted date:
2 August 2019 19 January 2020 28 January 2020
Please cite this article as: Tianzhen Wang , Jiao Li , Xue Gao , Wenqin Song , Chengbin Chen , Dawei Yao , Jing Ma , Lingyang Xu , Yi Ma , Genome-wide association study of milk components in Chinese Holstein cows using single nucleotide polymorphism, Livestock Science (2020), doi: https://doi.org/10.1016/j.livsci.2020.103951
This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2020 Published by Elsevier B.V.
Highlights
A total of 18 SNPs located in BTA14 were detected which were associated with at least one milk traits using mixed linear model after Bonferroni correction.
One 2.5 Mb region (range from 1.5 to 4.0 Mb) with high linkage disequilibrium (LD) on BTA14 were identified for fat percentage trait.
Five novel SNPs as candidate variants for fat percentage, protein percentage and cheese merit traits, also our haplotype-based analysis confirmed the region of potential association for these traits.
Genome-wide association study of milk components in Chinese Holstein cows using single nucleotide polymorphism Tianzhen Wang1,2, Jiao Li3, Xue Gao4, Wenqin Song2, Chengbin Chen2, Dawei Yao1, Jing Ma1,2, Lingyang Xu4,*, Yi Ma1,*
1. Tianjin Institute of Animal Science and Veterinary, Tianjin, 300381, China 2. College of Life Science, Nankai University, Tianjin, 300071, China 3. National Animal Husbandry Station, 100125, China 4 Innovation Team of Cattle Genetic Breeding, Institute of Animal Sciences, Chinese Academy of Agricultural Sciences, Beijing 100193, China
Email addresses: Tianzhen Wang: [email protected] Jiao Li: [email protected] Xue Gao: [email protected] Wenqin Song: [email protected] Chengbin Chen: [email protected] Dawei Yao: [email protected] Jing Ma: [email protected] Lingyang Xu: [email protected] Yi Ma: [email protected] *Corresponding Authors: Lingyang Xu: [email protected] & Yi Ma: [email protected] Running title: Genome-wide association for milk traits Address: 17 km, Jinjing Road, Xiqing District, Tianjin, 300381, P.R.China Tel: +86-022-83726967
Fax: +86-022-83726967
Abstract：Genome-wide association study (GWAS) has been widely utilized to identify candidate genes and causal loci for important traits in farm animals. In this study, we performed a GWAS for five milk- related traits using the GeneSeek Genomic Profiler Bovine 50K SNP chip in 300 Chinese Holstein cows. Using a mixed linear model, we detected 28 candidate SNPs surpassing the threshold level (P-value <1.2-6). A total of 18 SNPs located in Bos Taurus autosome (BTA)14 were significantly associated with at least one milk traits after Bonferroni correction. Interestingly, we observed a 2.5 Mb region (range from 1.5 to 4.0 Mb) with high linkage disequilibrium (LD) on BTA14 for fat percentage trait. Among them, four SNPs were identified for five traits including fat percentage, fat yield, protein percentage, milk yield and cheese merit, one SNP was associated with four traits fat percentage, fat yield, protein percentage and cheese merit, and nine SNPs were associated with fat percentage and fat yield traits. Moreover, our haplotype-based analysis confirmed the region of potential association for these traits. Notably, despite most identified SNPs have been reported as previous studies, we also found five novel SNPs as candidate variants for fat percentage, protein percentage and cheese merit traits.
Keyword: Milk traits, GWAS, haplotype-based analysis, mixed linear model, Chinese Holstein cows
1. Introduction Milk is an important source of nutrients as it plays a significant role in human nutrition all over the world. It is popular for its components that are also an important index for milk quality. Researchers have paid much attention to improve milk performance-related traits including fat yield (FY), fat percentage (FP), milk yield (MY), protein percentage (PP), and cheese merit (CM). CM is a composite index, which reflects the economic value of components like milk fat and protein. Genome-wide association study (GWAS) has been widely used to detect candidate genes for complex disease and quantitative traits in humans (Cordell, 2009) and farm animals (Goddard and Hayes, 2009). This approach has enabled us to understand genetic mechanisms underlying complex traits, and further promoted the rapid genetic improvement for these traits (Meuwissen et al., 2001; Zimin et al., 2009). With the development of high throughput technology, SNP array provides more opportunities for GWAS in farm animals. Several GWAS have detected loci and genes for milk
production traits in dairy cows. For instance, a previous study (Maxa et al., 2012) detected two SNPs for milk yield on BTA4, two SNPs for fat yield on BTA14 and BTA23, and one SNP was associated with fat percent on BTA1. Minozzi et al. (2013) identified 21 SNPs between 1.2 Mb and 2.8 Mb on BTA14 associated with fat yield and milk yield. Several studies identified a list of candidate SNPs or regions that were associated with multiple milk traits including milk production (Jiang et al., 2010), milk protein composition (Schopen et al., 2011), fat production, protein production and fat and protein deviation (Nayeri et al., 2016). Cole et al. (2011) identified a number of candidate genes, including GNAS for milk, fat and protein yields, DGAT-NIBP for fat percentage, FBP2 for protein yields and percentage, MGMT and PDGFRA for protein percentage. Moreover, many studies have suggested DGAT1 and TRPS1 located on BTA14 were associated with milk production, fat (Minozzi et al., 2013; Fang et al., 2017; Pausch et al., 2017; Pegolo et al., 2017; Ning et al., 2018) and lactation persistency (Do et al., 2017). For Chinese dairy cows, several studies in whole-genome level have investigated the genetic basis for economically important traits including milk fatty acid (Li et al., 2015), milk production (Ning et al., 2018), mastitis resistance (Cai et al., 2018), milk composition (Gao et al., 2017), pigmentation (Fan et al., 2014), female fertility (Liu et al., 2017) and body size (Zhang et al., 2017). Several candidate regions on BTA5, BTA6, BTA14 and BTA20 for milk production traits, fat percentage, milk yield, protein percentage have been identified (Ning et al., 2018). Previous findings also suggested several quantitative trait loci (QTLs) have pleiotropic effects for milk-related traits (Jiang et al., 2010; Pryce et al., 2010; Buitenhuis et al., 2014). However, the investigation of genetic mechanisms for new traits such as like cheese merit, milk composition and identification of the causal variants with pleiotropic effects have not been explored in Chinese dairy cows. The objectives of this study were to 1) detect the novel associated SNPs for five milk traits using the Mix Linear Model in Chinese Holstein cows, 2) evaluate the contribution of candidate variants with pleiotropic effects, and 3) identify candidate haplotypes for milk- related traits in Chinese Holstein cows.
2. Materials and Methods 2.1. Ethics statement All procedures strictly followed the guidelines developed by the China Council on Animal Care, and all protocols were approved by the Science Research Department of the Institute of Animal Science,
Tianjin Academy of Agricultural Science. We collected all the samples during regular quarantine inspections on the farms while all farm owners agreed to the use of these animals. 2.2. Sample genotyping and quality control Individuals (300 samples) from two farms in Tianjin were genotyped using the GeneSeek Genomic Profiler Bovine 50k SNP chip (Neogen Corporation, http://www.neogenchina.com.cn/). Genomic DNA was extracted from blood samples. The GeneSeek Genomic Profiler Bovine 50K SNP chip with 47,843 SNPs was utilized for individual genotyping. These SNPs were uniformly distributed on the genome with a mean inter-marker space of 59 kb. For quality control, individuals and their autosomal SNP were processed using PLINK v1.07 software (Purcell et al., 2007). First, all individuals with more than 10% missing SNP genotypes were excluded. Then, SNPs with either less than 95% call rate, more than 1% missing genotype, minor allele frequencies (MAF) less than 5%, and Hardy-Weinberg equilibrium (HWE) test P-value <10E-7 were also removed. 2.3. Predicted Transmitting Ability The association analysis was conducted by Predicted Transmitting Ability (PTA) which is the parent average and correction with the phenotype of offspring. It represents the expected performance of the daughters of a particular bull in relation to the average genetic herds (Cole et al., 2011; de Oliveira et al., 2014). PTA for five traits fat yield (FY), fat percentage (FP), protein percentage (PP), milk yield (MY) and cheese merit (CM) were calculated according to the standard method of CDCB (https://www.uscdcb.com/) from the U.S. Department of Agriculture (Beltsville, MD). 2.4. SNP-based GWAS using EMMAX Population stratification was assessed based on the inflation factor. Principal component analysis (PCA) was used to adjust the population stratification in GWAS for milk-related traits (Klei et al., 2008; Zhang et al., 2018). We performed GWAS for five traits (FP, FY, PP, MY and CM) using 40,501 high quality SNPs, the threshold for P-value was set to 0.05/36,532= 1.37E-06 (-log10(P) = 5.86) after Bonferroni correction. SNPs surpassing this threshold were considered as candidate variants. In this study, the mixed linear model implemented in EMMAX software was used to detect significant SNPs (Kang et al., 2008), and the equitation was as follow: y=Xβ+Zu+e
(1)
where, y is an n × 1 vector of PTA adjusted by non-genetic effects, and X is an n × q matrix of fixed
effect including SNPs and confounding variables. β is a q × 1 vector representing coefficients of the fixed effects. Z is an n × t incidence matrix relating phenotypes to the corresponding random polygenic effect. u is the random effect of the mixed model following a N (0, Iσ2g K), where K is the t × t kinship matrix. e is an n × n matrix of residual effect following a N (0, Iσ2g ). This approach can avoid repetitive variance component estimation procedure for association mapping using a mixed model (Kang et al., 2010). PCA was estimated within the population using PLINK (v1.9) (Chang et al., 2015). The statistical model for milk-related traits included the first two PCA components as covariates. Bonferroni correction for multiple tests was adopted for the value of the P-value at the genome significant of 0.05. Manhattan plot and quantile-quantile (Q-Q) plot were generated using qqman package version 0.1.4. All the statistical analyses were carried out based on R version 3.5.1 (The R Foundation, Vienna, Austria) 2.5. Candidate genes The SNPs coordinates were determined according to UMD3.1 (Zimin et al., 2009). We searched genes within 100 kb window size at upstream and downstream of associated SNPs (Table S1). In addition, we queried the significant SNPs for studied traits based on the QTLdb database (www.animalgenome.org/cgi-bin/QTLdb). Moreover, to evaluate the linkage disequilibrium (LD) pattern around the candidate regions, we estimated the parameter r2 using Haploview v4.2 software (Barrett et al., 2005). 2.6 Haplotype-based analysis The haplotype-based analysis was conducted for the target region in studied population. In this study, first, haplotype block was estimated using a standard expectation-maximum (EM) algorithm as described elsewhere (Barrett et al., 2005). Then, haplotype-based association analysis for five traits (FP, FY, PP, MY and CM) was performed using the command (--hap-assoc option) implemented in PLINK v1.07 software (Purcell et al., 2007). 3. Results 3.1. Quality control After quality control, one individual was removed due to more than 10% missing genotypes. A total of 4,523, 2,780 and 39 SNPs were deleted with missing genotyping rate > 0.01, MAF < 0.05 and HWE< 10E-7 respectively and 3969 SNPs in BTA0 and BTA30 were removed. In addition, 25 samples without PTA were removed due to incompletely pedigree information. Finally, 274 individuals with
36,532 high-quality SNPs were left for the subsequent GWAS analysis. 3.2. Population stratification assessment PCA based on the variance-standardized relationship matrix showed that the first two principal components (PC1 and PC2) contributed more than 93% of the variance. All individuals were clustered into three groups, as displayed in Fig S1, therefore, the population stratification based on PCA analysis result was considered and incorporated into the mixed linear model. 3.3. GWAS analysis In this study, 18 SNPs were significantly associated with five traits. The summary statistics of GWAS results are presented in Table 1 and the Manhattan plot for five traits is shown in Fig 1. We found the several most significant SNPs including Chr14_1757935 for FP (-log10(P) = 34.69), FY (-log10(P) = 17.57), PP (-log10(P) = 11.63), CM (-log10(P) = 7.84) and Chr14_1765835(-log10(P) = 8.54) for MY respectively. The region plot of P-values of 18 SNPs for each trait is shown in Fig 2a. For FP, we found 18 SNPs were significant associated with FP on BTA14.
Of those,
the
ARS-BFGL-NGS-94706, Chr14_1765835, BovineHD1400000246, Hapmap-24715-BTC-001-971, ARS-BFGL-NGS-26520, UA-IFASA-9288 located in the genes of VSP28, SLC52A2, MROH1, MAPK15, VAMP2, PTK2 (Table 1). For FY, nine SNPs were significantly associated with FY on BTA14.
These
SNPs
including
ARS-BFGL-NGS-
94706,
Chr14_1765835,
and
Hapmap24715-BTC-001973 were located in three genes (VPS28, SLC52A2, and MAPK15). For PP, MY and CM, we found four significant SNPs associations on BTA14. Of those, five SNPs including ARS-BFGL-NGS-57820, Chr14_1653693, Chr14_1757935, Chr14_1765835 and Chr14_2022745 were embedded or neighbored with FOXH1, ADCK5, SLC52A2, GRINA genes respectively. All the 18 significant SNPs were found embedded with discovered QTLs as reported in previous studies (Table S2). (Jiang et al., 2010; Pryce et al., 2010; Cole et al., 2011; Marques et al., 2011; Minozzi et al., 2013; Wang et al., 2013; Buitenhuis et al., 2014; Jiang et al., 2014; Capomaccio et al., 2015; Lehnert et al., 2015; Buitenhuis et al., 2016; Iso-Touru et al., 2016). These identified SNPs associated with multiple milk traits (FP, FY, MY, and PP), which may have pleiotropic effects related to milk traits in dairy cows. We found several novels QTLs for milk traits, for instance, the candidate region with 0.37 Mb (range 1.65 to 2.02 Mb) on BTA14 associated with three traits including FP, PP, CM, and ARS-BFGL-NGS-57820 was detected as a novel candidate for CM in the current study. 3.4. LD analysis of candidate region
To further determine correction between significant SNPs, we performed LD analysis in 2.5 Mb region (contains 31 SNPs) by computing pairwise r2 from Hapmap30383-BTC-005848 to UA-IFASA-9288 (1.5~4.0 Mb) on BTA14. As a result, we found two haploblocks (Fig 2b) with ~198 kb size
in
this region.
ARS-BFGL-NGS-57820,
Among them,
Chr14_1653693,
one
block had four
SNPs
BTA-34956-no-rs,
ARS-BFGL-NGS-94706,
while
the
five
SNPs
(Chr14_1757-935, Chr14_1765835, BovineHD1400000246, ARS-BFGL- NGS-71749, Chr14_202 2745) were found in other block. Moreover, we observed nine significant SNPs in candidate region with 35 candidate genes C14H8orf33, ZNF34, RPL8, COMMD5, ARHGAP39, MIR2306, C14H8orf82, LRRC24, LRRC14, RECQL4, MFSD3, GPT, PPIR16A, FOXH1, KIFC2, CYHR1, TONSL, VPS28, SLC39A4, CPSF1, ADCK5, SLC52A2, SCRT1, DGAT1, HSF1, MROH1, MIR1839, MAF1, SHARPIN, CYC1, GPAA1, EXOSC4, OPLAH, GRINA and MIR2309 respectively (Table S1). For CM and MY, regional plots of candidate region (1.33 to 2.60 Mb) on BTA14 were presented in Fig 3. The most two significant SNPs are Chr14_1757935 and Chr14_1765835 (-log10(P) = 7.84, 8.54), we also found nine SNPs with strong LD around them (CM: ARS-BFGL-NGS-57820, Chr14_1653693, Chr14_1765835, Chr14_2022745. MY: ARS-BFGL-NGS-57820, Chr14_1653693, Chr14_1757935). 3.5 Haplotype-based analysis The target region located at BTA14 (range from 1.5 to 4.0 Mb) were selected for haplotype –based analysis. Two blocks were identified with eight haplotypes across the studied population (Table S3), and these blocks were comprised of seven SNPs. Of those, the length of the large block was 190 kb, while 7.9 kb for the small block. In this study, we found three significant haplotypes (shown in Table 2) for each trait at the threshold level (P-value < 10-5). Among them, significant haplotypes H1-1 (AAGCT), H2-1 (TG), H2-3 (CC) were identified for FP, FY, MY, PP, while the most significant haplotype was found for FP trait (-log10(P) =21.88, 26.09,28.19). Notably, three haplotypes including H2-1 (-log10(P) = 4.29) and H2-3 (-log10(P) = 4.25), and H1-5 (AGGTC, -log10(P) = 4.25) showed significant association for CM. 4. Discussion In this study, we performed a GWAS for five important milk-related traits (fat percentage, fat yield, protein percentage, milk yield and cheese merit) using GeneSeek Genomic Profiler Bovine 50K SNP chip in Chinese Holstein cows. Some of these traits have been explored in previous studies; however, elucidation of the molecular mechanism underpinning milk traits in other populations can offer
additional insights for understanding the genetic basis of these traits in dairy cows. In the present study, we first estimated the population stratification in our population. Using PCA analysis implemented in PLINK v1.9 (Chang et al., 2015), we observed our data can be clustered into three groups, and the inflation factor (λ) suggested obvious population stratification for the studied population. Therefore, population stratification and SNP effects were considered as fixed effects, and the polygenic effects were considered as random factors in mix linear model (Banerjee et al., 2008; Ferreira and Purcell, 2009; Kim and Xing, 2009; O'Reilly et al., 2012). These methods not only can increase power for detecting pleiotropic genetic variants but also avoid false-positive result in association tests (Stephens, 2013). Previous studies found several regions and genes on BTA14 were widely associated with many important traits based on multiple statistical association tests. For instance, DGAT1 located at ~1.8Mb on BTA14 is generally accepted as a major gene affecting milk production traits (milk yield, protein content, fat content, and fatty acid composition) (Grisart et al., 2002; Schennink et al., 2007; Schennink et al., 2008; Conte et al., 2010). CYP11B1 located at ~1.33 Mb has significant effects on milk yield, protein yield, fat percentage, protein percentage (Kaupe et al., 2007). In this study, we detected 18 significant SNPs for milk-related traits on BTA14. Based these SNPs, we obtained 29 candidate genes including C14H8orf33, FOXH1, CYHR1, VPS28, SCRT1, DGAT1, OPLAH, MROH1, MAF1, GRINA, LY6D, LY6H, ZNF34, RPL8, GPT, PPP1R16A, CPSF1, HSF1, GPAA1, PUF60, MAPK15, EEF1D, LYNX1, PTK2, ARHGAP39, TONSL, LRRC14, ADCK5, which have been reported in previous studies. Among them, 22 genes (C14H8orf33, FOXH1, CYHR1, VPS28, DGAT1, OPLAH, MAF1, GRINA, LY6D, LY6H, ZNF34, RPL8, GPT, PPP1R16A, CPSF1, HSF1, GPAA1, PUF60, MAPK15, EEF1D, LYNX1, PTK2) have been detected for multiple milk traits in Chinese dairy cows (Jiang et al., 2010; Jiang et al., 2014). ARHGAP39 was identified as a candidate gene for mastitis susceptibility in Holsteins (Wang et al., 2015). TONSL, LRRC14, ADCK5 were detected as potential candidate genes for milk traits or mammary gland functions (Ibeagha-Awemu et al., 2016). MROH1 and SCRT1 were affected by FP (Capomaccio et al., 2015). In our study, we found several new candidate genes for milk-related traits. For instance, CPSF1, DGAT1, GPT, PPP1R16A, and CYHR1 were associated with CM. HSF1 was associated with FY, PP, MY, CM. ZNF34 was associated with FY, ARHGAP39, TSTA3 for FP and GRINA for FY and CM (Table S1). We also found several new candidate genes including SLC52A2, VAMP2 embedded with
three SNPs (Chr14_1765835 and ARS-BFGL-NGS-26520). Of these, SLC52A2 was associated with FP, FY, PP, MY, CM and VAMP2 was associated with FP (Table 1). In addition, we observed many significant SNPs (BTA-34956-no-rs, ARS-BFGL-NGS-57820, Chr14_1653693, Chr14_1765835, ARS-BFGL-NGS-71749,
Hapmap24715-BTC-001973,
ARS-BFGL-NGS-26520,
ARS-BFGL-NGS-22866) associated with milk traits nearside several candidate genes (Table S1). Our study identified 18 significant SNPs region located at ~2.5 Mb (from 1.65 to 2.02 Mb) on BTA14 associated with multi-traits. Of these, we found several new variants (ARS-BFGL-NGS-57820, Chr14_1653693, Chr14_1757935, Chr14_1765835, Chr14_2022745) for FP, PP and CM located within a 0.37Mb (between 1.65 to 2.02 Mb) on BTA14. Previous studies suggested many QTLs with pleiotropic effects related to milk-related traits in dairy cows (Olsen et al., 2011; Bolormaa et al., 2014; Saatchi et al., 2014). Our study identified 18 significant SNPs region located at ~2.5 Mb (from 1.65 to 2.02 Mb) on BTA14 associated with multi-traits. Of these, we found several new variants (ARS-BFGL-NGS-57820, Chr14_1653693, Chr14_1757935, Chr14_1765835, Chr14_2022745) for FP, PP and CM located within a 0.37Mb (between 1.65 to 2.02 Mb) on BTA14. Moreover, Haplotype-base analysis further confirmed the results of our SNP-based GWAS, we found three significant haplotypes for multi traits, and these haplotypes contain the identified SNPs, which have been detected by SNP-based GWAS. Haplotype H2-1, H2-3 were associated with the five traits (FP, FY, PP, MY, CM) and H1-1 was associated with four traits (FP, FY, PP, MY). The identification of candidate SNPs and genes can help improve our knowledge of milk traits and provide useful information for guarding the design of breeding programs. However, more precise genomic analysis and validation of functional variations with larger sample size are required in further study. 5. Conclusions Our study identified eighteen candidate SNPs for five traits in Chinese Holstein cows, and five of them were novel discovered in the current study. These findings suggested that nine candidate SNPs may have pleiotropic effects for milk-related traits. Our results may further help to promote the design of marker-assisted selection for a breeding program. But this candidate genes must be a further
investigated for the improvement of milk traits.
Author Contributions: Conceptualization, Wenqin Song, Lingyang Xu and Yi Ma; Data curation, Tianzhen Wang,
Chengbin Chen, Dawei Yao and Jing Ma; Investigation, Tianzhen Wang and Jiao Li; Methodology, Xue Gao and Lingyang Xu; Project administration, Yi Ma; Resources, Yi Ma; Writing – original draft, Tianzhen Wang; Writing – review & editing, Tianzhen Wang and Lingyang Xu.
Funding: This research was funded by the Innovation Team of Tianjin Cattle Research System (ITTCRS2017005) and Tianjin Seed Industry Special Project (16ZXZYNC00220).
Acknowledgments: The authors would like to thank the breeders and breeding associations who provide us samples. The SNPs data and PTA data were provided by Neogen Corporation ([email protected])
Declarations Ethics and consent to participate All animals used in the study were treated following the guidelines for the experimental animals established by the Council of China Animal Welfare. Animal experiments were approved by the Science Research Department of the Institute of Animal Science, Tianjin Academy of Agricultural Science. All farm owners agreed on the use of the animals.
Consent to publish Not applicable
Competing interests The authors declare that they have no competing interests. Availability of data The data and compute programs used in this manuscript are available from the corresponding author on request.
Figure legends
Figure 1. a) Manhattan plot of -log10(P) for FY in Chinese Holsteins dairy cows. Chromosomes are color-coded separately. The horizontal reference line indicated the genome-wide significance
levels (-log10(P) = 5.98). b) Q-Q plots for FY. Deviations from the slope lines correspond to loci that deviate from the null hypothesis. c) Manhattan plot of -log10(P) for FP. d) Q-Q plots for FP. e) Manhattan plot of -log10(P) for MY. f) Q-Q plots for five MY. g) Manhattan plot of -log10(P) for PP. h) Q-Q plots for PP. i) Manhattan plot of -log10(P) for CM. j) Q-Q plots for CM.
Figure 2. a) Association test with -log10(P) for 31 SNPs for five traits including FP, FY, PP, MY and CM. The horizontal reference line indicated the genome-wide significance levels (-log10(P) = 5.98), 18 of 31 were significant. b) Haplotype analysis for the region around BTA14: 1489496– 3956956 bp. Two haplotypes with 182 and 264 kb were estimated using Haploview software.
Figure 3. a) Regional plots of candidate region at 1.33 to 2.60 Mb on BTA14 for CM. b) Regional plots of candidate region at 1.33 to 2.60 Mb on BTA14 for MY.
All authors agree to publish.
Conflict of interest The authors declare that they have no competing interests.
References Archer-Nicholls, S., Lowe, D., Utembe, S., Allan, J., Zaveri, R.A., Fast, J.D., Hodnebrog, O., van der Gon, H.D., McFiggans, G., 2014. Gaseous chemistry and aerosol mechanism developments for version 3.5.1 of the online regional model, WRF-Chem. Geoscientific Model Development 7, 2557-2579. Banerjee, S., Yandell, B.S., Yi, N., 2008. Bayesian quantitative trait loci mapping for multiple traits. Genetics 179, 2275-2289. Barrett, J.C., Fry, B., Maller, J., Daly, M.J., 2005. Haploview: analysis and visualization of LD and haplotype maps. Bioinformatics 21, 263-265. Bolormaa, S., Pryce, J.E., Reverter, A., Zhang, Y., Barendse, W., Kemper, K., Tier, B., Savin, K., Hayes, B.J., Goddard, M.E., 2014. A multi-trait, meta-analysis for detecting pleiotropic polymorphisms for stature, fatness and reproduction in beef cattle. PLoS Genetics 10, 3. Buitenhuis, B., Janss, L.L., Poulsen, N.A., Larsen, L.B., Larsen, M.K., Sorensen, P., 2014. Genome-wide association and biological pathway analysis for milk-fat composition in Danish Holstein and Danish Jersey cattle. BMC Genomics 15, 1112. Buitenhuis, B., Poulsen, N.A., Gebreyesus, G., Larsen, L.B., 2016. Estimation of genetic parameters and detection of chromosomal regions affecting the major milk proteins and their post translational modifications in Danish Holstein and Danish Jersey cattle. BMC Genetics 17, 114. Cai, Z., Guldbrandtsen, B., Lund, M.S., Sahana, G., 2018. Prioritizing candidate genes post-GWAS using multiple sources of data for mastitis resistance in dairy cattle. BMC Genomics 19, 656. Capomaccio, S., Milanesi, M., Bomba, L., Cappelli, K., Nicolazzi, E.L., Williams, J.L., Ajmone-Marsan, P., Stefanon, B., 2015. Searching new signals for production traits through gene-based association analysis in three Italian cattle breeds. Animal Genetics 46, 361-370. Chang, C.C., Chow, C.C., Tellier, L.C., Vattikuti, S., Purcell, S.M., Lee, J.J., 2015. Second-generation PLINK: rising to the challenge of larger and richer datasets. Gigascience 4, 7. Cole, J.B., Wiggans, G.R., Ma, L., Sonstegard, T.S., Lawlor, T.J., Jr., Crooker, B.A., Van Tassell, C.P., Yang, J., Wang, S., Matukumalli, L.K., Da, Y., 2011. Genome-wide association analysis of thirty one production, health, reproduction and body conformation traits in contemporary U.S. Holstein cows. BMC Genomics 12, 408. Conte, G., Mele, M., Chessa, S., Castiglioni, B., Serra, A., Pagnacco, G., Secchiari, P., 2010. Diacylglycerol acyltransferase 1, stearoyl-CoA desaturase 1, and sterol regulatory element binding protein 1 gene polymorphisms and milk fatty acid composition in Italian Brown cattle. Journal of Dairy
Science 93, 753-763. Cordell, H.J., 2009. Detecting gene-gene interactions that underlie human diseases. Nature Reviews. Genetics 10, 392-404. de Oliveira, F.C., Borges, C.C., Almeida, F.N., e Silva, F.F., da Silva Verneque, R., da Silva, M.V., Arbex, W., 2014. SNPs selection using support vector regression and genetic algorithms in GWAS. BMC Genomics 15(7), 1-15. Do, D.N., Bissonnette, N., Lacasse, P., Miglior, F., Sargolzaei, M., Zhao, X., Ibeagha-Awemu, E.M., 2017. Genome-wide association analysis and pathways enrichment for lactation persistency in Canadian Holstein cattle. Journal of Dairy Science 100, 1955-1970. Fan, Y., Wang, P., Fu, W., Dong, T., Qi, C., Liu, L., Guo, G., Li, C., Cui, X., Zhang, S., Zhang, Q., Zhang, Y., Sun, D., 2014. Genome-wide association study for pigmentation traits in Chinese Holstein population. Animal Genetics 45, 740-744. Fang, L., Sahana, G., Su, G., Yu, Y., Zhang, S., Lund, M.S., Sorensen, P., 2017. Integrating Sequence-based GWAS and RNA-Seq Provides Novel Insights into the Genetic Basis of Mastitis and Milk Production in Dairy Cattle. Scientific Reports 7, 45560. Ferreira, M.A., Purcell, S.M., 2009. A multivariate test of association. Bioinformatics 25, 132-133. Gao, Y., Jiang, J., Yang, S., Hou, Y., Liu, G.E., Zhang, S., Zhang, Q., Sun, D., 2017. CNV discovery for milk composition traits in dairy cattle using whole genome resequencing. BMC Genomics 18, 265. Goddard, M.E., Hayes, B.J., 2009. Mapping genes for complex traits in domestic animals and their use in breeding programmes. Nature reviews. Genetics 10, 381-391. Grisart, B., Coppieters, W., Farnir, F., Karim, L., Ford, C., Berzi, P., Cambisano, N., Mni, M., Reid, S., Simon, P., Spelman, R., Georges, M., Snell, R., 2002. Positional candidate cloning of a QTL in dairy cattle: identification of a missense mutation in the bovine DGAT1 gene with major effect on milk yield and composition. Genome Research 12, 222-231. Ibeagha-Awemu, E.M., Peters, S.O., Akwanji, K.A., Imumorin, I.G., Zhao, X., 2016. High density genome- wide genotyping-by-sequencing and association identifies common and low-frequency SNPs, and novel candidate genes influencing cow milk traits. Scientific Reports 6. Iso-Touru, T., Sahana, G., Guldbrandtsen, B., Lund, M.S., Vilkki, J., 2016. Genome-wide association analysis of milk yield traits in Nordic Red Cattle using imputed whole genome sequence variants. BMC Genetics 17, 55. Jiang, L., Liu, J., Sun, D., Ma, P., Ding, X., Yu, Y., Zhang, Q., 2010. Genome wide association studies for milk production traits in Chinese Holstein population. PloS One 5, 10. Jiang, L., Liu, X., Yang, J., Wang, H., Jiang, J., Liu, L., He, S., Ding, X., Liu, J., Zhang, Q., 2014. Targeted resequencing of GWAS loci reveals novel genetic variants for milk production traits. BMC Genomics 15, 1105. Kang, H.M., Sul, J.H., Service, S.K., Zaitlen, N.A., Kong, S.Y., Freimer, N.B., Sabatti, C., Eskin, E., 2010. Variance component model to account for sample structure in genome-wide association studies. Nature Genetics 42, 348-354. Kang, H.M., Zaitlen, N.A., Wade, C.M., Kirby, A., Heckerman, D., Daly, M.J., Eskin, E., 2008. Efficient control of population structure in model organism association mapping. Genetics 178, 1709-1723. Kaupe, B., Brandt, H., Prinzenberg, E.M., Erhardt, G., 2007. Joint analysis of the influence of CYP11B1 and DGAT1 genetic variation on milk production, somatic cell score, conformation, reproduction, and productive lifespan in German Holstein cattle. Journal of Animal Science 85, 11-21. Kim, S., Xing, E.P., 2009. Statistical estimation of correlated genome associations to a quantitative trait
network. PLoS Genetics 5, 8. Klei, L., Luca, D., Devlin, B., Roeder, K., 2008. Pleiotropy and principal components of heritability combine to increase power for association analysis. Genetic Epidemiology 32, 9-19. Lehnert, K., Ward, H., Berry, S.D., Ankersmit-Udy, A., Burrett, A., Beattie, E.M., Thomas, N.L., Harris, B., Ford, C.A., Browning, S.R., Rawson, P., Verkerk, G.A., van der Does, Y., Adams, L.F., Davis, S.R., Jordan, T.W., MacGibbon, A.K., Spelman, R.J., Snell, R.G., 2015. Phenotypic population screen identifies a new mutation in bovine DGAT1 responsible for unsaturated milk fat. Scientific Reports 5, 8484. Li, X., Buitenhuis, A.J., Lund, M.S., Li, C., Sun, D., Zhang, Q., Poulsen, N.A., Su, G., 2015. Joint genome-wide association study for milk fatty acid traits in Chinese and Danish Holstein populations. Journal of Dairy Science 98, 8152-8163. Liu, A., Wang, Y., Sahana, G., Zhang, Q., Liu, L., Lund, M.S., Su, G., 2017. Genome-wide Association Studies for Female Fertility Traits in Chinese and Nordic Holsteins. Scientific Reports 7, 8487. Marques, E., Grant, J.R., Wang, Z., Kolbehdari, D., Stothard, P., Plastow, G., Moore, S.S., 2011. Identification of candidate markers on bovine chromosome 14 (BTA14) under milk production trait quantitative trait loci in Holstein. Journal of animal breeding and genetics = Zeitschrift fur Tierzuchtung und Zuchtungsbiologie 128, 305-313. Maxa, J., Neuditschko, M., Russ, I., Forster, M., Medugorac, I., 2012. Genome-wide association mapping of milk production traits in Braunvieh cattle. Journal of Dairy Science 95, 5357-5364. Meuwissen, T.H., Hayes, B.J., Goddard, M.E., 2001. Prediction of total genetic value using genome-wide dense marker maps. Genetics 157, 1819-1829. Minozzi, G., Nicolazzi, E.L., Stella, A., Biffani, S., Negrini, R., Lazzari, B., Ajmone-Marsan, P., Williams, J.L., 2013. Genome wide analysis of fertility and production traits in Italian Holstein cattle. PloS One 8, 1-10. Nayeri, S., Sargolzaei, M., Abo-Ismail, M.K., May, N., Miller, S.P., Schenkel, F., Moore, S.S., Stothard, P., 2016. Genome-wide association for milk production and female fertility traits in Canadian dairy Holstein cattle. BMC Genetics 17, 75. Ning, C., Wang, D., Zheng, X., Zhang, Q., Zhang, S., Mrode, R., Liu, J.F., 2018. Eigen decomposition expedites longitudinal genome-wide association studies for milk production traits in Chinese Holstein. Genetics, Selection, Evolution : GSE 50, 12. O'Reilly, P.F., Hoggart, C.J., Pomyen, Y., Calboli, F.C., Elliott, P., Jarvelin, M.R., Coin, L.J., 2012. MultiPhen: joint model of multiple phenotypes can increase discovery in GWAS. PloS One 7, e34861. Olsen, H.G., Hayes, B.J., Kent, M.P., Nome, T., Svendsen, M., Larsgard, A.G., Lien, S., 2011. Genome-wide association mapping in Norwegian Red cattle identifies quantitative trait loci for fertility and milk production on BTA12. Animal Genetics 42, 466-474. Pausch, H., Emmerling, R., Gredler-Grandl, B., Fries, R., Daetwyler, H.D., Goddard, M.E., 2017. Meta-analysis of sequence-based association studies across three cattle breeds reveals 25 QTL for fat and protein percentages in milk at nucleotide resolution. BMC Genomics 18, 853. Pegolo, S., Dadousis, C., Mach, N., Ramayo-Caldas, Y., Mele, M., Conte, G., Schiavon, S., Bittante, G., Cecchinato, A., 2017. SNP co-association and network analyses identify E2F3, KDM5A and BACH2 as key regulators of the bovine milk fatty acid profile. Scientific Reports 7, 17317. Pryce, J.E., Bolormaa, S., Chamberlain, A.J., Bowman, P.J., Savin, K., Goddard, M.E., Hayes, B.J., 2010. A validated genome-wide association study in 2 dairy cattle breeds for milk production and fertility traits using variable length haplotypes. Journal of Dairy Science 93, 3331-3345. Purcell, S., Neale, B., Todd-Brown, K., Thomas, L., Ferreira, M.A., Bender, D., Maller, J., Sklar, P., de
Bakker, P.I., Daly, M.J., Sham, P.C., 2007. PLINK: a tool set for whole-genome association and population-based linkage analyses. American Journal of Human Genetics 81, 559-575. Saatchi, M., Schnabel, R.D., Taylor, J.F., Garrick, D.J., 2014. Large-effect pleiotropic or closely linked QTL segregate within and across ten US cattle breeds. BMC Genomics 15, 442. Schennink, A., Heck, J.M., Bovenhuis, H., Visker, M.H., van Valenberg, H.J., van Arendonk, J.A., 2008. Milk fatty acid unsaturation: genetic parameters and effects of stearoyl-CoA desaturase (SCD1) and acyl CoA: diacylglycerol acyltransferase 1 (DGAT1). Journal of Dairy Science 91, 2135-2143. Schennink, A., Stoop, W.M., Visker, M.H., Heck, J.M., Bovenhuis, H., van der Poel, J.J., van Valenberg, H.J., van Arendonk, J.A., 2007. DGAT1 underlies large genetic variation in milk-fat composition of dairy cows. Animal Genetics 38, 467-473. Schopen, G.C., Visker, M.H., Koks, P.D., Mullaart, E., van Arendonk, J.A., Bovenhuis, H., 2011. Whole-genome association study for milk protein composition in dairy cattle. Journal of Dairy Science 94, 3148-3158. Stephens, M., 2013. A unified framework for association analysis with multiple related phenotypes. PloS One 8, e65245. Wang, H., Jiang, L., Liu, X., Yang, J., Wei, J., Xu, J., Zhang, Q., Liu, J.F., 2013. A post-GWAS replication study confirming the PTK2 gene associated with milk production traits in Chinese Holstein. PloS One 8, e83625. Wang, X., Ma, P.P., Liu, J.F., Zhang, Q., Zhang, Y., Ding, X.D., Jiang, L., Wang, Y.C., Zhang, Y., Sun, D.X., Zhang, S.L., Su, G.S., Yu, Y., 2015. Genome-wide association study in Chinese Holstein cows reveal two candidate genes for somatic cell score as an indicator for mastitis susceptibility. BMC Genetics 16. Zhang, W., Gao, X., Shi, X., Zhu, B., Wang, Z., Gao, H., Xu, L., Zhang, L., Li, J., Chen, Y., 2018. PCA-based multiple-trait GWAS analysis: A powerful model for exploring pleiotropy. Animal 8, 12. Zhang, X., Chu, Q., Guo, G., Dong, G., Li, X., Zhang, Q., Zhang, S., Zhang, Z., Wang, Y., 2017. Genome-wide association studies identified multiple genetic loci for body size at four growth stages in Chinese Holstein cattle. PloS One 12, e0175971. Zimin, A.V., Delcher, A.L., Florea, L., Kelley, D.R., Schatz, M.C., Puiu, D., Hanrahan, F., Pertea, G., Van Tassell, C.P., Sonstegard, T.S., Marcais, G., Roberts, M., Subramanian, P., Yorke, J.A., Salzberg, S.L., 2009. A whole-genome assembly of the domestic cow, Bos taurus. Genome Biology 10(4), 1-10.
Table 1. Candidate SNPs list associated with milk traits fat percentage (FP), fat yield (FY), protein percentage (PP), milk yield (MY), cheese merit (CM)) on BTA14. SNP
Position
Gene
Distance
MAF
-log10(P)
Traits
Hapmap30383-B
1489496
C14H8orf33
78
0.3571
14.55
FP
BTA-34956-no-rs
1514056
RPL8
-6430
0.4891
6.75
FP
ARS-BFGL-NGS
1651311
FOXH1
3390
0.2719
27.13
FP
FY
PP
MY
CM
Chr14_1653693
1653693
FOXH1
1008
0.2646
27.76
FP
FY
PP
MY
CM
ARS-BFGL-NGS
1696470
VPS28
in
0.4088
13.77
FP
FY
Chr14_1757935
1757935
ADCK5
-1634
0.2527
33.49
FP
FY
PP
MY
CM
Chr14_1765835
1765835
SLC52A2
in
0.2792
32.55
FP
FY
PP
MY
CM
BovineHD140000
1880378
MROH1
in
0.4487
8.61
FP
1954317
OPLAH
3427
0.1752
6.94
FP
FY
PP
FY
TC-005848
-57820
-94706
0246 ARS-BFGL-NGS -71749 Chr14_2022745
2022745
GRINA
-1036
0.2719
28.36
FP
Hapmap25384-B
2217163
MAPK15
17870
0.4744
9.02
FP
2239085
MAPK15
in
0.4707
9.52
FP
BTA-35941-no-rs
2276443
MIR193A-2
4374
0.4635
10.23
FP
ARS-BFGL-NGS
2386688
VAMP2
in
0.2482
8.92
FP
2524432
LY6H
38728
0.4234
11.16
FP
2553525
LY6H
9635
0.3704
9.85
FP
2826632
LYNX1
-9946
0.3485
6.18
FP
3956956
PTK2
in
0.293
7.68
FP
TC-001997 Hapmap24715-B
FY
TC-001973
-26520 Hapmap30086-B TC-002066 Hapmap30646-B TC-002054 ARS-BFGL-NGS -22866 UA-IFASA-9288
MAF: minor allele frequencies
FY
CM
Table 2. Significant haplotypes from the haplotype-based analysis. H1-1, H1-5 is the first block with Hapmap30381-BTC-005750,
Hapmap30383-BTC-005848,
BTA-34956-no-rs,
ARS-BFGL-NGS-57820,
Chr14_1653693, and H2-1, H2-3 is the second block with Chr14_1757935, Chr14_1765835 Trait
PP
MY
FY
FP
CM
Hapname
Haplotype
Start
End
Length/bp.
P-value
H1-1
AAGCT
1463676
1653693
190017
1.18E-08
H2-1
TG
1757935
1765835
7900
1.41E-09
H2-3
CC
1757935
1765835
7900
2.56E-10
H1-1
AAGCT
1463676
1653693
190017
4.34E-06
H2-1
TG
1757935
1765835
7900
2.25E-06
H2-3
CC
1757935
1765835
7900
1.11E-07
H1-1
AAGCT
1463676
1653693
190017
6.73E-08
H2-1
TG
1757935
1765835
7900
8.13E-10
H2-3
CC
1757935
1765835
7900
3.21E-09
H1-1
AAGCT
1463676
1653693
190017
1.31E-22
H2-1
TG
1757935
1765835
7900
8.06E-27
H2-3
CC
1757935
1765835
7900
6.47E-29
H1-5
AGGTC
1463676
1653693
190017
5.66E-05
H2-1
TG
1757935
1765835
7900
5.16E-05
H2-3
CC
1757935
1765835
7900
5.64E-05