Genomic investigation of milk production in Italian buffalo

Abstract

The aim of this study was to test the feasibility of genomic selection in the Italian Mediterranean water buffalo, which is farmed mainly in the south Italy for milk, and mozzarella, production. A total of 498 animals were genotyped at 49,164 loci. Test day records (80,417) of milk (MY), fat (FY) and protein (PY) yields from 4127 cows, born between 1975 and 2009, were analysed in a three-trait model. Cows born in 2008 and 2009 with phenotypes and genotypes were selected as validation animals (n = 50). Variance components (VC) were estimated with BLUP and ssGBLUP. Heritabilities for BLUP were 0.25 ± 0.02 (MY), 0.16 ± 0.01 (FY) and 0.25 ± 0.01 (PY). Breeding values were computed using BLUP and ssGBLUP, using VC estimated from BLUP. ssGBLUP was applied in five scenarios, each with a different number of genotypes available: (A) bulls (35); (B) validation cows (50); (C) bulls and validation cows (85); (D) all genotyped cows (463); (E) all genotypes (498). Model validation was performed using the LR method: correlation, accuracy, dispersion, and bias statistics were calculated. Average correlations were 0.71 ± 0.02 and 0.82 ± 0.01 for BLUP and ssGBLUP-E, respectively. Accuracies were also higher in ssGBLUP-E (0.75 ± 0.03) compared to BLUP (0.57 ± 0.03). The best dispersions (i.e. closer to 1) were found for ssGBLUP-C. The use of genotypes only for the 35 bulls did not change the validation values compared to BLUP. Results of the present study, even if based on small number of animals, showed that the inclusion of genotypes of females can improve breeding values accuracy in the Italian Buffalo.

Highlights

• The genotypes of males did not improve the predictions.

• Genotypes of females improve breeding values accuracy.

• Slight increase in prediction accuracy with weighted ssGBLUP.

Keywords:

• Validation
• single-step
• LR method

Introduction

The domestic water buffalo is the second dairy species of the world, accounting for about 15% of the total milk production in 2018 (FAOSTAT 2020). Largest buffalo milk producers are India, Pakistan, and China. The world buffalo stock exhibited a relevant increase in the last six decades, from 88 million in 1961 to 206 million in 2018 (FAOSTAT 2020), respectively. Reasons for such an increase can be found in the high market value of buffalo milk and of its dairy products. Moreover, it should be remembered the great relevance of this species in many rural areas of the world as a source of meat and of traction power.

Two types of water buffalo exist, the river (Bubalus bubalis bubalis) and the swamp (Bubalus bubalis carabanensis) buffalo. The two types have different number of chromosomes, 50 for the river and 48 for the swamp, respectively (di Berardino and Iannuzzi 1981), geographical distribution, and they are characterised by phenotypic and genetic differences. River buffaloes are farmed mainly in the west, from India to Europe, whereas the swamp type can be found mostly in eastern Asian Countries (Iamartino et al. 2017). The river type consists mostly of distinct breeds as a result of selection whereas swamp consists of local unselected populations that have evolved by adaptation to a specific environment (Zhang et al. 2016).

River buffalo breeds have been selected for milk production traits, although organised breeding plans are currently running only in few countries (Zhang et al. 2020). Genetic selection in buffalo is hampered by the lack of reliable pedigree information (Rosati and Van Vleck 2002; Abdel-Shafty et al. 2020; Ghoreishifar et al. 2020), the difficult implementation of performance recording due to the farming structure that in several countries is based on small holders (Mokhber et al. 2018; Abdel-Shafty et al. 2020), and to the poor reproductive performances. As a consequence, the potential of buffalo has not been completely exploited (Mokhber et al. 2018).

The availability of a medium density (90 K) SNP platform for the buffalo (Iamartino et al. 2017) has opened the possibility to implement Genomic Selection (GS) programs also for this species. In particular the possibility of calculating the true realised relationships, even when pedigree information is poor or totally lacking (Hayes et al. 2009), could help in addressing a main issue of buffalo breeding. Moreover, the use of the single step GS approach (Aguilar et al. 2010) that combines phenotypic, pedigree and genomic information in the genomic BLUP framework could further help to improve breeding values accuracy in buffalo also when only a limited portion of the population is genotyped.

