Genetic correlation and genome-wide association study (GWAS) of the length of productive life, days open, and 305-days milk yield in crossbred Holstein dairy cattle

Abstract

In this study, we estimated the genetic parameters and identified the putative quantitative trait loci (QTL) associated with the length of productive life (LPL), days open (DO), and 305-day milk yield for the first lactation (FM305) of crossbred Holstein dairy cattle. Data comprising 4,739 records collected between 1986 and 2004 were used to estimate the variance-covariance components using the multiple-trait animal linear mixed models based on the average information restricted maximum likelihood (AI-REML) algorithm. Thirty-six animals were genotyped using the Illumina BovineSNP50 Bead Chip [>50,000 single nucleotide polymorphisms (SNPs)] to identify the putative QTL in a genome-wide association study. The heritability of the production trait FM305 was 0.25 and that of the functional traits, LPL and DO, was low (0.10 and 0.06, respectively). The genetic correlation estimates demonstrated favorable negative correlations between LPL and DO (−0.02). However, we observed a favorable positive correlation between FM305 and LPL (0.43) and an unfavorable positive correlation between FM305 and DO (0.1). The GWAS results indicated that 23 QTLs on bovine chromosomes 1, 4, 5, 8, 15, 26, and X were associated with the traits of interest, and the putative QTL regions were identified within seven genes (SYT1, DOCK11, KLHL13, IL13RA1, PRKG1, GNA14, and LRRC4C). In conclusion, the heritability estimates of the LPL and DO were low. Therefore, the approach of multiple-trait selection indexes should be applied, and the QTL identified here should be considered for use in marker-assisted selection in the future.

Introduction

Functional traits, i.e., the length of productive life, health, and fertility, have become increasingly important to the dairy cattle industry. They have a potential effect on the cost of milk production (Groen et al., 1997). Consequently, the breeding objectives for dairy cattle have changed drastically in several countries in the last decade (Veerkamp et al., 2013). Various production and functional traits have been integrated into selection indexes in each country (Miglior et al., 2005). In this study, we investigated functional traits [length of productive life (LPL) and days open (DO)], and a production trait [305-days milk yield of the first lactation (FM305)]. Several studies have reported that high milk yield had negative effects on LPL and fertility (Pryce et al., 2004). The high-yielding Holstein cows have poor health and fertility and higher risk of culling. The LPL of high-yielding cows decreased from 3.4 to 2.8 (Knaus, 2009). Van Raden et al. (2004) developed a national genetic evaluation for cow fertility in the United State of America. They used a multiple-trait linear animal model, which included DO, LPL, milk, fat, protein yield, and somatic cell score. The genetic correlation between LPL and milk yield was found to be 0.03. This positive genetic correlation was in agreement with the results of Irano et al. (2014), who studied the genetic association between three economically important traits - milk yield, stayability, and mastitis - in Holstein cows under tropical conditions using the threshold linear multiple-trait animal model. Therefore, it is important to investigate the genetic parameters of traits before combining them in the selection index.

Recently, a genome-wide association study (GWAS) approach that utilizes all genotypes, phenotypes, and pedigree information jointly in one step has been proposed by Wang et al. (2012). It is also referred to as the single-step GWAS (ssGWAS). It can identify the quantitative trait loci (QTL) regions and candidate genes that affect the traits of interest. In dairy cattle, the GWAS of production traits, such as milk yield, fat yield, and protein yield, (Maxa et al., 2012; Minozzi et al., 2013; Nayeri et al., 2016; van den Berg et al., 2016); functional or fertility traits, such as calving to the first service interval, days open, non-return rate, and conception rate (Olsen et al., 2011; Minozzi et al., 2013; Nayeri et al., 2016; Parker Gaddis et al., 2016); health traits, such as clinical mastitis and udder type (Flury et al., 2014; Sahana et al., 2014); and longevity (Zhang et al., 2016), have been reported. A previous GWAS showed that the DGAT1 gene on bovine chromosome BTA14 was associated with milk production and the FAM1818 gene on BTA21 was associated with days open (Minozzi et al., 2013; Nayeri et al., 2016). Moreover, Kühn et al. (2003) reported the effect of QTL region on the length of productive life and the effect of the SIGLEC12 gene on BTA18 on the longevity of the German Holstein population. However, many studies focused on the milk production and fertility trait and only a few have examined longevity or the length of productive life. In addition, those results are studied in purebred dairy cattle, which might not be applicable to crossbred Holstein cattle. Therefore, the aims of this study were to: 1) estimate the genetic parameters of the LPL, DO, and FM305 in crossbred dairy cattle for their selection in terms of genetic performance, and 2) investigate the putative QTL associated to LPL, DO, and FM305 using GWAS.

