Single-step genome-wide association analyses for milk urea concentration in Walloon Holstein

The aims of this study were to conduct a single-step genome-wide association to identify genomic regions associated with milk urea ( MU ) and to estimate genetic correlations between MU and milk yield ( MY ), milk composition (calcium content ( CC ), fat percentage ( FP ), protein percentage ( PP ), and casein percentage ( CNP )), and cheese-making properties ( CMP ) (coagulation time ( CT ), curd firmness after 30 min from rennet addition ( a30 ), and titratable acidity ( TA )). The used data have been collected from 2015 to 2020 on 78,073 first-parity (485,218 test-day records), and 48,766 s-parity (284,942 test-day records) Holstein cows distributed in 671 herds in the Walloon Region of Belgium. Data of 565,533 single nucleotide polymorphisms ( SNP ), located on 29 Bos taurus autosomes ( BTA ) of 6,617 animals (1,712 males) were used. Random regression test-day models were used to estimate genetic parameters through the Bayesian Gibbs sampling method. The SNP solutions were estimated using a single-step genomic BLUP approach. The proportion of the total additive genetic variance explained by windows of 50 consecutive SNPs (with an average size of ~216 Kb) was calculated, and the top-3 genomic regions explaining the largest rate of the genetic variance were considered promising regions and used to identify potential candidate genes. Mean (standard deviation) MU was 25.38 (8.02) mg/dl and 25.03 (8,06) mg/dl in the first and second lactation, respectively. Mean heritability estimates for daily MU were 0.21 and 0.23 for the first and second lactation, respectively. The genetic correlations estimated between MU and MY, milk composition, and CMP were quite low (ranged from −0.10 (CC) to 0.10 (TA) and −0.05 (CT) to 0.13 (TA) for the first and second lactations, respectively). The top-3 regions associated with MU were located from 80.61 to 80.74 Mb on BTA6, 103.26 to 103.41 Mb on BTA11, and 1.59 to 2.15 Mb on BTA14. Genes including KCNT1 , MROH1 , SHARPIN , TSSK5 , CPSF1 , HSF1 , TONSL , DGAT1 , SCX , and MAF1 were identified as positional candidate genes for MU. The findings of this study can be used for a better understanding of the genomic architecture underlying MU in Holstein cattle.


