Genome-wide association study for selected cheese-making properties in Dual-Purpose Belgian Blue cows

This study aimed to estimate genetic parameters and identify genomic region(s) associated with selected cheese-making properties (CMP) in Dual-Purpose Belgian Blue (DPBB) cows. Edited data were 46,301 test-day records of milk yield, fat percentage, protein percentage, casein percentage, milk calcium content (CC), coagulation time (CT), curd firmness after 30 min from rennet addition (a30), and milk titratable acidity (MTA) collected from 2014 to 2020 on 4,077 first-parity (26,027 test-day records), and 3,258 second-parity DPBB cows (20,274 test-day records) distributed in 124 herds in the Walloon Region of Belgium. Data of 28,266 SNP, located on 29 Bos taurus autosomes (BTA) of 1,699 animals 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 25 consecutive SNPs (with an average size of ~2 Mb) was calculated, and regions accounting for at least 1.0% of the total additive genetic variance were used to search for candidate genes. Heritability estimates for the included CMP ranged from 0.19 (CC) to 0.50 (MTA), and 0.24 (CC) to 0.41 (MTA) in the first and second parity, respectively. The genetic correlation estimated between CT and a30 varied from −0.61 to −0.41 and from −0.55 to −0.38 in the first and second lactations, respectively. Negative genetic correlations were found between CT and milk yield and composition, while those estimated between curd firmness and milk composition were positive. Genome-wide association analyses results identified 4 genomic regions (BTA1, BTA3, BTA7, and BTA11) associated with the considered CMP. The identified genomic regions showed contrasting results between parities and among the different stages of each parity. It suggests that different sets of candidate genes underlie the phenotypic expression of the considered CMP between parities and lactation stages of each parity. The findings of this study can be used for future implementation and use of genomic evaluation to improve the cheese-making traits in DPBB cows.


INTRODUCTION
In most milk-producing countries, a large and growing fraction of the milk produced is used for making cheese (Ikonen et al., 1999;Wedholm et al., 2006;Cassandro et al., 2008).
In 1980, the world cheese production was 8.7 million tonnes, and this expanded to 11.4 million tonnes in 1990, 20 million tonnes in 2011, and 21 million tonnes in 2014. The world cheese production is expected to continue to increase and the expansion from 2014 to 2023 is expected to be 25%, equalizing 5.3 million tonnes of cheese (PM Food & Dairy Consulting, 2014). The composition and physicochemical properties of milk have marked effects on its processing properties and the quality of resulting dairy products. The quality, especially the milk coagulation properties (MCP), and the composition of milk, in particular the fat percentage (FP) and protein percentage (PP), used for cheese making have significant effects on the amount and the quality of produced cheese through their influences on rennet coagulation of milk and cheese texture (Guinee et al., 2006;Abd El-Gawad and Ahmed, 2011;Bittante et al., 2012;Pretto et al., 2013;Panthi et al., 2017). Commonly, the main MCP studied are milk rennet coagulation time (RCT, min), which is the time from the addition of coagulant to milk to the beginning of coagulation, the time to a curd firmness of 20 mm (k20), and curd firmness at 30 min after coagulant addition (a30, mm) (Ikonen et al., 2004;Cassandro et al., 2008;De Marchi et al., 2009;Bittante, 2011). Many factors including milk fat, protein, and minerals concentrations and MTA can influence all phases of milk coagulation as well as the cheese yield and texture Penasa et al., 2016;Panthi et al., 2017). In addition, the cheese-making traits can be affected by environmental factors such as feeding, udder health, season, physiological stage (e.g., parity, lactation stage), but they are also genetically influenced (Cassandro et al., 2008;Cecchinato et al., 2011;Tiezzi et al., 2013). Therefore, genetic selection can be used as a cost-effective and permanent strategy to improve the cheese-making traits, and consequently, to improve the quality and the amount of cheese yield per volume of milk.
Although Holstein is the most important cattle breed reared in the Walloon Region of Belgium, the Dual-Purpose Belgian Blue (DPBB), which is popular on organic and direct marketing farms in Wallonia, is known as the second most important cattle breed in the Walloon Region of Belgium. However, Atashi et al. (2022a) showed a continuous decrease in the effective population size of DPBB, indicating the need for effective conservation programs in this local breed. Compared with Holstein cows, the DPBB cows are considered less demanding, having better fertility and longevity, and higher milk quality, but are still less productive (Colinet, 2010;Mota et al., 2017;Bastin et al., 2020). Therefore, the improvement of economically important traits is a critical step to convince more farmers to use DPBB cattle in their farms, and consequently, to preserve the breed. In this sense, making targeted modifications to the cheese-making properties (CMP) can significantly contribute to the production of dairy products with improved added value; however, there is no information about the genetic variability of CMP in DPBB cows. Furthermore, identification of genomic regions and individual genes responsible for genetic variation in CMP will improve the understanding of biological pathways involved and can be used for improving cheese-making traits, and therefore, the quality and the amount of cheese yield per volume of milk. The aim of this study was to estimate genetic parameters and identify the genomic regions associated with cheese-making traits including milk casein percentage (CNP), milk calcium content (CC), milk coagulation time (CT), milk curd firmness after 30 min from rennet addition (a30), and MTA in DPBB cows.