The genetic improvement of buffalo is a topic of great interest for the Italian livestock industry. Italy is the sixth Country in the world for Buffalo milk production and the sixteenth for number of farmed animals. The milk is mainly processed into mozzarella Campana cheese, that is recognised by Protected Designation of Origin (PDO) certification, and that has a high market value. Together with Brazil, Italy has experienced the largest relative increase (>200%) of the buffalo stock in the last sixty years. Buffaloes farmed in Italy are of the Mediterranean river breed, a group that includes also the breed from Mozambique and Romania (Colli et al. 2018). The Italian buffalo stock consists of 403,093 animals (www.vetinfo.it), and the dairy performance recording plan carried out by the Italian National Association of Buffalo Breeders (ANASB) involved 34,490 buffaloes in 2019. In the last 15 years, the average milk yield per lactation showed an improvement of about 8%, from kg 2184 in 2014 to kg 2356 in 2019, whereas fat (8.06 vs. 8.01) and protein (4.68 vs. 4.63) percentages remained quite constant (www.anasb.it).

In this work, the usefulness of genomic models for the estimation of variance components and the prediction of breeding values for milk production traits in the Italian Buffalo breed is evaluated. In particular, the traditional pedigree based BLUP was compared with the single-step genomic BLUP.

Materials and methods

Data used in this study were obtained from pre-existing databases and therefore the animal care approval was not needed.

Data

The whole data set (Table 1) consisted of 80,417 test day (TD) records of milk (MY, kg), fat (FY, kg) and protein (PY, kg) yields from 4127 buffalo cows born from 1975 to 2009 and farmed in 8 commercial herds.

Table 1. Basic statistics and phenotypic correlations.

	MY (kg)	FY (kg)	PY (kg)
Statistics
Mean	8.74	0.76	0.41
SD	3.27	0.27	11
Max	29.8	3.35	1.65
Phenotypic correlation
MY		0.83	0.97
FY			0.82
PY

MY: Milk yield; FY: Fat yield; PY: Protein yield.

Phenotypes were provided by ANASB. A sample of 463 buffalo cows of the whole data set and 35 bulls were genotyped with the Axiom Buffalo Genotyping Array 90 K (Iamartino et al. 2017). Quality control was applied both on animals and SNPs with the following criteria: call rate of animals and markers ≥95% (27,167 SNPs discarded); minor allele frequency ≥5% (8322 SNPs discarded), and P-value for the Hardy Weinberg equilibrium test <10⁻⁶ (167 SNPs discarded). Moreover, markers mapped on sexual chromosomes or without known position were removed. After quality control, all animals and 49,164 markers were retained for the analysis.

A pedigree with 7730 animals was provided by ANASB from which 3 generations were traced back from animals with records or genotypes, resulting in a total of 5404 animals which were included in the analyses.

Analysis

Data were analysed with the following 3-trait repeatability animal model: $${\mathrm{y}}_{\mathrm{t}}={\mathrm{Xb}}_{\mathrm{t}}+\mathrm{ }{\mathrm{W}}_1{\mathrm{g}}_{\mathrm{t}}+{\mathrm{W}}_2{\mathrm{u}}_{\mathrm{t}}+{\mathrm{W}}_3{\mathrm{p}}_{\mathrm{t}}+{\mathrm{e}}_{\mathrm{t}}$$ y_t is the vector of TD records for trait t (MY, FY, or PY) b_t is the vector of fixed effects of calving year (23 levels), parity (5 levels), birth year (27 levels), milking events per day (2 levels), calving month (12 levels) and days in milk (linear covariate); g_t, u_t and p_t are the vectors of herd-test-day, additive genetic and permanent environmental random effects, respectively; and e_t is the vector of random residuals. The X, W₁, W₂ and W₃ are the incidence matrices associating phenotypic records to the effects in b_t, g_t, u_t and p_t, respectively. (Co)variance structures of random effects were Iσ²_HTD, RELσ²_a, and Iσ²_pe, where I is an identity matrix, REL is the relationship matrix, σ²_HTD, σ²_a, and σ²_pe, are variance components associated to random effects of herd-test-day, additive genetic, and permanent environment, respectively.