Materials and Methods

Animals and phenotypes

The raw data were collected from the crossbred dairy cows, which were born between 1986 and 2004 and raised under hot-humid environmental conditions in a commercial farm management in Thailand. After data editing, all available functional and production traits collected from 4,739 records were obtained for the calculation of LPL (interval from first calving to culling), DO (interval from calving to conception), and FM305 traits. Data from a pedigree file with a total of 6,734 animals, both with and without record, were included.

The LPL was defined as the number of days or months from the first calving to the last calving date of the last known lactation or culling date (Ducrocq, 1994; Figure 1).

Genotype data

A total of 36 animals were genotyped using the Illumina BovineSNP50 Bead Chip (Illumina Inc., San Diego, CA, USA), which contains probes to test more than 50,000 single nucleotide polymorphisms (SNPs). A commercial laboratory performed the SNP genotyping. The following quality control (QC) criteria were used for excluding SNPs from the data set: 1) monomorphic SNPs, 2) SNPs with minor allele frequency (MAF) less than 5%, 3) SNPs with call rate less than 90%, 4) animal genotypes with call rate less than 90%, and 5) parent-progeny Mendelian conflicts. Finally, after genotype quality control, a total 34 individuals and 43,218 SNPs were selected for the GWAS analysis.

Statistical analyses

The following multiple-trait animal linear mixed model was applied in this study:

y = Xβ + Za +e (Equation 1)

where y = vector of observations for LPL, DO, and FM305; β = vector of fixed effects for an interaction between month-calving and year-calving with 205 levels, Holstein genetics fractions effect with three levels (HFG1 < 87.50%, HFG2 ≥ 87.50 and < 93.75%, and HFG3 ≥ 93.75 and < 100%), and age at calving effect covariate; a = vector of additive genetic effects; ε = vectors of random residual effects; and X and Z = incidence matrices for the corresponding effects.

Genetic parameters estimation

The data were analyzed by using the multiple-trait animal linear mixed models based on the average information-restricted maximum likelihood (AI-REML) algorithm (Johnson and Thompson, 1995). The variance-covariance components (VCE) of all three traits were estimated using AIREMLF90 (Misztal et al., 2002). The variance-covariance structure of the random effects can be described as follows:

(Equation 2)

where A = numerator relationship matrix; = additive genetic variance; I = identity matrix; and = random residual variance.

Genome-wide association analysis (Single-step method)

Misztal et al. (2009) demonstrated that a numerator relationship matrix (A) can be modified to a matrix (H) that includes both pedigree and genomic relationships as shown below:

(Equation 3)

The inverse H matrix could be derived from the following equation (Aguilar et al., 2010):

(Equation 4)

where $A_{22}^{-1}$ is the inverse of a pedigree for genotyped animals. G is the genomic relationship matrix (Legarra et al., 2009).

G = ZDZ′(Equation 5)

where D is a diagonal matrix with elements containing the inverse of the expected marker (Van Raden, 2008)

A single-step genome-wide association study (ssGWAS) was described by Wang et al. (2012). The GEBV (genomic estimated breeding values) solutions were used to estimate the marker effects through an iterative process. A detailed description of the iterative algorithm was outlined in Wang et al. (2012).

The equation for predicting SNP effects uses the weighted genomic relationship matrix (Wang et al., 2012):

(Equation 6)

Where $\widehat{\mu }$ is a vector of SNP marker effects, D is a diagonal matrix of the weights for variances of SNP effects, Z is a matrix relating genotypes of each locus and ${\widehat{a}}_g$ is the animal effect of genotyped animals.

