Quantitative Trait Loci Affecting Milk Production Traits on BTA 14 in Iranian Holstein Dairy Cattle: A Confirmation

INTRODUCTION

The most important economic traits such as milk production and disease resistance have a continuous underlying distribution of measurement and the genes responsible for genetic variation in these traits have been termed quantitative trait loci (QTL) (Georges et al. 1995). A major objective of the QTL studies in livestock breeding schemes is to detect the markers that can be implemented in breeding programs via marker-assisted selection (MAS). Milk production is one of the economically important sectors in dairy cattle industries. In Iran, milk consumption per capita is about half the world average (85-90 liter). Iranian dairy production has undergone significant and considerable structural changes during the last two decades due to creating larger herds. In 1966, total cattle population of Iran was 4.9 million but it has reached 7.9 million cattle (mainly Holstein) in the last three decades (Ebrahimi-Hoseinzadeh et al. 2015a). In 2017 year, cattle population of Iran 1 million mainly and 5/4 million blended reported (Association of animal feed and aquatic animals, http://t.me/irfia_channel). The QTL mapping has mainly focused on milk production traits. Simulation studies have shown that MAS has a significant effect on the genetic progress of these traits (Kashi et al. 1990; Meuwissen and Goddard, 1996). Many other important QTLs have been reported by linkage information since the first QTL experiment on milk production was reported by Georges et al. (1995). Even though QTLs for milk production traits have been found in almost every bovine autosome chromosomes, but BTA14 is one of the most concerned chromosomes in QTL mapping linked with the milk production traits such as milk yield, fat percentage, fat content, fat yield, protein percentage and protein yield (Arranz et al. 1998; Ashwell et al. 1998; Khatkar et al. 2004). In early 2003-2018, via genome scan based on marker-QTL linkage analysis and reported 9239 QTLs were located on chromosome 14 (http://www.animalgenome.org/cgi-bin/QTLdb/index). It has also been shown that BTA14 region was found to harbor of diacylglycerol acyltransferase 1 (DGAT1) and cytochrome P450, family 11, subfamily B member 1 (CYP11B1) genes (Kaupe et al. 2007). DGAT1 is a microsomal enzyme that catalyzes the synthesis of triglycerides from diglycerides and acyl-coenzyme A (Cases et al. 1998; Furbass et al. 2006). Considering the role of DGAT1 in lipid metabolism and also QTL studies for milk production traits, DGAT1 has been regarded as a good positional candidate gene affecting milk production traits in dairy cattle (Grisart et al. 2002; Furbass et al. 2006). CYP11B1 encodes an enzyme called 11-beta-hydroxylase that is found in the adrenal glands and is positively associated with fat content (Kawamoto et al. 1990; Curnow et al. 1991). In Iranian Holstein dairy cattle compared to other countries, applying different breeding goals and selection criteria over generations may have resulted in differences in the genetic background of cattle compared to other populations. Although many studies have been conducted on Iranian cattle QTLs affecting milk yield, fat yield and fat percentage of these cattle have not been studied yet using random regression models (Alinaghizadeh et al. 2007; Javanmard et al. 2008; Mohammadi et al. 2009; Mohammad-Abadi et al. 2010; Kharrati-Koopaei et al. 2012a; Kharrati-Koopaei et al. 2012b; Pasandideh et al. 2015; Barazandeh et al. 2016). Hence, the main purpose of this study was to confirm the suggestive location and maps of the possible QTLs affecting milk yield, fat yield and fat percentage on chromosome 14 in Iranian Holstein dairy cattle using linkage (LA) and linkage disequilibrium (LDLA) information.

MATERIALS AND METHODS

Population and phenotypic data

In the current study, all the cows were in similar intensive dairy farms. The daughters were randomly selected from 10 sire family groups of Iranian Holstein dairy cattle across 58 herds located in the different regions of Iran (Tehran, Isfahan, Markazi, Khorasan, Zanjan, Fars, Yazd and Qazvin). Family sizes varied between 9 to 41 dairy cows with 24 daughters on the average. Semen samples of all sires were obtained from National Animal Breeding Center of Iran and blood samples were collected from 233 cows based on a daughter design approach (Table 1). Phenotypic data were collected for milk yield (MY), fat yield (FY) and fat percentage (FP). All phenotypic data were collected based on twice daily milking by National Animal Breeding Center and Promotion of Animal Products of Iran in 2009-2010.

DNA extraction and genotyping

Genomic DNA was purified from semen sires by using salting out method (Weyrich, 2012) and blood samples by using AccuPrep® Genomic DNA Extraction Kit (Bioneer, Korea), respectively. For genotyping 10 sires and their daughters, on chromosome 14 with a distance 0 to 63 cM, ten microsatellite markers (ILSTS011 (blue,6-FAM^TM), DIK2598 (blue,6-FAM^TM), DIK4884 (yellow,NED^TM), DIK5080 (green,VIC), CBDIKM002 (blue,6-FAM^TM), ILSTS039 (blue,6-FAM^TM), BM1508 (green,VIC), DIK4361 (yellow,NED^TM), CBDIKM004 (yellow,NED^TM) and CSSM066 (green,VIC) respectively label color and fluorescent lable) were used (Table 2). The average distance between markers in the linkage groups was 6.5 cM. PCR amplifications were performed in 15 μL reactions each containing 1.5 ng genomic DNA, 2.5 mM MgCl₂ (Fermentas, USA), 0.2 mM dNTPs mix (Roche Applied Science, Germany), 10 pmol each primer (Metabion, Germany), 1X PCR buffer (Fermentas, USA) and 0.5 U Taq DNA polymerase (Fermentas, USA). Amplifications were performed using the Gene Amp PCR 9700 (Applied Biosystems, USA). The PCR conditions were as follows: an initial denaturation for 5 min at 95 ˚C, followed by 30 cycles of 30 sec at 95 ˚C, 45 sec at 59 ˚C or 60 ˚C, 30 sec at 72 ˚C and a final extension of 10 min at 72 ˚C. Genotyping was performed using Applied Biosystems 3130 Genetic Analyzer. Allele sizes were determined using Gene Mapper v.4.0 (Applied Biosystems, USA).

Statistical analysis

Milk production records were analyzed using a mixed model in which QTL was considered as a random effect. The variance components of QTL effects were estimated to examine the putative QTL positions.

Table 1 Distribution of the informative samples of 10 half-sib sires families

Table 2 Details of the ten microsatellite markers used in current study

Gametic relationship matrix for QTL alleles of any two founder haplotypes was computed based on two methods: linkage analysis (LA) within sire groups (Wang et al. 1995) and joint analysis across sire groups using linkage (LA) and linkage disequilibrium (LDLA) information (Meuwissen and Goddard, 2001). The following equation was used as a generalized linear model:

y= Xb + Za + Wg + e

Where:

y: vector of observations (records of the traits).

b: vector of fixed effects (herd-year-season).

a: vector of random animal effects.

g: vector of random QTL effects.

e: vector of random residual effects.

X, Z, W: design matrix which related records to fixed, random animal and random QTL effects, respectively.

It was assumed that the genetic effects (a), QTL effects (g) and residual effects (e) were normally distributed as a ~ N (0, A σ²_a), g ~ N (0, HP σ²_g) and e ~ N (0, I σ²_e), respectively:

Where:

A: additive relationship matrix.

HP: identical by descent (IBD) matrix that contains the IBD probabilities for the QTL at location p.

I: an identity matrix.

The IBD probabilities for the QTLs based on the haplotypes were calculated using the analytical method suggested by Meuwissen and Goddard (2001). The DMU package (Madsen and Jensen, 2002) was used to estimate the (co)variance components based on the average information of restricted maximum likelihood. This algorithm approached the restricted likelihood to maximize with respect to variance components associated with the random effects (Sorensen et al. 2003). The parameters were estimated at the marker brackets along the chromosome region. Hypothesis tests to detect QTLs were performed with the likelihood ratio test (LRT) as follows:

LRT= -2 × [log likelihood_{(reduced) –} log likelihood_(full)]

Where:

log likelihood_(reduced): log likelihood of a model excluding the QTL effects (no QTL model).

log likelihood _(full): log likelihood of a model with a QTL effects.

It was assumed that LRT was a statistic with chi-squared distribution and degree of freedom (DF)= 1. According to 10 markers used in this study and with respect to possibility to calculate probability of inheritance of favorable allele of QTL in each location between two markers, 9 locations were assumed to be located between two markers. Then, LRT was computed using the DMU (Madsen and Jensen, 2002) and the results were compared with chi-square test and finally, the most likely location of QTLs was identified (P<0.05).

RESULTS AND DISCUSSION

Different family sizes in Table 1 indicated that factors such as the number of samples with genotypes and imbalance of allele segregation in markers could be more effective and informative (Esmailizadeh et al. 2006). The information content of the individuals is the proportion of animals in which the allele inherited from the sire can be unambiguously identified (Table 2). Three traits were evaluated to refine and confirm the detected QTL affecting milk yield and composition. This study successfully confirmed the mapped QTL for MY, FY and FP in Iranian Holstein dairy cattle. The genomic relationship information for mapping of the reported QTLs was based on both LA and LDLA approaches. QTL mapping using LD information across half-sib families and using LDLA can decrease the required number of phenotyped and genotyped individuals, compared to LA method (Ansari-Mahyari et al. 2009). Also, compared with LA, LDLA utilizes all historical genomic information in gametic relationship matrix between founders using linkage disequilibrium information and would be more informative for the inference (Hu et al. 2010). This method has been used widely in dairy cattle genomic studies. The significant QTLs in the most probable location (P<0.05) are shown in Table 3. Although the results from LA method indicated several QTLs for fat percentage (FP) could be located at 3 to 20 and 36 to 50 cM, the most likely location was at 3 cM (P<0.05). In addition, for fat yield (FY) a location was significant at 20 to 60 cM, while the most likely position was at 54 cM (P<0.05). For milk yield (MY), significant location was at 54 to 60 cM but the most likely position was revealed to be at 54 cM (P<0.05) (Figures 1 and 3). Using microsatellite markers in dairy cattle, many studies have found that segregated QTLs for milk production traits in several autosome bovine chromosomes (Khatkar et al. 2004). For example in different Holstein dairy cattle, several studies have reported the segregation of at least one QTL in the middle of BTA6, affecting one or more of the milk production traits, close to marker BM143 (Kühn et al. 1999; Velmala et al. 1999; Nadesalingam et al. 2001; Ron et al. 2001). To date, various QTLs affecting milk production traits have been identified on BTA14 in Holstein dairy cows. A major QTL affecting milk yield and composition on the centromeric end of BTA14 (DGAT1 K232A mutation) was first reported by Coppieters et al. (1998). The QTL mapping studies on BTA14 has showed that the markers were linked to the genes controlling important production traits. Other studies suggested thatILSTS039, CSSM066 and BM1508 markers were nearest markers to the fat percentage QTL (Heyen et al. 1999; Biochard et al. 2003; Ashwell et al. 2004). Most of the identified QTL regions for milk production traits in the current study have been previously reported in the other literature. Recently, Aliloo et al. (2015) confirmed genome regions with validated additive effects on milk yield on Bos taurus autosome 14 in Holstein and Jersey cows. Iso-Touru et al. (2016) reported that Bos taurus BTA14 is significantly (P<0.05) associated with fat yield in the DGAT1 gene, (region: from 1448510 bp to 2271832 bp) and milk yield (region: from 1448510 bp to 2271832 bp) in Nordic Red cattle using imputed whole genome sequence variants. They claimed that the strongest association with milk and fat yield was found on BTA14. Also, the chromosome 14 polymorphisms have been scanned for milk production traits in commercial half-sib families of French dairy goats. Martin et al. (2017) found 5 polymorphisms (positions: 0.033, 0.037, 0.111, 0.123 and 0.155 MB) that are significantly related to fat content and fat yield, base Illumina Goat SNP 50 Bead Chip. The results from the current study confirmed precisely those QTLs obtained by others in Iranian Holstein dairy cows. Besides, the QTL obtained by LDLA method showed that for milk yield a detection was observed at 12 to 20 and 60 cM, but significant positions were at 12 and 60 cM (P<0.05). Although, the most likely position affecting fat yield was detected at 60 cM on chromosome 14 (P<0.05) (Figures 2 and 3). Therefore, it could be concluded that current results on QTL mapping for milk and fat yield traits using LDLA information confirmed the result reported by Bennewitz et al. (2004) (marker:DGAT1). In addition, the results for milk yield QTL locations were in accordance with Bagnato et al. (2008) (marker:BMS1748) and Schnabel et al. (2005) (marker:BMS1899). Estimation of QTL, polygenic and residual variance components in 9 significant locations between each of the two markers (Table 4) indicated that QTL and polygenic variances for the most likely locations were greater than their residual variances. When applying LDLA, all variance components (genetic and non-genetic) were close to the true values and the estimates indicated lower standard error than LA (Ansari-Mahyari et al. 2009).

Table 3 Most probable locations (P<0.05) for MY, FY and FP on chromosome 14 using linkage analysis (LA) and linkage disequilibrium (LDLA) methods

Figure 1 Significant and likely positions of the QTL for MY, FY, and FP using LA method (MY: milk yield; FP: fat percentage; FY: fat yield; LRT: likelihood ratio test and LA: linkage analysis)

Figure 2 Significant and likely positions of the QTL for MY, FY, and FP using LDLA method (MY: milk yield; FP: fat percentage; FY: fat yield; LRT: likelihood ratio test and LDLA: linkage disequilibrium)

Figure 3 Likely positions estimated for milk yield and milk compositions traits mapped chromosome 14using LA and LDLA methods (MY: milk yield; FP: fat percentage; FY: fat yield; LA: linkage analysis and LDLA: linkage disequilibrium)

Table 4 Estimated variance components of the most likely locations and percentages of phenotypic variance of MY, FY and FP on chromosome 14 using linkage analysis (LA) and linkage disequilibrium (LDLA) methods

Based on LA method, the estimated QTL variance were 3877209 to 43575.5, 47.92 to 38.68, and 3.34 to 2.29 for MY, FY and FP, respectively. Based on the LDLA information, the estimated QTL FP respectively. Estimation of polygenic variance were 81331.1 to 88921.3, 108.007 to 88.65 and 6.8 to 5.8 for MY, FY, and FP, respectively. The polygenic variance for MY, FY and FP were 81331.1 to 88921.3, 108.007 to 88.65 and 6.8 to 5.8 respectively based on LA method, and 38044 to 41889.4, 52.67 to 61.07 and 1.92 to 1.77 respectively based on LDLA information. Residual variance size for MY, FY and FP ranged from 40901.3 to 39840.3, 56.86 to 48.41, and 2.58 to 2.50, respectively based on LA and 31754 to 32175, 45.41 to 46.63 and 2.55 to 2.52, for LDLA method. In Table 4, estimated variance components of the most likely locations for MY, FY and FP in two methods are shown. The results indicated that a greater proportion of the total variance can be explained by the effects, as also reported by Szyda et al. (2005). However, MY and FP had the largest trait percentages (ratio of QTL variance to phenotypic variance; Table 4), which can improve the power of QTL detection in QTL fine-mapping. The power of this analysis ranged from 0.15 to 0.35. The proportional QTL variance showed the contribution of specific trait loci to the total phenotypic variance of the trait. In each position mapping QTL using the desirable model determines whether a significant amount of the variance in a quantitative trait can be associated with the QTL in proposed position. QTL variance determined as the reduction of the residual variance acquired by fitting the QTL in the location was relatively small for the detection (Sohrabi et al. 2012). In the current study, according to variance component estimated in LA method, the largest ratio of polygenic and residual variances to phenotypic variance, QTL variance to residual variance and QTL variance to polygenic variance was observed for FP position at 60 cM, 50 cM and 60 cM, respectively. Also, in LDLA method, the largest ratio of polygenic variance to residual variance, QTL variance to residual variance and QTL variance to residual variance estimated for FP was refined at 54 cM, 54 cM and 60 cM, respectively.

CONCLUSION

In conclusion, the two methods considered in this study (LA and LDLA) showed a highly significant QTL located within 20 to 60 and 54 to 60 cM interval and within 12 to 60 and 60 cM interval for milk and fat yield traits, respectively. Also, a highly significant QTL for fat percentage was detected at 3, 12, 20, 36, 44 and 50 cM with LA method. The results suggest a possible use of these candidate markers in gene-assisted selection programs for improvement of milk production traits in Iranian dairy cattle, although large-scale studies in different breeds are required.

ACKNOWLEDGEMENT

This study was supported by genomics lab of Agricultural Biotechnology Research Institute of Iran (ABRII) and North Region Branch (ABRINI). Thus, their efforts are sincerely appreciated.