According to the relationship matrix used, two models were used: (i) the pedigree-based (BLUP) with the numerator relationship matrix (A); (ii) the single step genomic BLUP (ssGBLUP) where A and the genomic relationship matrix (G) are blended into H, which inverse (${\mathrm{H}}^{-1}$) was built according to Aguilar et al. (2010): $${\mathrm{H}}^{-1}={\mathrm{A}}^{-1}+\left[\begin{array}{cc}0& 0\\ 0& {\mathrm{G}}^{-1}-{\mathrm{A}}_{22}^{-1}\end{array}\right]$$ where ${\mathrm{A}}^{-1}$ and ${\mathrm{G}}^{-1}$ are the inverses of the pedigree and genomic relationship matrices, respectively; ${\mathrm{A}}_{22}^{-1}$ is the inverse of the pedigree relationship matrix for genotyped animals only. The G matrix was created according to VanRaden (2008), with the following equation: $$\mathrm{G}\mathrm{ }=\frac{\mathrm{Z}\mathrm{D}{\mathrm{Z}}^{\mathrm{\prime }}}{k}$$ where Z is the matrix of centred gene contents, D is a diagonal matrix of SNP weights, equal to I in ssGBLUP, and k is the scaling parameter defined as $2\sum {\mathrm{p}}_{\mathrm{i}}\mathrm{(}1\mathrm{–}{\mathrm{p}}_{\mathrm{i}}\mathrm{)}$ where p_i is the allele frequency of the i-th SNP. To avoid singularity, G was blended with 5% of A₂₂ (VanRaden, 2008). Therefore, in the ssGBLUP model, SNP information is used to construct the G matrix among all the genotyped individuals.

Variance components were estimated with both methods using the aiREML algorithm implemented in the airemlf90 software (Misztal et al. 2014).

Breeding values were estimated with BLUP and ssGBLUP models. The genomic model was tested in five scenarios that differ in the number of genotyped animals:

A = genotypes available only for 35 bulls;

B = genotypes available only for the 50 candidates;

C = genotypes available for 50 candidates + 35 bulls;

D = genotypes available for 463 cows;

E = genotypes available for 463 cows + 35 bulls.

In all scenarios, candidate animals were 50 females born in 2008 and 2009 with phenotypes and genotypes (except scenario A) available.

To evaluate prediction accuracy and bias, the LR method was applied (Legarra and Reverter, 2018). Breeding values estimation was run twice: in the first run, candidate animals had their phenotypes available (whole data), whereas in the second they had their phenotypes masked (partial data). The following four statistics were computed to evaluate the models:

Dispersion: $${\mathrm{b}}_{\mathrm{w},\mathrm{p}}\mathrm{ }=\mathrm{ }\frac{\mathrm{cov}({\widehat{\mathrm{u}}}_{\mathrm{w}},{\widehat{\mathrm{u}}}_{\mathrm{p}})}{\mathrm{var}\left({\widehat{\mathrm{u}}}_{\mathrm{p}}\right)}$$ where ${\widehat{\mathrm{u}}}_{\mathrm{w}}$ and ${\widehat{\mathrm{u}}}_{\mathrm{p}}$ are the vectors of breeding values estimated for the candidates in the whole and partial datasets, respectively.

Prediction accuracy: $$\widehat{\mathrm{acc}}\mathrm{ }=\mathrm{ }\sqrt{\frac{\mathrm{cov}({\widehat{\mathrm{u}}}_{\mathrm{w}},{\widehat{\mathrm{u}}}_{\mathrm{p}})}{(1-\overline{\mathrm{F}}){\widehat{\mathrm{\sigma }}}_{\mathrm{u}}^2}}$$ where $\overline{\mathrm{F}}$ is the average inbreeding of the candidates and ${\widehat{\mathrm{\sigma }}}_{\mathrm{u}}^2$ is the additive genetic variance.