The equation to estimate individual variance of SNP effect (Zhang et al., 2010):

(Equation 7)

Where ${\widehat{\sigma }}_{u,i}^2$ is the genetic additive variance of each SNP marker, ${\widehat{\mu }}_i^2$ is the square of the ith SNP marker effect, P_i is the allele frequency of the second allele of the i^th marker in the current population.

Gene search and identification of QTL

The variance of 10-SNP windows was computed for each individual. For this, (43,218-(9*30)/10) = 4,319 10-SNP windows were tested in the whole cattle genome (the method followed Schneider et al., 2012). Therefore, the expected proportion of variance accounted for by 1 window was 2.3 E^-04 (1/4,319). A QTL is a combination of consecutive 10-SNP windows with a greater proportion of genetic variance than the expected proportion of variance accounted for by 1 window (2.3 E^-04). Subsequently, we used the National Center for Biotechnology Information (NCBI) database to pinpoint possible causative genes that might relate to the functional and production traits in cattle.

Results and Discussion

The summary of the basic statistics for each trait (LPL, DO, and FM305) are presented in Table 1.

Table 1: Data summary for the length of productive life (LPL), days open (DO), and 305-days milk yield (M305) traits.

Estimation of genetic parameters

The heritability and genetic correlation estimates for LPL, DO, and FM305 are presented in Table 2.

Table 2: Estimates of the heritability (bold face, diagonal), additive genetic correlation (off the diagonal), and their standard errors for LPL, DO, and FM305 trait, respectively.

The heritability estimates obtained in this study revealed that the heritability of the production trait FM305 was higher (0.25) than that of the other two functional traits. The heritability estimates of both LPL (h² = 0.10) and DO (h² = 0.06) were in the range as reported by Van Raden et al. (2004) and Eghbalsaied (2011). However, the milk heritability was slightly lower than the values reported by Irano et al. (2014), Tsuruta et al. (2005), Pritchard et al. (2013), which were 0.28, 0.4, and 0.3, respectively. The heritability of LPL was also slightly lower than the value (0.13) reported by Wasana et al. (2015).

The genetic correlation estimates showed a favorable negative correlation between LPL and DO, and found to be low ( $r_{g,LPL,DO}$ = -0.04). A favorable positive correlation was observed between the production trait FM305 and functional trait LPL ( $r_{g,LPL,FM305}$ = 0.43) mainly because the dairy cattle producers prefer to keep high-milk-yielding cattle in the herd. The positive genetic correlation between FM305 and LPL was also reported by Irano et al. (2014) and Wasana et al. (2015). Our results also demonstrate an unfavorable positive correlation between FM305 and DO ( $r_{g,DO,FM305}$ = 0.14), indicating that the genetic selection for high milk yield could reduce the fertility of dairy cattle as previously pointed out by Van Raden et al. (2004).

Investigation of QTL using GWAS

In this study, the GWAS results for all traits with variances explained for the 10-SNP windows used to identify the putative QTL regions. The QTL regions were mapped by a high proportion of variances (higher than 2.3E-04, which is an expected proportion of the variance accounted for by one window). The diffuse peak spread mainly on chromosomes 1, 4, 5, 8, 15, 26, and X (Figure 2A to 2C).