INTRODUCTION
Nitrogen (N) emissions from livestock can have harmful effects on biodiversity, soil, water and air quality, as they account for 23% of global nitrous oxide and 60% of global ammonia emissions (Uwizeye et al., 2020).Nitrogen is released to the environment through the feces and urine of livestock as undigested N fractions and as the end products of N metabolism, respectively (Müller et al., 2021).Nitrogen, as oxides or ammonia, is one of the greenhouses gases contributing to air pollution and through leaching to rivers and ground water resources.Therefore, reducing N excretion by dairy cattle is desirable due to global concerns about the agricultural contribution to environmental pollution by N. The first requirement for conducting a successful genetic selection is to establish a method to measure the traits or phenotypes of interest on large number of animals with a low cost.
The breakdown of protein, both in the rumen and in the small intestine, results in the production of ammonia, which is then detoxified by conversion to urea in the liver or kidneys.The urea can then be recycled into the rumen via saliva or by direct diffusion across the rumen wall, be excreted in urine (urinary urea; UU), or be secreted in milk (milk urea; MU) (Alio et al., 2000).Although most urea is eliminated from blood with urine excretion, UU is difficult to measure.MU can be measured at a low cost and is included as a standard part in most milk recording systems and is positively correlated with UU (Burgos et al., 2007, Burgos et al., 2010).Furthermore, the content of urea in milk represents approximately 48% of NPN in milk and is known as an important parameter of milk quality.Moreover, MU is associated with milk yield composition (e.g., fat and protein contents) (Chen et al., 2022a).It has been reported that MU is associated with cheese-making properties (CMP) and that higher MU is associated with greater coagulation time required for milk (Guinot-Thomas, 1992).In addition, MU has been shown to be correlated with nitrogen utilization efficiency and is considered as an on-farm indicator to guide nutritional strategies in dairy cows (Schepers and Meijer, 1998, Nousiainen et al., 2004, Chen et al., 2022a).There is a strong positive relationship between blood urea (BU) nitrogen and MU (Butler et al., 1996, van den Berg et al., 2021).Furthermore, it has been shown that an increased BU is associated with a decreased uterine pH and a decreased in the probability of conception in dairy cows (Elrod and Butler, 1993, Elrod et al., 1993, Butler et al., 1996).It has been reported that an increased MU is associated with a decreased reproductive performance in dairy cattle (Rajala-Schultz et al., 2001, Guo et al., 2004, Hojman et al., 2004).
Reducing MU can have many benefits, such as reducing N emissions, improving feed efficiency and animal health, milk quality and reproductive performance; therefore, this trait could be a potential selection criterion for breeding dairy cattle (Ariyarathne et al., 2021).However, a comprehensive knowledge of the genetic background of MU is needed before it is included into the breeding program.Furthermore, identification of genomic regions and individual genes responsible for genetic variation in MU will improve our understanding about the biological pathways involved and can be used for decreasing MU.
Although MU has been included in the regular milk recording of dairy cows in the Walloon Region of Belgium for more than 2 decades, the potential contribution of this trait in the breeding program of Walloon dairy cows is still under investigation.The genetic parameters for MU and its genetic correlations with traditional traits have previously been investigated in Walloon Holstein (Chen et al., 2021(Chen et al., , 2022b)); however, no comprehensive study has been conducted to investigate the genomic architecture of MU for the whole period of lactation and the potential genetic relationships between MU and milk technological traits (e.g., CMP) has not been investigated in this population of dairy cattle.Therefore, the aims of this study were to identify genomic regions associated with MU concentration and to estimate genetic correlation between MU and selected cheese-making traits in Walloon Holstein cows.

Phenotypic Data
The data used consisted of test-day records of traits including milk yield (MY), milk urea (MU), midinfrared (MIR) prediction of fat percentage (FP), protein percentage (PP), casein percentage (CNP), calcium content (CC), titratable acidity (TA), coagulation time (CT), and curd firmness at 30 min after coagulant addition (a30).MIR prediction equations used to predict CNP, CC, TA, CT, and a30 were obtained from various studies (Soyeurt et al., 2009, Colinet et al., 2010, Colinet et al., 2013, Colinet et al., 2015).The data were collected from 2015 to 2020 on 78,073 animals distributed in 671 herds in the Walloon Region of Belgium.Data were edited to include only cows with known birth date, calving date, and parity number.Only records from the first 2 parities that had data for all included traits on a given test-day were kept.Records from days in milk (DIM) lower than 5 d and greater than 365 d were eliminated.Age at the first calving (AFC) was calculated as the difference between birth date and first calving date and restricted to the range of 540 to 1200 d.Daily MY, FP, and PP were restricted to range from 3 to 99 kg, 1.5 to 9% and 1 to 7%, respectively (ICAR, 2022).Test-day records of the rest of the considered traits were edited to remove records outside the range of mean ± 3 standard deviations (SD).The number of test-day records in the first-and second-parity cows were 485,218 (on 78,073 cows) and 284,942 (on 48,766 cows), respectively.On average across the data set, 5.51 test-day records were available per cow per lactation.Pedigree depth of the animals was traced back to 5 generations.Full pedigree records included 186,548 females and 10,076 males.