Correlation: $${\mathrm{\rho }}_{\mathrm{w},\mathrm{p}}=\text{cor}({\widehat{u}}_{\mathrm{w}},{\widehat{u}}_{\mathrm{p}})$$

Bias: $${\mu }_{wp}\mathrm{ }=\mathrm{ }{\stackrel{-}{\widehat{\mathrm{u}}}}_{\mathrm{p}}-\mathrm{ }{\stackrel{-}{\widehat{\mathrm{u}}}}_{\mathrm{w}}$$

Variance components estimated with the BLUP model were used in the validation. Convergence criterion was set to 10⁻¹⁵.

A further development is the weighted single-step (WssGBLUP) where different weights are assigned to SNP in the calculation of the genomic relationship matrix (Wang et al. 2012). In this procedure, the optimum weights are determined iteratively, with the first iteration being equivalent to ssGBLUP. In the first iteration, p-values for SNPs effects were computed according to Aguilar et al. (2019); Bonferroni correction, i.e. 0.05 divided by the number of SNPs (49,164), was set as significance threshold. Usually, 3–5 iterations are needed for obtaining the maximum prediction accuracy. Fragomeni et al. (2019) suggested that non-linear weights would be a better option in WssGBLUP compared to linear weights proposed by Wang et al. (2012). In this study, the non-linear A weights were calculated as follows (VanRaden 2008; Legarra et al. 2018): $${\mathrm{d}}_{\mathrm{i}}={\mathrm{CT}}^{\frac{\left|{\widehat{\mathrm{a}}}_{\mathrm{i}}\right|}{\mathrm{sd}\left(\widehat{\mathrm{a}}\right)}-2}$$ where d_i is the weight for SNP i; CT is a constant the determines the departure from normality (usually assumed for the SNP effects); $\left|{\widehat{\mathrm{a}}}_{\mathrm{i}}\right|$ and $\mathrm{sd}(\widehat{\mathrm{a}})$ are the absolute effect for SNP i and the standard deviation of the SNP effects, respectively. A CT value of 1.0 is equivalent to the unweighted ssGBLUP and deviations from it determine the departure from normality. Because this value is empirically determined, we tested three different values (1.025, 1.125, and 1.25) and the same LR validation was applied over 5 iterations of WssGBLUP to evaluate which weight and iteration combination would maximise accuracy of predictions with the least bias.

Results and discussion

Genetic parameters

Compared to values in Table 1, slightly lower average daily milk yield (7.33 ± 3.03 kg) was reported for the Nili-Ravi buffalo breed (Abdel-Shafy et al. 2020).

Heritability and repeatability estimates were very similar in both considered models (Table 2); lower values of genetic parameters were obtained for FY. As far as the kind of model is concerned, the ssGBLUP provided smaller estimates (Table 2) and smaller standard errors for genetic correlations. These differences were basically due to the larger additive genetic variance and the lower permanent environmental variance estimated with BLUP. On the contrary, residual variances were the same across methods, and variance associated to HTD were very similar (Table 2).

Table 2. Genetic parameters estimated using two different methods, BLUP and single-step GBLUP.

	BLUP			ssGBLUP
	MY	FY	PY	MY	FY	PY