A total of 23 QTL regions were associated with all three traits, and the genes located within the QTL regions are shown in Table 3. We detected 10 SNPs in the QTL regions associated with LPL, and they were mostly located on chromosomes BTA5 and BTAX. Four genes were located within the QTL regions on BTA5 (SYT1) and BTAX (DOCK11, KLHL13, and IL13RA1). Eight QTLs associated with DO were located on chromosomes BTA1, BTA4, and BTA26; only one gene (PRKG1) was located within the QTL region on chromosome BTA26. Five QTLs associated with FM305 contained two genes on chromosomes BTA8 (GNA14) and BTA15 (LRRC4C). In this study, no overlapping regions were detected between LPL, DO, and FM305 traits. Recently, QTL database ( http://www.animalgenome.org/cgi-bin/QTLdb/BT/index ) reported that QTL for LPL were found on every chromosome; QTL for DO were found on chromosome 2, 4, and 18; and QTL for FM305 were found on chromosomes 3, 4, 5, 6, 7, 12, 13, 14, 16, 18, 19, 23, and 29.

Table 3: Summary of the QTL regions and genes associated with the length of productive life (LPL), days open (DO), and 305-days milk yield of first lactation (FM305) in crossbred Holstein dairy cattle.

LPL is also referred to as the longevity of cattle. Meszáros et al. (2014) reported that the longevity of cattle is affected by the SYT10 gene, which is essential for the release of insulin-like growth factor IGF1. However, in this study, we found that the major QTL affecting LPL variation was the synaptotagmin I (SYT1) gene, which is located on BTA5. This gene is involved in calcium metabolism (Flori et al., 2009). According to Ubach et al. (2001), the C2B domain of SYT1 gene functions as a calcium sensor in neurotransmitter release. A previous study reported that this QTL region was related to feed intake in pigs (Do et al., 2013). The genes DOCK11 (dedicator of cytokinesis 11) and IL13RA1 (interleukin 13 receptor, alpha1) located on chromosome X region lied within in the gene cluster DOCK11-IL13RA1-AF074402-LOC616260 that could be targeted for fat and protein yields (Cole et al., 2011). The KLHL13 (kelch-like 13) gene belongs to the KLHL gene family that comprises 42 KLHL genes classified by the HUGO Gene Nomenclature Committee (HGNC). The KLHL genes are associated with cancer in humans. However, the function of KLHL13 is not clear yet (Dhanoa et al., 2013). Longevity is highly influenced by poor reproductive performance and high milk production. Therefore, the genes SYT1, DOCK11, and IL13RA1 might affect longevity in dairy cattle.

Our analysis identified only one gene associated with DO. The protein kinase, cGMP-dependent, type I (PRKG1) gene located on chromosome BTA26 is known to play a role in the key genes and causal mutations affecting milk fatty acid traits in dairy cattle (Li et al., 2014). In contrast, another study found that the FAM181A gene was associated with DO (Nayeri et al., 2016).

Two genes were associated with FM305. The GNA14 (guanine nucleotide binding protein, G protein) is the gene associated with calcium signaling pathway, as identified by the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway (Noyes et al., 2011). The LRRC4C (leucine-rich repeat containing 4C) gene has a homologous gene, netrin-G1 (NGL-1), in humans. The function of the NGL-1 receptor was to promote the outgrowth of thalamocortical axons (TCA), and the membrane-bound netrins can participate in the receipt of signals from axonal signaling pathways (Lin et al., 2003). Jensen et al. (2013) found that LRRC4C is the major gene relating to the different expressions between the controls, S. aureus infected, and E. coli infected quarters. We found no association between QTLs and milk production trait indicating that the DGAT1 gene on chromosome BTA14 affects production traits, as has been suggested in several reports (Sun et al., 2009; Minozzi et al., 2013; Nayeri et al., 2016).

In conclusion, the heritability estimates of LPL and DO were low in this study. Thus, genetic progress by means of single-trait selection alone would be slow. However, fertility traits, such as DO, had a very low genetic correlation (r_g) with LPL. In future, other fertility traits must be investigated to determine the appropriate selection index. Differences in breed and population might cause GWAS results to differ from those reported in other studies. The putative QTL regions identified in this study are novel and related to important functional traits in the dairy cattle industry. Hence, these QTLs should be considered for use in marker-assisted selection in Holstein dairy cattle.

Conflicts of interest

The authors declare that they have no conflict of interest.

Acknowledgments

Research supported by the Higher Education Research Promotion and National Research University Project of Thailand, Office of the Higher Education Commission, through the Food and Functional Food Research Cluster of Khon Kaen University and the Thermo- Tolerant Dairy Cattle Research Group at Khon Kaen University. The authors gratefully thank Dr. Chok Bullakul, Chokchai Dairy Farm, Nakhonratchasima, Thailand, for supporting the analysis data.