MATERIALS AND METHODS
The data used in the present study were obtained from a preexisting database managed by the Walloon Breeders Association (elevéo; 5590 Ciney, Belgium); therefore, ethics approval was not required in advance of conducting the study.

Phenotypic Data
The data used consisted of test-day records of milk yield (MY), FP, PP, mid-infrared (MIR) prediction of CNP, CC, MTA, CT, and a30. These traits have already been reported by the Walloon Breeders Association to their farmers who are interested in monitoring the milk quality for making cheese and butter at the farm (Bastin et al., 2020).
The MIR prediction equations used in this study were obtained from various studies (Soyeurt et al., 2009;Colinet et al., 2010Colinet et al., , 2015. Table 1 shows the calibration and the cross-validation statistics of the used MIR prediction equations. The MIR prediction equations of the studied traits were applied on standardized spectra from individual milk samples collected in the frame of the Walloon milk recording scheme. Full cross-validation was performed to assess the accuracy of the developed equations. The equation for predicting milk CC using milk MIR data was developed by Soyeurt et al. (2009). Briefly, the CC of milk samples were measured by inductively coupled plasma atomic emission spectrometry. Then, using partial least squares, a calibration equation was built to estimate the contents of calcium in milk. The cross-validation coefficients of determination (R 2 ) and the ratio performance/deviation (RPD) of cross-validation defined, as the ratio of standard deviation (SD) to standard error of cross-validation, was 0.82 and 2.34, respectively. The partial least squares regression was used to develop the equation to predict MTA using milk MIR data . The R 2 and the RPD of cross-validation was 0.68 and 1.77, respectively. The R 2 and the RPD of cross-validation estimated for the equation used to predict milk CNP were 0.95 and 5.47, respectively. The equations used to predict milk coagulation parameters using MIR data were developed by Colinet et al. (2015). The R 2 and the RPD of cross-validation were, respectively, 0.63 and 1.64 for RCT, and 0.42 and 1.31 for a30.
The data were collected from 2014 to 2020 on 4,077 animals distributed in 124 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 and second parities that had data for all included traits on a given test-day were kept. Records from DIM lower than 5 d and greater than 365 d were eliminated. For cow, if parity 2 was present, parity 1 was also present. Age at the first calving was calculated as the difference between birth date and first calving date and restricted to the range of 640 to 1,500 d. Daily MY, FP, and PP were restricted to range from second-parity cows were 26,027 (on 4,077 lactations) and 20,274 (on 3,258 lactations), respectively.

Genotypic Data
Genotypic data were available for 1,699 animals (639 males and 1,060 females). Individuals were genotyped using the BovineSNP50 Beadchip v1 to v3 and EuroG MD (SI) v9 (Illumina). The SNP in common among the 4 chips were kept. Non-mapped SNP, SNP located on sexual chromosomes, and triallelic SNPs were excluded. The minimum GenCall Score of 0.15 and minimum GenTrain Score of 0.55 were used to keep SNP (Wilmot et al., 2022). The SNP with Mendelian conflicts, and minor allele frequency less than 5% were excluded. The difference between observed and expected heterozygosity was estimated, and if the difference was greater than 0.15, the SNP was excluded (Wiggans et al., 2009). In total, 28,266 SNP located on 29 BTA were used in the genomic analysis.

Variance Component Estimation
The (co)variance components and breeding value for CMP were estimated based on the integration of the random regression test-day model (RR-TDM) into the single-step genomic BLUP procedure (SS RR-TDM) using a model adapted for DPBB data structure (Atashi et al., 2021). The following SS RR-TDM through multiple-trait (4 traits), multiple-lactation (first 2 lactations) was used: where y ijklmn is the test-day record (MY, FP, PP, CNP, CC, MTA, CT, and a30) on DIM n of cow m in parity l, belonging to ith class of HTDp, jth class of AS, and kth class of HY; HTDp is the fixed effect of herd-testday-parity; AS is the fixed effect of age-season of calving defined as follows: age at calving class (6 and 4 classes of calving age were created for the first and second parity, respectively) × season of calving (4 seasons: from January to March, from April to June, from July to September, and from October to December), t are, respectively, the random regression coefficients of herd-year of calving, permanent environmental, and additive effects modeled using Legendre polynomials of order 2, and e ijklmn is the residual effect. The herd-year of calving, permanent environment, additive genetic, and residual variances were obtained as following: where HY is the 24 × 24 covariance matrix of the herd-year of calving regression coefficients; I is an identity matrix, ⊗ represents the Kronecker product function, P is the 24 × 24 covariance matrix of the permanent environmental regression coefficients; Ga is the 24 × 24 covariance matrix of the additive genetic regression coefficients, blocks within R R p = + ∑ contain diagonal matrices (4 × 4) of residual covariances among Atashi et  Coagulation time is defined here as the sum of the rennet coagulation time plus the time to a curd firmness of 20 mm (k20) measured by the computerized renneting meter.
traits with elements that depend on parity (p). Residual covariances among traits on the same test-day were therefore allowed to be different from zero, and residual covariances were the same within each parity. The H is a matrix that combines pedigree and genomic relationships, of which the inverse consists of 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 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 minor allele frequency 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 that 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 100,000 iterates. The first 20,000 iterates of the chain were regarded as a burn-in period to allow sampling from the proper marginal distributions. Convergence of the Gibbs chain and the length of burn-in was checked using postgibbsf90 (Misztal et al., 2002). Genetic (co)variances on each test-day were calculated using the equation described by Jamrozik and Schaeffer (1997). Daily heritability (h 2 ) was defined as the ratio of total additive genetic variance to the sum of additive genetic, permanent environmental, herd-year calving, and residual variances at a given DIM.

Prediction of Daily Genomic Breeding Values
The vector of genomic estimated breeding values (GEBV) of CMP for each animal i, which included daily GEBV from all DIM (5 to 365) in each parity, was estimated by multiplying the vector of predicted regression coefficients of the additive genetic by the matrix of Legendre orthogonal polynomial covariates; that is, where ĝ i is the vector of predicted regression coefficients of the additive genetic for animal i, and T is a matrix of orthogonal covariates associated with the Legendre orthogonal polynomial functions.

Genome-Wide Association Analyses
The GWAS analyses were performed for all included CMP in the first and second parities considering the following 3 lactation stages in each parity: (1) from 5 to 60 DIM, representing the ascending production stage and lactation peak; (2) from 61 to 200 DIM, representing the lactation persistency stage; and (3) from 201 to 365 DIM, representing the production decline until the end of the lactation (Oliveira et al., 2019). Therefore, GEBV for each lactation stage of each animal i (for each trait in each parity) was obtained by averaging the daily GEBV solutions of the specific DIM; that is, where GEBV î , 1 GEBV î , 2 and GEBV î 3 are the GEBV for the first, second, and third lactation stages of animal i obtained by averaging the GEBV from 5 to 60, 61 to 200, and 201 to 365 DIM, respectively. The SNP effects for each lactation stage were estimated individually for each trait in each parity 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), σ u 2 is the total 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 following equation: where λ was defined by VanRaden (2008) as a normalizing constant, as described below: The percentage of 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 25 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 results were presented as the proportion of variance explained by each window of 25 adjacent SNPs with an average size of ~2 Mb. To identify important genomic regions related to the analyzed traits, moving windows (SNPby-SNP) of 25 adjacent SNPs explaining at least 1.0% of the total additive genetic variance were selected for further analyses.