Genotypic Data
Genotypic data were available for 6,617 (1,712 males) phenotyped animals or those animals included in the pedigree.Individuals were genotyped using the BovineSNP50 Beadchip v1 to v3 and EuroG MD (SI) v9 (Illumina, San Diego, CA, USA).Single nucleotide polymorphisms (SNP) in common among the 4 chips were kept.Non-mapped SNP, SNP located on sexual chromosomes, and triallelic SNPs were excluded.A minimum GenCall Score of 0.15 and a minimum Gen-Train Score of 0.55 were used to keep SNP.Then, the genotypes were imputed to HD with a reference population of 4,352 HD individuals (1,046 males) using FImpute V2.2 software (Sargolzaei et al., 2014).SNP with Mendelian conflicts, and those with minor allele frequency (MAF) less than 5% were excluded.The Atashi et al.: Genomic analysis of milk urea in Holstein difference between the observed and expected heterozygosity was estimated, and if the difference was greater than 0.15, the SNP was excluded (Wiggans et al., 2009).Finally, 565,533 SNPs located on 29 Bos taurus autosomes (BTA) were used in the genomic analyses.

Variance Component Estimation
The (co)variance components and breeding values for MU were estimated based on the integration of the random regression test-day model (RR-TDM) into the single-step GBLUP procedure (SS RR-TDM) using the following single-trait, multiple-lactation (first 2 lactations) model (Paiva et al., 2022):   Φ are, respectively, the random regression coefficients of herd-year of calving, permanent environmental, and additive effects modeled using Legendre polynomials of order 3; and e ijklmn is the residual effect.The herd-year of calving, permanent environment, additive genetic, and residual variances were assumed as follows: Var where HY is the 8 × 8 covariance matrix of the herdyear of calving regression coefficients; I is an identity matrix; ⊗ represents the Kronecker product function; P is the 8 × 8 covariance matrix of the permanent environmental regression coefficients; Ga is the 8 × 8 covariance matrix of the additive genetic regression coefficients; blocks within R r p = + ∑ contain residual variance (r) that depends on parity (p).Residual variance was assumed to be the same within each parity.The H is a matrix that combines pedigree and genomic relationships, which its inverse consists on the integration of additive and genomic relationship matrices, A and G, respectively (Aguilar et al., 2010): where A is the numerator relationship matrix based on the pedigree for all animals; A 22 is the numerator relationship matrix for genotyped animals; and G is the weighted genomic relationship matrix obtained using the following function: The G* is the genomic relationship matrix obtained using the following function described by VanRaden (2008): where Z is a matrix of gene content adjusted for allele frequencies (0, 1 or 2 for aa, Aa and AA, respectively); D is a diagonal matrix of weights for SNP variances (D = I); M is the number of SNPs, and p i is the MAF of the ith SNP.The H matrix was built scaling G based on A 22 considering that the average of the diagonal of G is equal to the average of the diagonal of A 22 and, the average of the off-diagonal G is equal to the average of the off-diagonal A 22 .The (co)variance components were estimated by Bayesian inference using the GIBBS3F90 software (Aguilar et al., 2018).Gibbs sampling was used to obtain marginal posterior distributions for the various parameters using a single chain of 200,000 iterates with a sampling interval of 20 samples.The first 50,000 iterates of the chain were regarded as a burn-in period to allow sampling from the proper marginal distributions.Genetic (co)variances on each test-day were calculated using the equation described by Jamrozik and Schaeffer (1997).Daily heritability was defined as the ratio of genetic variance to the sum of the additive genetic, permanent environmental, herd-year calving, and residual variances at a given DIM.
The vector of genomic estimated breeding values (GEBV) of MU for each animal i, which included daily GEBV from all DIM (5 to 365) in each parity, was estimated by multiplying the vector of additive genetic predicted regression coefficients by the matrix of Legendre orthogonal polynomial covariates; that is, , where ĝi is the vector of additive genetic predicted regression coefficients for animal i; and T is a matrix of orthogonal covariates associated with the Legendre orthogonal polynomial functions.
In addition, the same RR-TDM through multipletrait (2 traits), multiple-lactation (first 2 lactations) was used to estimate correlation between MU and the studied traits of milk yield, milk composition, and CMP.Moreover, phenotypic trend for MU was obtained by linear regression of phenotypic values over calving year.