Genetic parameters^a
HTD	0.24 ± 0.02	0.41 ± 0.02	0.06 ± 0.02	0.23 ± 0.01	0.40 ± 0.02	0.06 ± 0.01
Additive	1.87 ± 0.13	1.01 ± 0.08	0.38 ± 0.03	1.71 ± 0.11	0.91 ± 0.06	0.34 ± 0.02
Permanent	1.05 ± 0.10	0.73 ± 0.06	0.23 ± 0.02	1.11 ± 0.08	0.77 ± 0.04	0.24 ± 0.01
Residual	4.21 ± 0.02	4.08 ± 0.02	0.86 ± 0.02	4.21 ± 0.02	4.08 ± 0.02	0.86 ± 0.01
Heritability	0.25 ± 0.02	0.16 ± 0.01	0.25 ± 0.01	0.23 ± 0.01	0.15 ± 0.01	0.23 ± 0.01
Repeatability	0.40 ± 0.01	0.28 ± 0.01	0.40 ± 0.01	0.39 ± 0.01	0.27 ± 0.01	0.39 ± 0.01
Genetic correlation
MY		0.92 ± 0.02	0.98 ± 0.04		0.92 ± 0.01	0.97 ± 0.01
FY			0.96 ± 0.04			0.95 ± 0.01
PY

^aValues for FY and PY were multiplied by 10.

BLUP: best linear unbiased prediction; ssGBLUP: single-step genomic best linear unbiased prediction; MY: Milk yield; FY: Fat yield; PY: Protein yield; HTD: herd-test-day.

Genetic parameters of milk production traits have been estimated in buffalo with traditional pedigree-based BLUP methods, both on cumulated yields and on TD records. Lower h² values for lactational traits (0.14 for MY, 0.11 for FY, and 0.14 for PY) were reported for Italian Buffalo, whereas genetic correlations were similar to those of the present study (Rosati and Van Vleck 2002). Heritability values for lactational milk traits similar to those found in the present study with BLUP and sGBLUP models were estimated in Brazilian Murrah and Philippine buffaloes (Tonhati et al. 2008; Aspilcueta-Borquis et al. 2010; Flores et al. 2013). Lower genetic correlation values among yield traits were reported for Murrah buffaloes (Aspilcueta-Borquis et al. 2010). The modelling of TD records for estimating variance components resulted in lower heritability estimates. For example, Flores et al. (2013) reported values of 0.15, 0.09, and 0.03 for MY, FY and PY, respectively. Genetic correlations reported by these authors were 0.99 between MY and FY and 0.97 between MY and PY, respectively.

Genomic information has already been used to estimate variance components for milk traits in buffaloes. Liu et al. (2018) reported h² of 0.33, 0.35 and 0.27 for lactational MY, FY, and PY, respectively. Abdel-Shafy et al. (2020) estimated a slightly lower h² (0.20) for daily milk yield in Egyptian buffalo using a ssGBLUP approach. Some comparisons between BLUP and GBLUP were performed. Aspilcueta-Borquis et al. (2015) found similar heritabilities using the two methods, whereas Iamartino et al. (2017), reported h² for milk yield of 0.38 using the pedigree relationship matrix and 0.45 using a genomic relationship matrix.

In the present study, genetic parameters estimated with BLUP and ssGBLUP were quite similar. Unbiased variance components can be estimated when the whole data about the studied population is used in the computation. The unbiased estimations with complete data can be realised also in population under selection if the data used to make the selection decision are included in the analysis (Cantet et al. 2000). However, sometimes the use of complete dataset is not possible, either because the data is too big which makes variance components computationally expensive or because of incomplete and or incorrect pedigrees in the historical dataset. The ssGBLUP seems to be better than BLUP with subsets of data in simulated populations under selection (Cesarani et al. 2019). Recently, Aldridge et al. (2020) in a pig study (with 10 different traits, including mean litter birth weight, total number born, number still born and litter mortality) found very similar variance components using BLUP and different ssGBLUP models. However, because genomic estimates exhibited lower standard errors, they concluded that H is a more informative than A. Moreover, they also suggested that the genomic approach is to be preferred for more complex models, as those with repeated records where a permanent environment effect is fitted. Similar results for variance components estimated with BLUP and ssGBLUP were recently reported for two functional traits in Italian Simmental (Cesarani, Gaspa, et al. 2020). Hidalgo et al. (2020) investigated the changes in genetic parameters for fitness and growth traits in pigs under GS. These authors reported that the genetic parameters estimated with the traditional BLUP are biased in populations undergoing genomic selection and the differences between BLUP and ssGBLUP could be caused by genomic preselection unaccounted for by BLUP.

Breeding value prediction