Candidate Genes and Functional Analyses
In total, 1,699 (639 males and 1,060 females) of the included animals were genotyped using the BovineSNP50 Beadchip v1 to v3 and EuroG MD (SI) v9 (Illumina), 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. The Lift Genome Annotations tool, available through a simple web interface (https: / / genome .ucsc .edu/ cgi -bin/ hgLiftOver) was used to convert genomic coordinates of the identified regions from the UMD3.1 to the ARS-UCD1.2 assembly. Then, to identify possible candidate genes associated with considered milk's cheese-making traits, 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.2 assembly) were further investigated. The candidate genes were identified using the National Center for Biotechnology Information Map Viewer tool for the ARS-UCD1.2 assembly as the reference map. Moreover, the list of the identified candidate genes was uploaded to Enrichr for gene ontology analysis (Chen et al., 2013;Kuleshov et al., 2016). Significantly enriched terms with at least 3 genes from the input gene list were identified based on the retrieved adjusted P value.
The lactation curves for the considered traits are presented in Figure 1. The peak of MY occurred at 4 wk in milk at both parities (15.6 and 19.6 kg in the first and second parity, respectively). The minimum values of FP, PP, CNP, CC, a30, and MTA were found at 6 wk in milk at both parities. The studied milk composition (FP, PP, CNP, and CC) reached the maximum level at the end part of the lactation.