Genome-Wide Association Study
GEBV of animal i through the entire lactation (DIM 5 to 365) were obtained by averaging the daily GEBV solutions of all DIM; that is, where GEBV i ˆ is the GEBV of animal i through the entire lactation, obtained by averaging the GEBV from 5 to 365.
The SNP effects were estimated using the postGSf90 software (Aguilar et al., 2014).The animal effects were decomposed into those for genotyped (a g ) and ungenotyped animals (a n ).The animal effects of genotyped animals are a function of the SNP effects, a g = Zu, where Z is a matrix relating genotypes of each locus and u is a vector of the SNP marker effect.The variance of animal effects was assumed as: where D is a diagonal matrix of weights for variances of markers (D = I) and σ u 2 is the additive genetic variance captured by each SNP marker when the weighted relationship matrix (G) was built with no weight.The SNP effects were obtained using the following equation:

Z
where λ was defined by VanRaden ( 2008) as a normalizing constant, as described below: The percentage of the total additive genetic variance explained by the ith genomic region was estimated as following: where a i is the genetic value of the ith region that consists of 50 adjacent SNPs; σ a 2 is the total additive genetic variance; Z j is the vector of the SNP content of the jth SNP for all individuals; and ûj is the marker effect of the jth SNP within the ith region.The proportion of the total additive genetic variance explained by windows of 50 consecutive SNPs (with an average size of ~216 Kb) was calculated, and the top-3 genomic regions explaining the largest rate of the total additive genetic variance were considered promising regions and used to identify potential candidate genes (Soares et al., 2021).The concept of grouping SNP into windows was adopted as a way to better capture the genetic information such as the extent of linkage disequilibrium (LD) in neighboring SNPs (Habier et al., 2011).In addition, we performed a GWAS with frequentist p-values to detect single SNP associated with MU.The thresholds of the Bonferroni corrected p values for 5% and 1% genome-wide significance associations were set as 8.84E −8 (0.05 divided by the number of SNPs) and 1.74E −8 (0.01 divided by the number of SNPs), respectively.

Identification of positional candidate genes for the studied traits
The animals included in this study were genotyped using the BovineSNP50 Beadchip v1 to v3 and EuroG MD (SI) v9 (Illumina, San Diego, CA, USA); then, the genotypes were imputed to the BovineHD Beadchip which are based on the bovine reference genomes assembly UMD3.1.However, new bovine reference genome assembly ARS-UCD1.2,assembled using long sequencing reads, filled gaps and resolved repetitive regions of the UMD3.1 assembly, and has more credible annotation information (Rosen et al., 2020).The Lift Genome Annotations tool, available through a simple web interface (https: / / genome .ucsc.edu/cgi Atashi et al.: Genomic analysis of milk urea in Holstein -bin/ hgLiftOver) was used to convert coordinate ranges of the identified genomic regions from the UMD3.1 to the ARS-UCD1.2assembly.Then, to identify possible positional candidate genes associated with MU, genes located within the identified genomic regions (i.e., between the start and end of genomic coordinates of the identified regions based on the ARS-UCD1.2assembly) were further investigated.We identified genes using the National Center for Biotechnology Information (NCBI) Map Viewer tool for the ARS-UCD1.2assembly as the reference map.The list of identified positional candidate genes was uploaded to Enrichr to perform gene ontology (GO) enrichment analyses (Kuleshov et al., 2016).Significantly enriched terms were identified based on the retrieved adjusted P value.