Correlation between breeding values estimated for candidate cows with BLUP and ssGBLUP were 0.96, 0.95, and 0.95 for MY, FY, and PY, respectively. Correlations of BLUP and ssGBLUP breeding values estimated for the other animals were 0.99 for all the three considered traits. Also, the rank correlations between BLUP and ssGBLUP breeding values was high (0.98), indicating that the methods are comparable.

Correlations and prediction accuracies increased as the number of genotypes included in the analysis increased (Table 3). MY showed better values compared to the other two traits, probably because its higher h². Indeed, prediction accuracy is influenced by the trait heritability (Hayes et al. 2009). Average correlations among traits in the scenario with the largest number of genotypes available (E) were 0.11 larger in ssGBLUP compared to BLUP. The highest correlation between the two sets of BVs (i.e. stability of the model) was found with the largest number of genotypes: values of correlation were 0.83, 0.81, and 0.82 for MY, FY, and PY. Average predictions accuracies increased from 0.57 to 0.75 moving from BLUP to ssGBLUP.

Table 3. LR validation results with BLUP and single-step genomic BLUP (ssGBLUP).

		ssGBLUP^a
	BLUP	A	B	C	D	E
N genotypes	–	35	50	85	463	498
Correlation
Milk	0.72	0.72	0.75	0.77	0.82	0.83
Fat	0.71	0.70	0.75	0.76	0.81	0.81
Protein	0.69	0.69	0.73	0.75	0.82	0.82
Accuracy
Milk	0.60	0.60	0.67	0.68	0.77	0.77
Fat	0.55	0.55	0.62	0.63	0.71	0.72
Protein	0.57	0.57	0.65	0.66	0.76	0.76
Dispersion
Milk	0.93	0.93	0.92	0.95	0.89	0.90
Fat	0.98	0.98	0.99	1.01	0.92	0.93
Protein	0.92	0.92	0.93	0.96	0.92	0.92
Bias
Milk	−0.107	−0.106	−0.103	−0.099	−0.014	−0.019
Fat	−0.013	−0.013	−0.013	−0.015	−0.005	−0.006
Protein	−0.005	−0.005	−0.005	−0.005	−0.001	−0.001

^aA: genotypes only for 35 bulls; B: genotypes only for the 50 candidates; C: genotypes for 50 candidates + 35 bulls; D: genotypes for 463 cows; E: genotypes for 463 cows + 35 bulls.

BLUP: best linear unbiased prediction; ssGBLUP: single-step genomic best linear unbiased prediction.

Interesting results have been obtained for the dispersion (Table 3), even if a clear pattern across the different scenarios could not be detected. The lowest dispersion was observed for ssGBLUP in the scenario C (i.e. with genotypes for candidates and 35 bulls). Even if dispersions for ssGBLUP-E were worser than BLUP ones, they were within a 10% difference from the ideal value (1). Tsuruta et al. (2011) reported that values within 15% are acceptable. The lowest biases were found when genotypes of all cows (N = 463) were included in the analysis (i.e. ssGBLUP-D).

The differences in correlations and accuracies between the use of genotypes of only females or both females and males (i.e. scenarios B vs C and D vs E) were limited, suggesting that the number of genotypes available and not the source (males vs females) mainly affects genomic prediction performances. The inclusion of only bull genotypes (i.e. ssGBLUP-A) did not change prediction statistics in comparison with BLUP, apart from bias, probably because the sire contribution was already included in the analysis through the phenotypes of the cows.