Heritabilities and Genetic Correlations
The h 2 calculated for the considered traits in the first and second parity are presented in Tables 3 and  4, respectively. The mean (SD) estimates daily h 2 for MY, FP, and PP were, respectively, 0.37 (0.10), 0.29 (0.06), and 0.37 (0.02) in the first parity. The corresponding h 2 values estimated for the second parity were 0.36 (0.05), 0.28 (0.02), and 0.32 (0.01), respectively. The mean (SD) estimated daily h 2 for considered CMP varied between 0.19 (0.07) and 0.50 (0.07), and 0.24 (0.10) and 0.41 (0.04), respectively, in first and second parity. Among the considered CMP, MTA was the most heritable, and CC had the lowest h 2 .
Genetic correlations estimated among the considered traits in the first and second parity are presented in Tables 3 and 4, respectively. Low to medium negative genetic correlations were found between MY and milk composition (FP, PP, CNP, CC) in both parities. Mean genetic correlations between MY and MTA was low in both lactations. Low positive genetic correlations were found between CC and MTA and a30; however, the genetic correlations between CC and CT were weak and negative. The genetic correlation between CT and a30 was negative in the whole period of the lactation. Negative genetic correlations were found between CT and MY and milk composition, whereas genetic correlations between curd firmness and milk composition were positive.

Genome-Wide Association Study
General information about the results of single-step GWAS for the included CMP in the 3 defined lactation stages of the first and second parities are described in Data S1-S30 [5 traits (CC, CNP, CT, a30, and MTA) × 2 parities × 3 stages per parity]. The windows associated with the included CMP along with corresponding genes are presented in Table 5. In total, 4 genomic regions (BTA1, BTA3, BTA7, and BTA11) were identified to be associated with the studied CMP (Figures 2,  3, 4, 5, and 6). . This window explained 1.11 to 1.39% and 0.90 to 1.41% of the total additive genetic variance of CC in the first and second lactation, respectively.

BTA7
The genomic region located between 40.6 to 42.5 Mb on BTA7 was associated with a30 in both parities. This window was 1.88 Mb in size, contained 56 genes and  Coagulation time is defined here as the sum of the rennet coagulation time plus the time to a curd firmness of 20 mm (k20) measured by the computerized renneting meter. explained 1.15% of the total additive genetic variance of a30 in the second stage of the first parity. This region also explained more than 1.00% of the total additive genetic variance of a30 in the first and second stages of the second lactation.

BTA11
The genomic region located between 13.7 to 15.5 Mb on BTA11 was associated with CT in both lactations. This window was 1.78 Mb in size, contained 16 genes and explained 1.18 to 1.37 and 1.05 to 1.40% of the total additive genetic variance of CT in the first and second lactation, respectively. This window harbors the latent transforming growth factor β-binding protein 1 located in positions 15.47 to 15.93 Mb (LTBP1), xan-thine dehydrogenase located in positions 14.16 to 14.22 Mb (XDH), and tetratricopeptide repeat domain 27 located in positions 15.18 to 15.36 Mb (TTC27) genes. The gene ontology analyses showed the candidate genes identified for CMP enriched biological process terms related to the regulation of aspartic-type endopeptidase activity.