RESULTS
The descriptive statistics for studied traits are presented in Table 1.Mean (SD) MU was 25.38 (8.02) and 25.03 (8.06) mg/dl for the first and second lactation, respectively.Daily MY averaged 24.2 kg (4.01%fat, 3.38% protein, and 2.59% casein), and 27.9 kg (4.10% fat, 3.46% protein, and 2.63% casein) in the first and second lactations.The trend in concentrations of MU across DIM is given in Figure 1.The lowest value for MU was found at the beginning of lactation, increased rapidly by increasing DIM, reached its peak at the middle of lactation, then slightly decreased by increasing DIM to the end of the lactation.The coefficient of variation (CV%) of MU was 32% in both lactations.Environmental factors including herd, calving year, calving season, age at first calving, parity and lactation stage affected MU (P < 0.05).Phenotypic trend, obtained by linear regression of MU phenotypic values over calving year, showed a positive trend of 0.73 and 0.79 mg/dl milk for the first and second lactations, respectively (Figure S1) (https: / / github .com/hadiatashi/ Holstein -MU).
Mean heritability estimates for daily MU was 0.21 (SD = 0.02; range = 0.15 to 0.23), and 0.23 (SD = 0.02; range = 0.19 to 0.24) for the first and second lactation, respectively.Genetic correlations between MU in the first and second lactation ranged from 0.84 to 0.97 with a mean of 0.95.Genetic correlations estimated between MU and the studied traits are presented in Table 2. Genetic correlations between MU and milk yield traits ranged from −0.10 (CC) to 0.05 (MY) and −0.02 (CC) to 0.12 (FP) for the first and second lactation, respectively.Low genetic correlations ranged from 0.00 (a30) to 0.10 (TA) and from −0.05 (CT) to 0.13 (TA) were estimated between MU and cheese-making traits in the first and second lactation, respectively.
General information about the results of single-step GWAS for MU concentration are described in Data S1 and Data S2 (https: / / github .com/hadiatashi/ Holstein -MU).The Manhattan plots of the proportion of total additive genetic variance explained by 50-SNP windows are shown in Figures 2. The top 3 windows explaining the highest percentage of the genetic variance for MU along with corresponding genes are presented in Table 3.The identified regions were located from 80.61 to 80.74 Mb on BTA6, 103.26 to 103.41 Mb on BTA11, and 1.59 to 2.15 Mb on BTA14 (UMD3.1 assembly) and combined explained 0.99% and 1.04% of the total additive genetic variance for MU at the first and second lactations, respectively.The region identified on BTA6 was 0.13 Mb in size and explained 0.34% of the total additive genetic   method were the same as those identified using the window-based GWAS procedure (Figure S2) (https: / / github .com/hadiatashi/ Holstein -MU).