Higher breeding values accuracies (+8%), based on the prediction error variance (Gonzáles-Récio et al. 2006), were reported for ssGBLUP compared to BLUP in Brazilian buffaloes (Aspilcueta-Borquis et al. (2015), together with lower standard deviations for ssGBLUP.

More accurate predictions were reported by genomic methods (GBLUP and ssGBLUP) compared to BLUP in Philippine dairy buffaloes (Herrera et al. 2019). These authors compared accuracies for genotyped and for non-genotyped animals and they reported the highest increase in accuracy for the ungenotyped animals (+0.07 on average among milk, fat and protein) compared to the group of genotyped animals only (+0.03 on average among milk, fat and protein).

LR method was recently widely adopted to validate breeding values estimation in almost all livestock species, such as cattle (Porto Neto et al. 2019; Cesarani, Hidalgo, et al. 2020; Durbin et al. 2020), sheep (Granado-Tajada et al. 2020; Macedo et al. 2020), pigs (Aliakbari et al. 2020), chicken (Chu et al. 2019; Bermann et al. 2021) and fish (Silva et al. 2019). Results about the comparison between BLUP and ssGBLUP, by using LR method, can change according to the considered species. No clear advantages were found for genomic selection of Latxa dairy sheep breed (Granado-Tajada et al. 2020), whereas other studies reported better performances of ssGBLUP compared to BLUP (Chu et al. 2019; Silva et al. 2019). For example, Chu et al. (2019) reported that the use of ssGBLUP increased population accuracy of BV for both genotyped and non-genotyped birds compared to the use of BLUP. However, bias of BV prediction increased for non-genotyped birds.

The results for the WssGBLUP are presented in Table 4. Overall, there was an increase in prediction accuracy and a reduction in bias which are beneficial, however at the same time there was an increase in the dispersion of the breeding values, which is undesirable. The difference between iterations (2–5) was small, with the second iteration reaching the highest accuracy. On the other hand, there was a considerable increase in accuracy when the CT values were increased. For instance, for MY prediction accuracy went from 0.79 (CT = 1.025) to 0.82 (CT = 1.25) which represent an increase of 2.6 and 6.5% respectively, over ssGBLUP.

Table 4. Weighted ssGBLUP-E using the largest number of genotypes (463 cows + 35 bulls).

CT		1.025				1.125				1.250
Round	ssGBLUP*	2	3	4	5	2	3	4	5	2	3	4	5
Correlation
Milk	0.83	0.82	0.82	0.82	0.82	0.82	0.82	0.82	0.82	0.81	0.81	0.81	0.81
Fat	0.81	0.81	0.81	0.81	0.81	0.81	0.81	0.81	0.81	0.80	0.80	0.80	0.80
Protein	0.82	0.82	0.82	0.82	0.82	0.82	0.82	0.82	0.82	0.81	0.81	0.81	0.80
Accuracy
Milk	0.77	0.79	0.79	0.79	0.79	0.80	0.80	0.80	0.80	0.81	0.82	0.82	0.82
Fat	0.72	0.73	0.73	0.73	0.73	0.75	0.75	0.75	0.75	0.76	0.77	0.77	0.77
Protein	0.76	0.77	0.77	0.77	0.77	0.79	0.79	0.79	0.79	0.80	0.81	0.81	0.81
Dispersion
Milk	0.90	0.89	0.89	0.89	0.89	0.88	0.88	0.88	0.88	0.86	0.85	0.85	0.85
Fat	0.93	0.92	0.92	0.92	0.92	0.91	0.91	0.91	0.91	0.90	0.89	0.88	0.88
Protein	0.92	0.92	0.92	0.92	0.92	0.90	0.90	0.90	0.90	0.89	0.88	0.87	0.87
Bias
Milk	−0.019	−0.0172	−0.0172	−0.0172	−0.0172	−0.0113	−0.0099	−0.0096	−0.0095	−0.0047	−0.0005	0.0023	0.0028
Fat	−0.006	−0.0060	−0.0060	−0.0060	−0.0060	−0.0057	−0.0056	−0.0055	−0.0055	−0.0052	−0.0049	−0.0048	−0.0048
Protein	−0.001	−0.0012	−0.0012	−0.0012	−0.0012	−0.0010	−0.0009	−0.0009	−0.0009	−0.0007	−0.0005	−0.0004	−0.0004

*Results from the 1st iteration of WssGBLUP are the same as unweighted ssGBLUP.

CT: carcass traits; ssGBLUP: single-step genomic best linear unbiased prediction.

The performance of the WssGBLUP compared to ssGBLUP depends on the architecture of the trait and on the sample size. For instance, Lourenco et al. (2017), showed that for highly polygenic traits there were no improvements in prediction accuracy with WssGBLUP compared to ssGBLUP but under oligogenic traits, accuracy can be improved with SNP weighting. Additionally, with larger sample sizes, the benefit of using different weights for markers may be reduced, as pointed out by Karaman et al. (2016) who showed that the prediction accuracy of Bayesian variable selection methods and GBLUP converged to the similar prediction accuracy once the sample size was large enough.

As aforementioned, the WssGBLUP is a weighted version of ssGBLUP in which different weights are used for each SNP. Giving more weights to some SNPs allows the model to take into account the presence of major genes or QTL that affect the trait of interest. This is particularly appealing for species, like buffalo, where milk protein composition or cheesemaking traits are key features. In this regard, the integration of α_s1 casein gene in genomic evaluation of dairy goats are effective examples. Dagnachew et al. (2011) analysed 38 SNPs located within the four casein genes in Norwegian goats and showed significant additive effect of a single SNP within two casein variants (i.e. CSN1S1 and CSN3 genes) on fat and protein percentages, milk yield and milk taste. Pizarro Inostroza et al. (2020) analysed the effects of considering 48 casein loci-located SNP in the Murciano-Granadina goats. Teissier et al. (2018) investigated three weighted ssGBLUP methods to integrate information on the α_s1 casein in genomic evaluations of dairy goats. All authors concluded that the inclusion of detailed information on major genes (e.g. DGAT1 for fat content or αs1 casein for protein content) with additive, dominance, and epistatic effects in the genetic evaluation could improve both the statistical power of the model and the accuracy of breeding values for milk yield and composition compared to a model without these effects.

Indeed, the slight improvement in accuracy with WssGBLUP in our study may indicate the presence of important markers associated with these traits. However, only one SNP showed a p-value of its effect passing the significance threshold (Supplementary Figure 1). The low number of genotypes involved in the present study could have been influenced the results; anyway, further studies will be needed to determine the genomic make up of milk production in Italian buffalos. Liu et al. (2018) reported 4 SNP in 2 genomic regions significantly associated with fat yield and protein percentage in water buffalo and concluded that a more precise map of the buffalo genome as well as a larger sample size is needed in future studies to be able to identify important genes affecting the traits.

It should be also pointed out that the use of WssGBLUP resulted in an increase of prediction accuracy but together with an increase of breeding value dispersion. This latter statistic should be as much as closer to 1, which means that breeding values are not under or over estimated. In our results, the dispersion was 0.90, 0.93, and 0.92 for ssGBLUP (scenario E), compared to 0.86, 0.90, and 0.89 for WssGBLUP (CT = 1.25, iteration 2) for MY, FY, and PY, respectively (Table 4). As pointed out by Legarra and Reverter (2017), the bias and the slope of the regression (dispersion) need to be taken in to account specially when proven and young animals are mixed in the population, because they could affect selection decisions, beyond the accuracy of selection. In this scenario, one would have to establish a balance between the highest gains in accuracy and the least bias and dispersion of the breeding values when choosing the model for the official evaluations.

Conclusions

This study, even if based on a small number of genotyped animals, provided interesting insights about the use and future application of ssGBLUP in the Italian Buffalo. Variance components estimated using BLUP and ssGBLUP were similar, but the latter model showed lower standard deviations for genetic correlations. The results showed that the inclusion of genotypes only for bulls (i.e. ssGBLUP-A) did not lead to any improvements compared to BLUP. On the contrary, the inclusion of genotypes for cows led to an increase in almost all validation parameters. Thus, to implement a genomic evaluation in the Italian Buffalo, it is important to emphasise that female genotyping could increase the accuracies of the prediction compared to both pedigree-based and genomic evaluations with only male genotypes models.

Acknowledgements

We gratefully acknowledge the Associazione Nazionale Allevatori Specie Bufalina (ANASB) for providing data.

Disclosure statement

No potential conflict of interest was reported by the author(s).