DISCUSSION
Milk's CMP are among important breeding traits in both dairy and dual-purpose cattle breeds, especially in modern animal husbandry environments. This study is the first genomic analysis performed on milk's cheesemaking traits in a dual-purpose breed of cattle. Although many genetic and genomic studies have been   Coagulation time is defined here as the sum of the rennet coagulation time plus the time to a curd firmness of 20 mm (k20) measured by the computerized renneting meter.
3 Curd firmness is defined as the curd firmness measured 30 min after enzyme addition by the computerized renneting meter. 4 The GWAS analyses were performed for all considered traits for each parity, considering 3 stages of lactation: (1) from 5 to 60 DIM, representing the ascending production stage and lactation peak; performed on CMP in dairy cattle, no study has focused on dual-purpose breeds; therefore, we compared our results with studies of dairy cattle breeds. The mean estimated daily h 2 for CMP varied between 0.19 to 0.50 in the first parity and between 0.24 to 0.41 in the second parity. We found that our estimates of h 2 were globally in agreement with previously reported results for international breeds of dairy cattle (Vallas et al., 2010;Tiezzi et al., 2013Tiezzi et al., , 2015. This study showed that mean daily h 2 estimates for CT, a30, and MTA were, respectively, 0.48, 0.36, and 0.50 in the first parity and 0.40, 0.39, and 0.41 in the second parity. The high heritability estimated for MCP supports possible genetic improvement of milk coagulation ability in this local breed. Tiezzi et al. (2013) reported that h 2 of RCT and a30 were, respectively, 0.21 and 0.24, in Italian  Holstein cows. Cecchinato et al. (2011) reported that h 2 of RCT, a30, and MTA ranged, respectively, from 0.22 to 0.58, 0.05 to 0.32, and 0.11 to 0.42 in Holstein, and from 0.12 to 0.32, 0.08 to 0.27, and 0.11 to 0.30 in Brown Swiss dairy cows. The variation found for h 2 of CMP in the literature can be explained by the differences in the studied breeds, structure of the data, number of records, statistical models, and the length of the period of data collection.
Curd firmness and coagulation time showed weak negative genetic correlations with MY; therefore, selection for improved curd firmness may be associated with a slightly reduced milk production. This is consistent with the results of previous studies in which weak genetic correlations among coagulation properties (CT and a30) and MY were reported for Finnish Ayrshire (Ikonen et al., 1999(Ikonen et al., , 2004, Italian Holstein (Cassandro et al., 2008), and Estonian Holstein (Vallas et al., 2010).   The composition of milk, in particular the concentration of fat and protein, have significant effects on the amount and the quality of produced cheese through their influences on rennet coagulation of milk (Guinee et al., 2006;Abd El-Gawad and Ahmed, 2011;Bittante et al., 2012;Pretto et al., 2013;Panthi et al., 2017;Troch et al., 2017). This study showed moderate to high positive genetic correlations between a30 and fat, protein, and casein contents. Some previous studies reported moderate to high positive genetic correlations between a30 and protein or casein contents (Cassandro et al., 2008;Tiezzi et al., 2015). However, Ikonen et al. (1999) reported negative genetic correlation between a30 and protein (casein) contents and Ikonen et al. (2004) reported no genetic correlation between milk protein (casein) content and a30. This study showed negative genetic correlations between CT and milk composition. Some studies reported that short RCT was genetically associated with high protein (casein) content (Cassandro et al., 2008;Tiezzi et al., 2015), however some other studies reported that long RCT is genetically correlated with high protein content (Ikonen et al., 1999;Vallas et al., 2010). Several factors such as sample size, the studied breeds, models and methods of estimations, and variation across laboratories might explain these inconsistencies.
The genetic correlation estimated between CC and a30 was weak and positive, whereas a low negative genetic correlation was between CC and CT. Milk calcium concentration is a factor affecting the cheese texture through influencing the RCT, the size of the casein aggregates, strength of the clot and body and texture of cheese (Santos et al., 2013).
The results showed that desirable MCP (i.e., short coagulation time and high curd firmness) were strongly associated with acidity of milk, measured as titratable acidity, which agrees with results from previous studies (Ikonen et al., 1999(Ikonen et al., , 2004Cassandro et al., 2008). The genetic correlation between CT and a30 was always negative which agrees with those reported by other studies (Cassandro et al., 2008;Tiezzi et al., 2015).
The GWAS results identified 4 genomic regions associated with the considered CMP. On BTA1, a window from 142.0 to 143.4 Mb was associated with MTA, CT, curd firmness, and CNP. Previous studies reported that this region was associated with MY, FY, PY, FP, PP in Holstein cows (Yang et al., 2015;Ibeagha-Awemu et al., 2016;Iung et al., 2019), milk lactalbumin percentage in Holstein, Normande, and Montbéliarde cows (Sanchez et al., 2017), and CNP, milk pH, and curd firmness in Montbéliarde cows (Sanchez et al., 2019). Atashi et al. (2022b) reported that this region was associated with PP in DPBB cow. This region harbors the SLC37A1 and PDE9A genes. The SLC37A1 gene encodes a glucose-6-phosphate transporter involved in the glucose blood homeostasis (Pan et al., 2011) and is highly expressed in the mammary gland (Kemper et al., 2015;Raven et al., 2016). Previous studies found that the SLC37A1was associated with MY, FP, and PP (Kemper et al., 2015;Raven et al., 2016;Pausch et al., 2017), and phosphorus concentration . Sanchez et al. (2019) reported that SLC37A1 was associated with milk's CMP, and milk mineral content in Montbéliarde cows (Sanchez et al., 2019). The SLC37A1 plays a role in the inorganic anion transport and is a good candidate for CMP and milk composition (Sanchez et al., 2019). The PDE9A gene is a cGMPspecific phosphodiesterase, which is highly expressed in mammary gland and was associated with MY, FY, and PY in Holstein cows (Yang et al., 2015).
The genomic region between 14.5 to 16.9 Mb on BTA3 was associated with CC. Previous studies showed that this region is associated with FY, PY, FP, PP, milk lactose percentage, milk fatty acids (FA) contents, and milk magnesium content in various cattle breeds (Kemper et al., 2018;Benedet et al., 2019;Jiang et al., 2019;Sanchez et al., 2019). This region harbors genes including ARHGEF2, DENND4B, GON4L, SYT11, ADAM15, GBA, and EFNA1. It has been reported that ARHGEF2, DENND4B, GON4L, SYT11, and ADAM15 are associated with PP (Jiang et al., 2019), GBA is associated with FP and PP (Raven et al., 2014;Jiang et al., 2019), and EFNA1 is associated with MY, FP, and PP in Holstein cows (Jiang et al., 2019). The results of a high-throughput analysis showed that AR-HGEF2 gene is associated with bovine mastitis (Bagnicka et al., 2021).
A genomic region located on BTA7 position 40.6 to 42.5 Mb was associated with curd firmness. The SNP located in this region has been showed to be associated with curd firmness, CC, CNP, and milk whey PP in Montbéliarde (Sanchez et al., 2019), MY, FP, fat yield, PP, protein yield, CNP, SCC, milk FA contents, and milk riboflavin content in Holstein cows (Schopen et al., 2009;Xu et al., 2014;Poulsen et al., 2015;Durán Aguilar et al., 2017). The NLR family pyrin domain containing 3 (NLRP3) gene is located in this region and was associated with the dairy capacity composite index and angularity in Holstein dairy cows (Kolbehdari et al., 2008).
The region located between 13.7 to 15.5 Mb on BTA11 was associated with CT. Previous studies showed that SNP located in this region is associated with MY (Meredith et al., 2012), FY, PY, FP, PP (Cole et al., 2011), milk FA contents (Pegolo et al., 2016), clinical mastitis (Schulman et al., 2009), and SCS (Wang et al., 2015) in dairy cows. This window harbors genes including LTBP1, XDH, and TTC27. The LTBP1 which regulates transforming growth factor-β activity (Saharinen et al., 1999) has been reported to be associated with milk FA content (Iung et al., 2019). It has been reported that TTC27 was associated with FY in Holstein (Cai et al., 2018), and XDH was associated with milk FA contents in Brown Swiss cows (Pegolo et al., 2016).
Caseins is the most abundant milk proteins and determine the quantity and quality of cheese (Wedholm et al., 2006). The 4 caseins (α S1 -CN, α S2 -CN, β-CN, and κ-CN) account for more than 75% of the whole bovine milk protein. They are encoded by 4 genes mapped to BTA6 in a tightly linked 250-kb cluster (Hayes and Petit, 1993;Popescu et al., 1996). The gene order is CSN1S1 (α S2 -CN-encoding gene), CSN2 (β-CNencoding gene), CSN1S2 (α S2 -CN-encoding gene), and CSN3 (κ-CN-encoding gene; Threadgill and Womack, 1990). Although, pervious studied found that milk protein content and composition are associated with SNP inside these genes (Schopen et al., 2009;Dadousis et al., 2016;Zhou et al., 2019;Mohan et al., 2021), we found no association between the region harboring CSN genes and the studied traits. The SNP file used in this study contained no SNP inside the CSN genes, which can explain why this region showed no association with the studied CMP.
The SNP located in DGAT1 gene (SNP rs109421300 located on position 1801116 bp of BTA14) was included in the used genomic data; however, windows harboring this SNP did not explain a considerable part of the total additive genetic variance of the studied traits. Although, the DGAT1 gene has been found to affect MY and composition (Ashwell et al., 2004;Maxa et al., 2012;Oliveira et al., 2019), the association between DGAT1 and production traits varied among breeds and populations. Factor such as breeds, genetic linkage phase, statistical methods, definition of the traits, as well as environmental conditions, such as season, nutrition, and management can contribute to the variability of GWAS results in the studied breeds and populations (Iung et al., 2019). Oliveira et al. (2019) reported that genomic regions associated with MY were located on BTA11, BTA16, and BTA28 for the Ayrshire; BTA6, BTA13, and BTA14 (many genes including DGAT1) for the Holstein; and BTA2, BTA5, BTA11, BTA20, BTA26, and BTA27 for the Jersey breed. Atashi et al. (2020) used data of 777,000 SNP on a large number of Holstein cows and found no association between windows harboring SNPs inside the DGAT1 gene and MY traits. Dadousis et al. (2016) showed no association between DGAT1 gene and MCP, curd firmness, PP, and acidity in Italian Brown Swiss cows.
In this study, we found traits with high genetic correlations showing different associated genomic regions. For example, CNP and CC had a genetic correlation of >0.50; however, CNP was associated with BTA1, whereas CC was linked to BTA3. Thus, perhaps, the different genomic regions identified on genetically correlated traits could be attributed to the noncorrelated part between the traits. However, MTA, CT, a30, and CNP were found to be genetically correlated and were all associated with a genomic region on BTA1.
Results of gene ontology analyses showed that biological process terms related to the regulation of aspartic-type endopeptidase activity were enriched by the candidate genes identified for CMP which is consistent with the role of aspartic proteinases in cheese-making process. Aspartic proteinases are present in all living organisms and paly very important roles including protein processing, maturation, and degradation, and are widely used as milk-coagulating agents in industrial cheese production (Yegin and Dekker, 2013).

CONCLUSIONS
This study is believed to be the first study to estimated genetic parameters and to identify genomic region associated with milk's cheese-making traits in a dualpurpose breed of cattle. The findings showed that the included CMP are moderately heritable and could be included into the breeding program currently used for DPBB breed. Genetic correlations of the CT with milk production and composition traits were negative and weak; and moderate to high positive genetic correlation was estimated between a30 and milk composition. Therefore, selection for shorter milk CT and higher curd firmness is expected to result in a higher milk composition. The MTA showed a moderate to high heritability, weak genetic correlations with milk production traits, strong genetic correlation with coagulation properties, and can be measured easily; therefore, enhancement of MCP could be achieved through indirect selection based on this indicator trait. Using available genotyping data, additional analyses were done to identified genomic regions associated with selected milk's cheesemaking traits. The identified regions showed contrasting results between parities and among different stages of each lactation; therefore, it might be worthwhile to analyze CMP per lactation stage as different traits. Relevant information supporting the results not presented here is provided in supplemental data (https: / / github .com/ hadiatashi/ Atashi -et -al -cheese -making -dpbb -2022). None of the data were deposited in an official repository because they are the property of the breeding organizations, and they may be available upon reasonable request. Authors are credited with the following contributions: Hadi Atashi, conceptualization, formal analysis, investigation, methodology, writing the original draft; Catherine Bastin, development of original trait definitions, reviewing and editing the draft; Hélène Wilmot, reviewing and editing the draft; Sylvie Vanderick, reviewing and editing the draft; Xavier Hubin, data curation, reviewing and editing the draft; Nicolas Gengler, data curation, project administration, reviewing and editing the draft. The authors have not stated any conflicts of interest.