DISCUSSION
Mean MU in the first and second lactations was 25.38 and 25.03 mg/dl, respectively, which is in close agreement with those previously reported for Holstein (Bastin et al., 2009, Rzewuska and Strabel, 2015, Atashi and Hostens, 2021, Chen et al., 2021) and the Italian Brown Swiss cows (Samoré et al., 2007), but higher than that reported for the Dual-Purpose Belgian Blue (DPBB) cows (Atashi et al., 2022).The coefficient of variation was 32%.The CV for MU was 26% in Belgian Holstein (Atashi and Hostens, 2021), 33% in Dutch Holstein (Stoop et al., 2007), and 44% to 45% in the DPBB cows (Atashi et al., 2022).The trend in concentration of MU across DIM showed that lowest value for MU was found at the beginning of lactation, then increased rapidly by increasing DIM, reached the peak at the middle of lactation, then slightly decreased by increasing DIM to the end of the lactation.It is well documented that dairy cows are in a negative energy balance state in the early stage of the lactation when the feed intake cannot meet their requirements.During this time, it can be assumed that cows use N as efficiently as possible which can explain, at least in part, the obtained results (Chen et al., 2021).
Similar lactation curve pattern has been reported for MU by previous studies (Mucha and Strandberg, 2011, Rzewuska and Strabel, 2013, Atashi and Hostens, 2021, Atashi et al., 2022).Mean heritability estimates for daily MU were 0.21 and 0.23 for the first and second lactation, respectively.The heritability reported in the literature for MU range from 0.09 to 0.47, depending on the method used to measure the MU content, the population considered, the parity, and the model used to estimate variance components.
High genetic correlation was found between MU in first and second lactation which is in a close agreement previous studies (Miglior et al., 2007, Atashi andHostens, 2021).The results showed that MU is weakly correlated with milk yield and composition.Genetic correlation between MU and milk yield traits has been investigated by several researchers, however the results are inconsistent.Atashi and Hostens (2021) reported that genetic correlations between MU with milk yield traits were close to zero in Holstein dairy cows.Godden et al. (2001) reported a positive genetic correlation between MY and MU.Diab and Hillers (1996) reported a negative relationship between MU and milk yield.Samoré et al. (2007) reported that MU had a positive genetic relationship with FY (0.12), null with PY (0.03) and a negative genetic correlation with MY (−0.17) in Italian Brown Swiss dairy cattle.Research on the relationship between MU and cheese-making traits is limited.This study showed low genetic correlations between MU and cheese-making traits.Mean genetic correlation estimates between MU and TA ranged from 0.10 to 0.13.Martin et al. (1997) reported that MU influences directly acidification kinetics, chemical composition and texture characteristics of cheeses.Guinot-Thomas (1992) reported that higher MU is associated with greater coagulation time required for milk.The variation found for genetic correlation of MU with productive traits in the literature can be explained by the differences in structure of the data, number of records, the statistical models used, and the length of the period of data collection.
Typically, GWAS methods are based on testing the significance of SNP effects on the traits of interest.However, SNPs within a genomic region can be highly correlated and jointly influence the phenotype.Furthermore, the genetic information in neighboring SNPs, such as the extent of LD, is not considered in single SNP-based association study (Bao and Wang, 2017).Therefore, window-based GWAS procedure have been proposed as an effective method to estimate the combined effect of several consecutive SNPs in a specific region and to identify genomic regions that explain a given amount of genetic variance (Aguilar et al., 2019).The common form for declaring significance is to use a threshold on the additive genetic variance explained by individual window (Aguilar et al., 2019).However, it is unclear what window size is optimal, and no standard presently exists to define the threshold on explained genetic variance.Therefore, determining the proper window size is usually subjective and researchers often have not justified their choices or sometimes have acknowledged that their choices are arbitrary (Beissinger et al., 2015).Medeiros de Oliveira Silva et al. ( 2017), using the BovineHD SNP panel, considered 50-adjacent SNP windows (with an average of 280 kb) that explained at least 0.50% of additive genetic variance as the threshold to declare significance.Atashi et al. (2020), using the BovineHD SNP panel, considered 50-adjacent SNP windows that explained more than 1% of the total additive genetic variance as the threshold to declare significance for milk production and lactation curve parameters.Han and Peñagaricano (2016) considered 1.5-Mb windows that explained at least 0.50% of the total genetic variance as the threshold to declare significance.Suwannasing et al. (2018) considered windows that explained more than 1% of the total genetic variance as the threshold to declare significance.Tiezzi et al. (2015) calculated the variance absorbed by 10-SNP moving windows and reported the Atashi et al.: Genomic analysis of milk urea in Holstein 10 windows explaining the largest amount of genomic variance as most important windows.In this study, a window-based GWAS through the single-step genomic best linear unbiased predictor (ssGBLUP) was used.The results were presented by the proportion of total genetic variance explained by window of 50 adjacent SNP with an average size of ~216 Kb and he top-3 genomic regions explaining the largest rate of the total additive genetic variance were considered promising regions and used to identify potential candidate genes.We used 1 SNP as the moving step of the window, which ensured that we do not miss genomic regions potentially associated with the traits due to the combination of SNPs.The top-3 identified genomic regions were located from 80.61 to 80.74 Mb on BTA6, 103.26 to 103.41 Mb on BTA11, and 1.59 to 2.15 Mb on BTA14 and combined explained 0.99% and 1.04% of the total additive genetic variance for MU at the first and second lactations, respectively.The region identified on BTA6 was 0.13 Mb in size and explained 0.34% of the total additive genetic variance of MU in the first and second lactation.It has been reported that SNPs inside this region are associated with traits including MY, fat yield (FY), protein yield (PY), milk fatty acid (FA) profile and CNP (Cole et al., 2011, Pausch et al., 2017, Sanchez et al., 2017b, Benedet et al., 2019).van den Berg et al. (2022) reported that this region is associated with MU nitrogen in Australian and New Zealand dairy cattle.Ma et al. (2023) reported that SNPs inside this region are associated with MU nitrogen in Holstein cows.
The genomic region located from 103.26 to 103.41 Mb on BTA11 was associated with MU in the first and second lactations.Previous studies showed that SNPs inside this region are associated with MY, FY, PY, FP, PP, CNP, somatic cell score (SCS), and CMP traits (Bennewitz et al., 2004, Chen et al., 2006, Kučerová et al., 2006, Schopen et al., 2009).This region harbors KCNT1 gene which encodes a sodium-activated potassium channel subunit which is thought to function in ion conductance and developmental signaling pathways.KCNT1 gene has been reported to be associated with MY, FY, and PY in US Holstein (Cole et al., 2011).
The genomic region located from 1.59 to 2.15 Mb on BTA14 was associated with MU in the first and second lactation.Previous studies showed that SNPs inside this region are associated with MY, FY, PY, FP, PP, SCS, milk FA profile, and CNP in breeds of dairy cattle (Conte et al., 2010, Buitenhuis et al., 2016, Pedrosa et al., 2021).Chen et al. (2022b) reported that this region is associated with nitrogen efficiency index and milk true protein nitrogen in Holstein cows.Clancey et al. (2019) reported that SNP inside this region are associated with MY in Holstein cows.This region was 0.56 Mb in size and harbors genes including the scratch family zinc finger 1 (SCRT1), diacylglycerol O-acyltransferase 1 (DGAT1), cleavage and polyadenylation specific factor 1 (CPSF1), tonsoku like, DNA repair protein (TONSL), and spermatogenesis and centriole associated 1 (SPATC1).The DGAT1, involves in the last step of the synthesis of triacylglycerol, has a major effect on milk production traits (Jiang et al., 2010, Maxa et al., 2012, Clancey et al., 2019, Cruz et al., 2019).Nayeri et al. (2016) reported that SNP located within CPSF1 and TONSL are associated with MY in Canadian Holstein cattle.The SHANK associated RH domain interactor gene (SHARPIN) is mapped inside this region.This product of this gene involved in the regulation of immune and inflammatory responses (Wang et al., 2012) and has been reported to be associated with the colostrum and serum albumin concentrations in Holstein cows (Lin et al., 2020).Sanchez et al. (2017a) reported DGAT1, maestro heat like repeat family member 1 (MROH1), and ribosomal biogenesis factor (BOP1) as the most important genes explaining the majority of the variability of milk protein composition in Montbéliarde, Normande and Holstein dairy cattle.

CONCLUSION
This study showed that MU is moderately heritable and is weakly correlated with MY, milk composition, and cheese-making traits; therfore, selection for lower MU concentration would not influnce the production performance.In general, windows explained less than 0.50% of the total additive genetic variance of MU and these low contributing regions were spread across the genome.This indicates that MU is highly polygenic, in which many regions across the genome contribute to its genetic variation.80,611,741,948 78,982,989,732 BTA11 103,264,409,247 103,309,336 The positions of the identified genomic regions based on the UMD3.1 assembly. 3 The positions of the identified genomic regions based on the ARS-UCD1.2assembly.

Figure 2 .
Figure 2. Additive genetic variance explained by windows of 50 adjacent SNP across chromosomes for milk urea concentration for the first (first row) and second (second row) lactation in Walloon Holstein cows.The red lines showed that the most important 3 windows explained more than 0.25% of the total additive genetic variance of MU in both lactations.
Atashi et al.: Genomic analysis of milk urea in Holstein

Table 1 .
Atashi et al.:Genomic analysis of milk urea in Holstein Descriptive statistics for milk yield traits and cheese making properties in Walloon Holstein cows 1 2Milk titratable acidity in Dornic degree (°D).

Table 2 .
Genetic correlation estimated between MU and milk yield, milk composition and cheese-making traits in Walloon Holstein cows1Coagulation time is defined here as the sum of the rennet coagulation time (RCT) plus the time to a curd firmness of 20 mm (k20) measured by the Computerized Renneting Meter.

Table 3 .
Genomic regions associated to milk urea concentration in the first and second lactations for Walloon Holstein cows 1