Advertisement

Genome-wide association studies to identify quantitative trait loci affecting milk production traits in water buffalo

  • J.J. Liu
    Affiliations
    Key Lab of Agricultural Animal Genetics, Breeding and Reproduction of Ministry of Education, Huazhong Agriculture University, Wuhan, Hubei, China 430070

    Hubei Province's Engineering Research Center in Buffalo Breeding and Products, Wuhan, Hubei, China 430070
    Search for articles by this author
  • A.X. Liang
    Affiliations
    Key Lab of Agricultural Animal Genetics, Breeding and Reproduction of Ministry of Education, Huazhong Agriculture University, Wuhan, Hubei, China 430070

    Hubei Province's Engineering Research Center in Buffalo Breeding and Products, Wuhan, Hubei, China 430070
    Search for articles by this author
  • G. Campanile
    Affiliations
    Department of Veterinary Medicine and Animal Productions, University of Naples “Federico II”, Naples, Italy 80137
    Search for articles by this author
  • G. Plastow
    Affiliations
    Department of Agricultural, Food, and Nutritional Sciences, University of Alberta, Edmonton, AB, Canada T6G 2C8
    Search for articles by this author
  • C. Zhang
    Affiliations
    Department of Agricultural, Food, and Nutritional Sciences, University of Alberta, Edmonton, AB, Canada T6G 2C8
    Search for articles by this author
  • Z. Wang
    Affiliations
    Department of Agricultural, Food, and Nutritional Sciences, University of Alberta, Edmonton, AB, Canada T6G 2C8
    Search for articles by this author
  • A. Salzano
    Affiliations
    Department of Veterinary Medicine and Animal Productions, University of Naples “Federico II”, Naples, Italy 80137
    Search for articles by this author
  • B. Gasparrini
    Affiliations
    Department of Veterinary Medicine and Animal Productions, University of Naples “Federico II”, Naples, Italy 80137
    Search for articles by this author
  • M. Cassandro
    Affiliations
    Department of Agronomy, Food, Natural Resources, Animal, and Environment, University of Padova, Agripolis, Legnaro, Italy 35020
    Search for articles by this author
  • L.G. Yang
    Correspondence
    Corresponding author
    Affiliations
    Key Lab of Agricultural Animal Genetics, Breeding and Reproduction of Ministry of Education, Huazhong Agriculture University, Wuhan, Hubei, China 430070

    Hubei Province's Engineering Research Center in Buffalo Breeding and Products, Wuhan, Hubei, China 430070
    Search for articles by this author
Open ArchivePublished:November 08, 2017DOI:https://doi.org/10.3168/jds.2017-13246

      ABSTRACT

      Water buffalo is the second largest resource of milk supply around the world, and it is well known for its distinctive milk quality in terms of fat, protein, lactose, vitamin, and mineral contents. Understanding the genetic architecture of milk production traits is important for future improvement by the buffalo breeding industry. The advance of genome-wide association studies (GWAS) provides an opportunity to identify potential genetic variants affecting important economical traits. In the present study, GWAS was performed for 489 buffaloes with 1,424 lactation records using the 90K Affymetrix Buffalo SNP Array (Affymetrix/Thermo Fisher Scientific, Santa Clara, CA). Collectively, 4 candidate single nucleotide polymorphisms (SNP) in 2 genomic regions were found to associate with buffalo milk production traits. One region affecting milk fat and protein percentage was located on the equivalent of Bos taurus autosome (BTA)3, spanning 43.3 to 43.8 Mb, which harbored the most likely candidate genes MFSD14A, SLC35A3, and PALMD. The other region on the equivalent of BTA14 at 66.5 to 67.0 Mb contained candidate genes RGS22 and VPS13B and influenced buffalo total milk yield, fat yield, and protein yield. Interestingly, both of the regions were reported to have quantitative trait loci affecting milk performance in dairy cattle. Furthermore, we suggest that buffaloes with the C allele at AX-85148558 and AX-85073877 loci and the G allele at AX-85106096 locus can be selected to improve milk fat yield in this buffalo-breeding program. Meanwhile, the G allele at AX-85063131 locus can be used as the favorable allele for improving milk protein percentage. Genomic prediction showed that the reliability of genomic estimated breeding values (GEBV) of 6 milk production traits ranged from 0.06 to 0.22, and the correlation between estimated breeding values and GEBV ranged from 0.23 to 0.35. These findings provide useful information to understand the genetic basis of buffalo milk properties and may play a role in accelerating buffalo breeding programs using genomic approaches.

      Key words

      INTRODUCTION

      Water buffalo (Bubalus bubalis) is an important livestock species for the agricultural economy, supplying milk, meat, and draft power (
      • Warriach H.M.
      • Mcgill D.M.
      • Bush R.D.
      • Wynn P.C.
      • Chohan K.R.
      A review of recent developments in buffalo reproduction - a review.
      ). The global buffalo population was recently estimated to be 194 million, 97% of which were reared in Asia (
      • FAOSTAT
      Food and Agriculture Organization the United Nations Statistics Division.
      ). Buffalo is well known for its high milk quality, with higher fat (6.4–8.0% vs. 4.1–5.0%) and protein (4.0–4.5% vs. 3.4–3.6%) contents than cow milk (
      • Khedkar C.
      • Kalyankar S.
      • Deosarkar S.
      Buffalo milk.
      ). Its compositional and functional properties make buffalo milk suitable for manufacture of dairy products, such as superior cream, butter, yogurt, and cheese, especially mozzarella cheese (
      • Michelizzi V.N.
      • Dodson M.V.
      • Pan Z.
      • Amaral M.E.J.
      • Michal J.J.
      • McLean D.J.
      • Womack J.E.
      • Jiang Z.
      Water buffalo genome science comes of age.
      ). As one of the most famous dairy buffalo breeds, the Italian Mediterranean buffalo has reached a high productivity standard due to the intense work of selection and study by the Italian buffalo-breeding program. Italian Mediterranean buffaloes are of the river type, whereas in China and some southeast Asian countries, most buffalo breeds belong to the swamp type and have poor milk production. To improve buffalo milk performance, a crossbreeding strategy is often used by the traditional breeding industry. Although a remarkable improvement has been achieved over the years, milk production in buffalo is still considerably lower than in cow. Presently, buffalo milk is ranked second for the world's total milk production but accounts for only about 13% of the total (
      • FAOSTAT
      Food and Agriculture Organization the United Nations Statistics Division.
      ). Therefore, understanding the genetic architecture of milk properties is essential to accelerate the genetic improvement in water buffalo-breeding programs.
      As is well known, QTL can be utilized to identify candidate genes and contribute to the dissection of genetic mechanisms underlying economic traits in animals. A large number of QTL have been detected for milk production traits in cattle, such as DGAT1, ABCG2, and SCD1 genes (
      • Lengi A.J.
      • Corl B.A.
      Identification and characterization of a novel bovine stearoyl-CoA desaturase isoform with homology to human SCD5.
      ;
      • Weller J.I.
      • Ron M.
      Invited review: quantitative trait nucleotide determination in the era of genomic selection.
      ). Polymorphisms for DGAT1 and ABCG2 were detected in buffalo although the alleles were found to be fixed in some breeds (
      • Tantia M.S.
      • Vijh R.K.
      • Mishra B.P.
      • Mishra B.
      • Kumar S.B.
      • Sodhi M.
      DGAT1 and ABCG2 polymorphism in Indian cattle (Bos indicus) and buffalo (Bubalus bubalis) breeds.
      ;
      • Shi D.S.
      • Wang J.
      • Yang Y.
      • Lu F.H.
      • Li X.P.
      • Liu Q.Y.
      DGAT1, GH, GHR, PRL and PRLR polymorphism in water buffalo (Bubalus bubalis).
      ). This implies that these 2 genes may be responsible for the high milk fat and milk protein in buffalo. To more precisely identify markers and genomic regions associated with quantitative traits, genome-wide association studies (GWAS) are increasingly used and successfully incorporated into dairy cattle (
      • Hayes B.J.
      • Bowman P.
      • Chamberlain A.
      • Goddard M.
      Invited review: Genomic selection in dairy cattle: Progress and challenges.
      ), pig (
      • Do D.N.
      • Janss L.L.G.
      • Jensen J.
      • Kadarmideen H.N.
      SNP annotation-based whole genomic prediction and selection: An application to feed efficiency and its component traits in pigs.
      ), and poultry (
      • Fulton J.E.
      Genomic selection for poultry breeding.
      ) breeding programs. However, genomic research on buffalo is still very limited. Without an available complete buffalo genomic reference map, buffalo genomic studies use information from the most homologous species available—cattle (
      • Amaral M.E.J.
      • Grant J.R.
      • Riggs P.K.
      • Stafuzza N.B.
      • Edson Filho A.
      • Goldammer T.
      • Weikard R.
      • Brunner R.M.
      • Kochan K.J.
      • Greco A.J.
      A first generation whole genome RH map of the river buffalo with comparison to domestic cattle.
      ;
      • Di Meo G.P.
      • Perucatti A.
      • Floriot S.
      • Hayes H.
      • Schibler L.
      • Incarnato D.
      • Di Berardino D.
      • Williams J.
      • Cribiu E.
      • Eggen A.
      An extended river buffalo (Bubalus bubalis, 2n= 50) cytogenetic map: assignment of 68 autosomal loci by FISH-mapping and R-banding and comparison with human chromosomes.
      ;
      • Venturini G.C.
      • Cardoso D.
      • Baldi F.
      • Freitas A.
      • Aspilcueta-Borquis R.
      • Santos D.
      • Tonhati H.
      Association between single-nucleotide polymorphisms and milk production traits in buffalo.
      ). In our previous study (
      • Wu J.J.
      • Song L.J.
      • Wu F.J.
      • Liang X.W.
      • Yang B.Z.
      • Wathes D.C.
      • Pollott G.E.
      • Cheng Z.
      • Shi D.S.
      • Liu Q.Y.
      Investigation of transferability of BovineSNP50 BeadChip from cattle to water buffalo for genome wide association study.
      , we investigated the transferability of Bovine SNP50 BeadChip (Illumina Inc., San Diego, CA) data from cattle to buffalo for GWAS and identified 7 SNP among 935 bovine SNP associating with buffalo milk performance. The development of a buffalo genotyping array opens new opportunities to explore key genes regulating buffalo milk properties and provides the possibility of utilizing genomic selection in the buffalo breeding industry (
      • Iamartino D.
      • Williams J.L.
      • Sonstegard T.
      • Reecy J.
      • v. Tassell C.
      • Nicolazzi E.L.
      • Biffani S.
      • Biscarini F.
      • Schroeder S.
      • de Oliveira D.A.
      The buffalo genome and the application of genomics in animal management and improvement.
      ). With the buffalo SNP array, 78 SNP (
      • Iamartino D.
      • Williams J.L.
      • Sonstegard T.
      • Reecy J.
      • v. Tassell C.
      • Nicolazzi E.L.
      • Biffani S.
      • Biscarini F.
      • Schroeder S.
      • de Oliveira D.A.
      The buffalo genome and the application of genomics in animal management and improvement.
      ;
      • de Camargo G.M.
      • Aspilcueta-Borquis R.
      • Fortes M.
      • Porto-Neto R.
      • Cardoso D.
      • Santos D.
      • Lehnert S.
      • Reverter A.
      • Moore S.
      • Tonhati H.
      Prospecting major genes in dairy buffaloes.
      ;

      El-Halawany, N., H. Abdel-Shafy, A. E. M. A. Shawky, M. A. Abdel-Latif, A. F. M. Al-Tohamy, and O. M. A. El-Moneim. 2015. Genome-wide association study for milk production in Egyptian buffalo. Pages 10–16 in Proc. Int. Symp. Animal Functional Genomics. The 6th International Symposium on Animal Functional Genomics, Piacenza, Italy.

      ) have been found to influence milk production traits among different buffalo breeds. Furthermore, several suggestive genomic regions on BTA1, 5, 6, and 27 have been found to be associated with daily milk yield using GWAS in Egyptian buffalo (
      • El-Halawany N.
      • Abdel-Shafy H.
      • Abd-El-Monsif A.S.
      • Abdel-Latif M.A.
      • Al-Tohamy A.F.
      • El-Moneim O.M.A.
      Genome-wide association study for milk production in Egyptian buffalo.
      ). Therefore, we hypothesize that a variety of novel genes and QTL may be identified to help with the genetic dissection of buffalo milk performance.
      The present study aimed to detect important markers and genomic regions affecting milk production traits, and to investigate the feasibility of genomic selection as a potential selection strategy in buffalo breeding programs for accelerating the genetic improvement of buffalo economic traits.

      MATERIALS AND METHODS

      Ethics Statement

      The data and samples of buffalo used in this study were provided by The Italian Buffalo Breeders Association (ANASB), which is responsible for the official herd book of buffalo population in Italy. The experimental design and animal treatments were approval by the Ethical Animal Care and Use Committee of University of Naples “Federico II.” Moreover, the farmers were previously informed and in agreement with purpose and methods used.

      Animal Resources and Phenotypic Data

      A total of 1,424 lactation records were collected from 489 Italian Mediterranean buffaloes born from 2000 to 2011 and reared in 4 herds in southern Italy. Pedigree data consisted of 937 animals over 3 generations. Six milk production traits, peak milk yield (PM), total milk yield (MY), fat yield (FY), fat percentage (FP), protein yield (PY), and protein percentage (PP), were recorded. All milk production traits were adjusted to 270 d in milk (
      • Baldi F.
      • Laureano M.M.M.
      • Gordo D.G.M.
      • Bignardi A.B.
      • Borquis R.R.A.
      • Albuquerque L.G.d.
      • Tonhati H.
      Effect of lactation length adjustment procedures on genetic parameter estimates for buffalo milk yield.
      ). The linear model used to adjust the records, including the factors herd-season (HS, 4 farms and 2 seasons), year of calving (<2005 and 2005–2014), parity (1 to 7 and ≥8), and calf sex (male and female), were tested through a fixed linear model. Only significant (P < 0.05) factors were included in the model as fixed effects to adjust the records of lactation length (LL) between 150 and 270 d using the LSM method. The records for LL >270 d were truncated at the 270-d milk production. Data from LL <150 d were excluded from this analysis. The adjusted formula for MY, FY, and PY was
      Y270=Yn×LSM270LSMn,


      where Y270 is the 270-d adjusted phenotypes MY270, FY270, and PY270; Yn is the observed phenotype at day n; LSM270 and LSMn are the least squares means of the observed phenotypes at d 270 and day n; n (150 < n < 270) is the days of LL. Then, FP270 and PP270 were adjusted as follows:
      FP270=FY270MY270andPP270=PY270MY270.


      Genotyping and Quality Control

      Genomic DNA was extracted from whole blood using a standard phenol-chloroform extraction protocol. Genotyping was conducted at Delta Genomics (Edmonton AB, Canada) using the 90K Axiom Buffalo SNP Array (Affymetrix/Thermo Fisher Scientific, Santa Clara, CA). Quality control (QC) for genotypes was performed using PLINK1.9 () for individuals with a call rate ≥97%. The SNP were selected with the criteria of genotyping call rate ≥95%, minor allele frequency ≥0.05, P-value of chi-squared test of Hardy-Weinberg equilibrium ≥10−5, and with map information.

      Breeding Values and Estimation of Genetic and Phenotypic Parameters

      The EBV for the 6 milk production traits (PM, MY270, FY270, PY270, FP270, and PP270) were estimated with a univariate animal model suing ASReml 3.0 (
      • Gilmour A.R.
      • Gogel B.
      • Cullis B.
      • Thompson R.
      • Butler D.
      ASReml user guide release 3.0.
      ):
      y=Xb+Z1a+Z2p+e,


      where y is a vector of phenotypes adjusted to 270-d records for all traits except PM; X is a incidence matrix associated with the fixed effects; b is a fixed vector containing HS, calving year, parity, and calf sex for 270-d MY270, FY270, and PY270, and containing HS, calving year, and parity for FP270, PP270, and PM traits; Z1 is a incidence matrix associated with individual additive genetic effects; a is a random vector of animal's breeding values or additive genetic effects [a(0,Aσa2)]; Z2 is a incidence matrix associated with animals' permanent environmental effects; p is a vector of permanent environmental effects of individual animal [p(0,Iσp2)]; and e is a random vector of residual errors [e(0,Iσe2)]. The σ2a, σ2p, and σ2e terms are the additive genetic, permanent environment, and residual error variances, respectively; and A and I are the additive genetic relationship and identity matrices, respectively.
      The genetic and phenotypic parameters, including heritability and genetic and phenotypic correlations, were estimated with the following pairwise bivariate animal model using ASReml 3.0 (
      • Gilmour A.R.
      • Gogel B.
      • Cullis B.
      • Thompson R.
      • Butler D.
      ASReml user guide release 3.0.
      ):
      [y1y2]=[X100X2][b1b2]+[Z100Z2][a1a2]+[Z300Z4][p1p2]+[e1e2],


      where y1 and y2 are the vectors of phenotypes for traits 1 and 2, respectively; X1 and X2 are the incidence matrices associating the fixed effects to vectors y1 and y2, respectively; b1 and b2 are the vectors of fixed effects for traits 1 and 2, respectively; Z1 and Z2 are the incidence matrices associating the random additive genetic effect of a1 and a2 to vectors of y1 and y2, respectively; Z3 and Z4 are the incidence matrices associating the random permanent effects p1 and p2 to vectors of y1 and y2, respectively; and e1 and e2 are the random vectors of residual errors for traits 1 and 2, respectively. The variance components definition for a, p, and e remained as the same as in the univariate model. The heritability was an averaged estimates using variance components obtained from the corresponding pairwise bivariate analyses (
      • Miar Y.
      • Plastow G.
      • Moore S.
      • Manafiazar G.
      • Charagu P.
      • Kemp R.
      • Van Haandel B.
      • Huisman A.
      • Zhang C.
      • McKay R.
      Genetic and phenotypic parameters for carcass and meat quality traits in commercial crossbred pigs.
      ).

      Genome-Wide Association Analyses

      Deregressed EBV (DEBV) were calculated according to the method reported by
      • Garrick D.J.
      • Taylor J.F.
      • Fernando R.L.
      Deregressing estimated breeding values and weighting information for genomic regression analyses.
      . The deregression procedure was performed by weighted EBV with the weight wi=(1h2)/{[c+(1ri2)/ri2]h2}, where c was the part of genetic variance not explained by markers, which was assumed to be 0.4 (
      • Saatchi M.
      • Schnabel R.D.
      • Taylor J.F.
      • Garrick D.J.
      Large-effect pleiotropic or closely linked QTL segregate within and across ten US cattle breeds.
      ); h2 was the heritability of the trait; and r2i was the reliability of EBV of the ith animal. The DEBV were considered a new phenotype for subsequent genomic association and prediction analyses.
      A total of 412 animals with both DEBV and genotype were used to perform the genome-wide association analysis using the ridge regression BLUP (rrBLUP) model in R (
      • Endelman J.B.
      Ridge regression and other kernels for genomic selection with R package rrBLUP.
      ):
      y=Xβ+u+e,


      where y is the vector of DEBV, X is a vector of coded SNP genotypes (−1, 0, or 1), β is the fixed additive genetic value attributed to SNP under evaluation, u is the vector of the background polygenic effects with normal distribution [u(0,Kσu2)], where K is the genomic relationship matrix obtained from pedigree information (
      • Endelman J.B.
      Ridge regression and other kernels for genomic selection with R package rrBLUP.
      ) and σ2u is the genetic variance, and e is the vector of residual errors with normal distribution [e(0,Iσe2)], where I is the identity matrix and σ2e is the residual variance. False discovery rate (FDR) was adopted to adjust P-value for all detected SNP, and the genome-wide significant threshold of SNP was defined as FDR <0.10 and the suggestive threshold as P < 10−4. The LSM of EBV for the 3 genotypes affecting milk production traits of significant SNP were calculated by a general liner model using R package lsmeans (
      • Lenth R.V.
      Least squares means: The R package lsmeans.
      ), and the significant threshold was set at P < 0.05.

      Haplotype Analyses

      The 0.5-Mb genomic window around significant SNP was identified to construct the haplotype blocks using Haploview 4.2 (
      • Barrett J.C.
      • Fry B.
      • Maller J.
      • Daly M.J.
      Haploview: Analysis and visualization of LD and haplotype maps.
      ). Association between each haplotype combination and 6 milk production traits were evaluated with Bonferroni t-test in SAS 9.4 (SAS Institute Inc., Cary, NC) as described by
      • Yang S.-H.
      • Bi X.-J.
      • Xie Y.
      • Li C.
      • Zhang S.-L.
      • Zhang Q.
      • Sun D.-X.
      Validation of PDE9A gene identified in GWAS showing strong association with milk production traits in Chinese Holstein.
      . The genomic region that covered haplotype blocks affecting (P < 0.05) buffalo milk production traits was identified as a candidate region.

      Bioinformatics Analyses of Candidate Genes

      Because the 90K Affymetrix Buffalo SNP Array was aligned to the bovine genome, identification of genes within candidate regions was based on the Bos taurus UMD3.1 genomic assembly (http://bovinegenome.org/?q=node/61). Genes were submitted to the database Bovinemine (http://bovinegenome.org/bovinemine) for pathway and gene ontology (GO) term annotation. The functions and related information about the genes were summarized using the available database of GeneCards (http://www.genecards.org/).

      Genomic Prediction

      A 5-fold cross validation was used to evaluate the accuracy of predicted genomic EBV (GEBV;
      • Lukić B.
      • Pong-Wong R.
      • Rowe S.
      • Koning D.
      • Velander I.
      • Haley C.
      • Archibald A.
      • Woolliams J.
      Efficiency of genomic prediction for boar taint reduction in Danish Landrace pigs.
      ). Briefly, the 412 buffaloes were randomly grouped into 5 folds. Each time, 4 folds (330 buffaloes) were set as the training population and the remaining fold (82 buffaloes) was used as the validation population. For each validation, the EBV and reliability were used as pseudo-phenotypes in a genomic BLUP model, and GEBV and its reliability were evaluated through the genomic relationship matrix using the gebv software (
      • Sargolzaei M.
      • Schenkel F.S.
      • Vanraden P.M.
      gebv: Genomic breeding value estimator for livestock. Technical report to the Dairy Cattle Breeding and Genetics Committee.
      ). Both the reliability of GEBV and the correlation between EBV and GEBV were used to evaluate the prediction accuracy (
      • Moser G.
      • Tier B.
      • Crump R.E.
      • Khatkar M.S.
      • Raadsma H.W.
      A comparison of five methods to predict genomic breeding values of dairy bulls from genome-wide SNP markers.
      ). The average reliability value of the 5 cross-fold validation was finally used to evaluate the feasibility of genomic prediction.

      RESULTS AND DISCUSSION

      Phenotype and Genotype Data Description

      Italian Mediterranean buffalo lactation usually lasts around 270 d. Hence, we adjusted the buffalo milk production traits for 270 d to estimate more reliable EBV (
      • Baldi F.
      • Laureano M.M.M.
      • Gordo D.G.M.
      • Bignardi A.B.
      • Borquis R.R.A.
      • Albuquerque L.G.d.
      • Tonhati H.
      Effect of lactation length adjustment procedures on genetic parameter estimates for buffalo milk yield.
      ). Six 270-d adjusted traits including PM, MY270, FY270, FP270, PY270, and PP270 were used for GWAS. In Table 1, the phenotypic descriptive statistics of analyzed samples are reported. Overall, the Italian Mediterranean buffalo in this study can be considered a highly selected breed with a higher milk production yield (2,965.6 ± 516.1 kg) than that reported in Murrah (1,712.5 kg;
      • Tonhati H.
      • Muñoz M.
      • Duarte J.
      • Reichert R.
      • Oliveira J.
      • Lima A.
      Estimates of correction factors for lactation length and genetic parameters for milk yield in buffaloes.
      ) and Nili-Ravi (1,984.4 kg;
      • Khan M.
      • Hyder A.
      • Bajwa I.
      • Rehman M.
      • Hassan F.
      Prediction of lactation yield from last-record day and average daily yield in Nili-Ravi buffaloes.
      ) buffaloes. The peak milk yield of total lactation showed an average of 14.9 ± 2.6 kg. Buffalo milk is usually rich in fat and protein (
      • Rosati A.
      • Van Vleck L.D.
      Estimation of genetic parameters for milk, fat, protein and mozzarella cheese production for the Italian river buffalo Bubalus bubalis population.
      ). The average fat and protein percentages in our population were 8.3 and 4.7%, respectively, similar to that of Murrah and Nili-Ravi buffalo (
      • Tonhati H.
      • Muñoz M.
      • Duarte J.
      • Reichert R.
      • Oliveira J.
      • Lima A.
      Estimates of correction factors for lactation length and genetic parameters for milk yield in buffaloes.
      ;
      • Khan M.
      • Hyder A.
      • Bajwa I.
      • Rehman M.
      • Hassan F.
      Prediction of lactation yield from last-record day and average daily yield in Nili-Ravi buffaloes.
      ) but higher than that reported by
      • Khedkar C.
      • Kalyankar S.
      • Deosarkar S.
      Buffalo milk.
      . The buffalo milk fat and protein yields were 244.9 ± 47.3 and 137.5 ± 22.9 kg, respectively. In the future, it could be interesting to assess coagulation ability of buffalo milk to improve cheese yield genetically.
      Table 1Descriptive statistics of the 270-d adjusted milk production traits in water buffalo
      Trait
      PM = peak milk yield; MY270 = 270-d total milk yield; FY270 = 270-d fat yield; FP270 = 270-d fat percentage; PY270 = 270-d protein yield; PP270 = 270-d protein percentage.
      nMeanSDMinimumMaximumCV (%)
      PM (kg)1,40814.92.66.628.017.7
      MY270 (kg)1,4082,965.6516.11,394.05,160.017.4
      FY270 (kg)1,401244.947.3102.042919.3
      FP270 (%)1,4018.30.95.511.611.2
      PY270 (kg)1,408137.522.963.0259.016.7
      PP270 (%)1,4084.70.33.85.55.6
      1 PM = peak milk yield; MY270 = 270-d total milk yield; FY270 = 270-d fat yield; FP270 = 270-d fat percentage; PY270 = 270-d protein yield; PP270 = 270-d protein percentage.
      The estimated genetic parameters for each trait are shown in Table 2. We found the heritability of buffalo milk production traits to be moderate (0.19–0.38), and the h2 values for FY270 (0.35) and PP270 (0.38) were higher than those of other traits (0.19–0.33). In agreement with previous studies (
      • Rosati A.
      • Van Vleck L.D.
      Estimation of genetic parameters for milk, fat, protein and mozzarella cheese production for the Italian river buffalo Bubalus bubalis population.
      ;
      • de Camargo G.M.
      • Aspilcueta-Borquis R.
      • Fortes M.
      • Porto-Neto R.
      • Cardoso D.
      • Santos D.
      • Lehnert S.
      • Reverter A.
      • Moore S.
      • Tonhati H.
      Prospecting major genes in dairy buffaloes.
      ), we also detected high and positive genetic correlations among PM, MY270, FY270, and PY270, ranging from 0.80 to 0.98, and a moderate genetic correlation between FP270 and PP270 (0.58). As the genetic parameters of studied traits underpin further GWAS (
      • de Camargo G.M.
      • Aspilcueta-Borquis R.
      • Fortes M.
      • Porto-Neto R.
      • Cardoso D.
      • Santos D.
      • Lehnert S.
      • Reverter A.
      • Moore S.
      • Tonhati H.
      Prospecting major genes in dairy buffaloes.
      ), there would be more power to detect SNP associated with FY270 and PP270 because they have higher heritability values.
      Table 2Estimates of heritability (bold, in diagonal), genetic correlation (above diagonal), and phenotypic correlation (below diagonal) among the studied traits
      Trait
      PM = peak milk yield; MY270 = 270-d total milk yield; FY270 = 270-d fat yield; FP270 = 270-d fat percentage; PY270 = 270-d protein yield; PP270 = 270-d protein percentage.
      PMMY270FY270FP270PY270PP270
      PM0.280.910.81−0.330.81−0.31
      MY2700.880.330.82−0.460.95−0.49
      FY2700.770.830.350.490.980.29
      FP270−0.25−0.280.400.270.130.58
      PY2700.840.960.87−0.100.190.17
      PP270−0.16−0.220.150.500.150.38
      1 PM = peak milk yield; MY270 = 270-d total milk yield; FY270 = 270-d fat yield; FP270 = 270-d fat percentage; PY270 = 270-d protein yield; PP270 = 270-d protein percentage.
      We filtered genotype data under a QC scenario, and several 462 individuals and 60,387 SNP remained. The SNP that passed QC were distributed uniformly across bovine chromosomes based on the bovine genome (Supplemental Table S1; https://doi.org/10.3168/jds.2017-13246) with an average density of 23 SNP/Mb on each chromosome, which provided powerful information for further GWAS and genomic prediction.

      GWAS

      Because of the special phenotypic data structure due to unbalanced repeated lactation records, DEBV of buffalo individuals were calculated and selected as a new response variable for GWAS analysis. The DEBV make good use of available information from genotyped animals as well as from their relatives, which can appropriately avoid bias introduced by simply pooling or averaging data information and account for heterogeneous variance (
      • Garrick D.J.
      • Taylor J.F.
      • Fernando R.L.
      Deregressing estimated breeding values and weighting information for genomic regression analyses.
      ). Moreover, use of DEBV is needed to generate a higher reliability of genomic breeding values than EBV and are widely accepted in genomic studies (
      • Ostersen T.
      • Christensen O.F.
      • Henryon M.
      • Nielsen B.
      • Su G.
      • Madsen P.
      Deregressed EBV as the response variable yield more reliable genomic predictions than traditional EBV in pure-bred pigs.
      ). With single marker and single trait association between DEBV and SNP for each trait (Figure 1), we identified 26 suggestive SNP (P < 10−4) associated with at least 1 of the 6 milk production traits (Supplemental Table S2; https://doi.org/10.3168/jds.2017-13246). After FDR adjustment, 4 SNP (AX-85148558, FDR = 0.006; AX-85106096, FDR = 0.006; AX-85073877, FDR = 0.08; AX-85063131, FDR = 0.09) were significantly associated with buffalo milk fat yield and protein percentage. According to the map information, the 4 SNP were located on BTA3 and BTA14, homologous to buffalo chromosomes 6 (BBU6) and 15 (BBU15), respectively (
      • El Nahas S.M.
      • de Hondt H.A.
      • Womack J.E.
      Current status of the river buffalo (Bubalus bubalis L.) gene map.
      ).
      Figure thumbnail gr1
      Figure 1Manhattan plots for genome-wide association studies of 6 milk production traits in water buffalo: PM = peak milk yield; MY270 = 270-d total milk yield; FY270 = 270-d fat yield; FP270 = 270-d fat percentage; PY270 = 270-d protein yield; PP270 = 270-d protein percentage. On the y-axis are −log10 (P-values) and the horizontal line represents P = 1.0 × 10−5. On the x-axis are the physical positions of SNP by chromosome based on the Bos taurus UMD3.1 genome assembly (http://bovinegenome.org/?q=node/61). Color version available online.
      For the significant SNP, we calculated the LSM of EBV for the 3 genotypes affecting the trait to investigate their genetic contribution shown in Figure 2. Three loci, AX-85148558, AX-85106096, and AX-85073877, were associated with PP270, and the individuals with CC or GG homozygous genotypes showed higher (P < 0.01 or P < 0.05) protein percentage among all 3 genotypes. Hence, we suggest that buffaloes with the C allele at AX-85148558 and AX-85073877 and the G allele at AX-85106096 could be selected to improve protein percentage in buffalo milk. Meanwhile, AX-85063131, affecting FY270, showed that the individuals with the GG genotype had higher (P < 0.01 or P < 0.05) fat yield compared with animals with AA or AG genotypes. Therefore, the G allele at the AX-85063131 locus could be selected as a favorable allele to improve milk fat yield in this buffalo-breeding program. The low frequency of the favorable allele provides more opportunity for improvement. However, it is also possible that the significant association was due to the small number of observations for the favorable allele. Compared with dairy cattle, dairy buffalo selection breeding programs started only recently. The lack of selection pressure in the buffalo industry and different buffalo breeds may also contribute to the low frequency of the favorable allele and genotype. We therefore recommend verifying the favorable allele effect for milk production traits in a larger sample size or by testing them in an independent buffalo population.
      Figure thumbnail gr2
      Figure 2Least squares means (and 95% CI) of EBV for the 3 genotypes affecting the trait of significant SNP (false discovery rate <0.10) detected from genome-wide association studies. PP270 = 270-d protein percentage; FY270 = 270-d fat yield. Individuals with different genotypes are shown for each SNP showing significant (P < 0.05 or P < 0.01) association with PP270 or FY270. Numbers under each genotype represent the total number of animals with that genotype; the SNP position under each SNP is based on the Bos taurus UMD 3.1 genome assembly (http://bovinegenome.org/?q=node/61).
      In previous studies, several suggestive SNP (P < 10−4) were identified and showed associations with milk production traits in Mediterranean buffalo (9 SNP;
      • Iamartino D.
      • Williams J.L.
      • Sonstegard T.
      • Reecy J.
      • v. Tassell C.
      • Nicolazzi E.L.
      • Biffani S.
      • Biscarini F.
      • Schroeder S.
      • de Oliveira D.A.
      The buffalo genome and the application of genomics in animal management and improvement.
      ), Murrah buffalo (22 SNP;
      • de Camargo G.M.
      • Aspilcueta-Borquis R.
      • Fortes M.
      • Porto-Neto R.
      • Cardoso D.
      • Santos D.
      • Lehnert S.
      • Reverter A.
      • Moore S.
      • Tonhati H.
      Prospecting major genes in dairy buffaloes.
      ), and Egyptian buffalo (8 SNP;
      • El-Halawany N.
      • Abdel-Shafy H.
      • Abd-El-Monsif A.S.
      • Abdel-Latif M.A.
      • Al-Tohamy A.F.
      • El-Moneim O.M.A.
      Genome-wide association study for milk production in Egyptian buffalo.
      ) populations. Unfortunately, overlapping SNP were not found among the 4 GWAS studies. The most likely reasons are the differences of genetic composition of breeds, different selection pressure in the populations, varying environmental conditions, inbreeding, effective population size, and allelic frequencies (
      • El-Halawany N.
      • Abdel-Shafy H.
      • Abd-El-Monsif A.S.
      • Abdel-Latif M.A.
      • Al-Tohamy A.F.
      • El-Moneim O.M.A.
      Genome-wide association study for milk production in Egyptian buffalo.
      ). In addition, buffalo SNP arrays are not specific for different species of buffalo. Furthermore, we adjusted all milk production traits to 270 d to account for different lactation lengths as a selection criterion (
      • Rosati A.
      • Van Vleck L.D.
      Estimation of genetic parameters for milk, fat, protein and mozzarella cheese production for the Italian river buffalo Bubalus bubalis population.
      ;
      • Tonhati H.
      • Cerón-Muñoz M.F.
      • Oliveira J.A.d.
      • El Faro L.
      • Lima A.L.F.
      • Albuquerque L.G.d.
      Test-day milk yield as a selection criterion for dairy buffaloes (Bubalus bubalis Artiodactyla, Bovidae).
      ) and used multiple testing (
      • Genovese C.R.
      • Lazar N.A.
      • Nichols T.
      Thresholding of statistical maps in functional neuroimaging using the false discovery rate.
      ) to control the FDR of significant SNP to make them more conservative.

      Candidate Genomic Regions

      To further detect QTL associated with milk production traits, we used a 0.5-Mb window around 4 significant SNP to determine the linkage disequilibrium (LD) relationships (Figure 3). Three identified SNP AX-85148558, AX-85106096, and AX-85073877 for PP270 were located on BTA3 at 43 Mb, and within the detected window, from 43,084,940 to 43,795,880 bp, 2 haplotype blocks were recognized. The first block (3_block1), spanning 149 kb, involved 5 SNP, and the second block (3_block2), spanning 14 kb, consisted of only 2 SNP. Moreover, 3 SNP with almost complete LD (0.97–1.00) were harbored in 3_block1. For AX-85063131, which had an effect on FY270, the region from 66,462,566 to 66,962,566 bp of BTA14 was identified, and 2 adjacent haplotype blocks were built among 13 SNP. The SNP AX-85063131 showed strong LD (0.96–1.00) with another 4 SNP to construct the first haplotype block on BTA14 (14_block1). The second block (14_block2), with a length of 219 kb, covered 8 SNP based on the 90K buffalo SNP panel.
      Figure thumbnail gr3
      Figure 3Haplotype block patterns for the significant SNP (false discovery rate <0.10) based on linkage disequilibrium (LD) within detected regions. The numbers on the top indicate the SNP order in the region; SNP in bold are significant SNP detected from genome-wide association studies. The SNP grouped in each triangle box mean they are grouped in one block based on LD (squared correlation coefficient, r2). Color version available online.
      For each block, we performed association between each haplotype and 6 buffalo milk production traits (Table 3). Compared with the single-marker association, analysis of haplotype association is more powerful for detecting significant effects on traits (
      • Yang S.-H.
      • Bi X.-J.
      • Xie Y.
      • Li C.
      • Zhang S.-L.
      • Zhang Q.
      • Sun D.-X.
      Validation of PDE9A gene identified in GWAS showing strong association with milk production traits in Chinese Holstein.
      ). Hence, 2 blocks within the 3_43 region were shown to be associated with buffalo milk protein percentage (P = 0.001 and 0.003) as well as fat percentage (P = 0.001 and 0.03). Two adjacent blocks on the 14_66 region were found to have an effect on fat yield (P = 0.002 and 0.0006), protein yield (P = 0.02 and 0.0004), and total milk yield (P = 0.007 and 0.0005) in buffalo. The results are consistent with the positive genetic correlations among these traits. According to GWAS and haplotype analyses, we identified 2 regions at 43.3–43.8 Mb on BTA3 (3_43) and 66.5–67.0 Mb on BTA14 (14_66) as QTL affecting milk production traits of water buffalo. Interestingly, the region on BTA3 at 43.29 Mb was also identified (
      • Hu Z.-L.
      • Park C.A.
      • Reecy J.M.
      Developmental progress and current status of the Animal QTLdb.
      ) as a bovine QTL affecting milk fat yield, fat percentage, protein yield, and protein percentage. From the BovineMine database (
      • Elsik C.G.
      • Unni D.R.
      • Diesh C.M.
      • Tayal A.
      • Emery M.L.
      • Nguyen H.N.
      • Hagen D.E.
      Bovine Genome Database: New tools for gleaning function from the Bos taurus genome.
      ), we discovered 2 milk protein yield–associated QTL on BTA14 at 66.5 and 66.9 Mb, respectively, which overlapped with the 14_66 region in the present study. In addition, previous studies have demonstrated that the region of 66.02–66.15 Mb on BTA14 was linked to a QTL mainly affecting milk fat yield in dairy cattle (
      • Harder B.
      • Bennewitz J.
      • Reinsch N.
      • Thaller G.
      • Thomsen H.
      • Kühn C.
      • Schwerin M.
      • Erhardt G.
      • Förster M.
      • Reinhardt F.
      Mapping of quantitative trait loci for lactation persistency traits in German Holstein dairy cattle.
      ;
      • Wibowo T.A.
      • Gaskins C.T.
      • Newberry R.C.
      • Thorgaard G.H.
      • Michal J.J.
      • Jiang Z.
      Genome assembly anchored QTL map of bovine chromosome 14.
      ).
      • de Camargo G.M.
      • Aspilcueta-Borquis R.
      • Fortes M.
      • Porto-Neto R.
      • Cardoso D.
      • Santos D.
      • Lehnert S.
      • Reverter A.
      • Moore S.
      • Tonhati H.
      Prospecting major genes in dairy buffaloes.
      detected 1 SNP (AX-85154407) on BTA14 close to the 14_66 region associated with buffalo milk production in Murrah buffalo. From GWAS analyses,
      • El-Halawany N.
      • Abdel-Shafy H.
      • Abd-El-Monsif A.S.
      • Abdel-Latif M.A.
      • Al-Tohamy A.F.
      • El-Moneim O.M.A.
      Genome-wide association study for milk production in Egyptian buffalo.
      discovered several suggestive genomic regions on BTA1, 5, 6, and 27 affecting daily milk yield in 94 Egyptian buffaloes, which have also been recognized in some dairy cattle breeds. The identification of same genomic regions in both cattle and buffalo using different experimental designs and analysis methods increase the confidence that 3_43 and 14_66 could be promising QTL associated with milk performance.
      Table 3Haplotype association analyses for 6 milk production traits in water buffalo
      Block (BTA_block)HaplotypeFrequency
      Frequency of individuals of each haplotype among population.
      (%)
      Trait
      PM = peak milk yield; MY270 = 270-d total milk yield; FY270 = 270-d fat yield; FP270 = 270-d fat percentage; PY270 = 270-d protein yield; PP270 = 270-d protein percentage.
      PMMY270FY270FP270PY270PP270
      3_block1TGTGT41.940.690.210.360.0001
      P < 0.001: significant association after Bonferroni multiple test.
      0.560.0001
      P < 0.001: significant association after Bonferroni multiple test.
      TATAT30.28
      CGGGC19.85
      3_block2CG59.820.330.050.270.03
      P < 0.05
      0.050.0026
      P < 0.001: significant association after Bonferroni multiple test.
      AA21.96
      CA18.09
      14_block1GGACT40.280.210.0071
      P < 0.01
      0.0022
      P < 0.01
      0.690.02
      P < 0.05
      0.57
      GAGAC18.62
      AAGAT18.00
      GAGAT10.41
      14_block2AGATGAAC27.850.04
      P < 0.05
      0.0005
      P < 0.001: significant association after Bonferroni multiple test.
      0.0006
      P < 0.001: significant association after Bonferroni multiple test.
      0.990.0004
      P < 0.001: significant association after Bonferroni multiple test.
      0.73
      GGATGAAC27.59
      GGGTAAGC21.16
      GGATGAAT7.55
      1 Frequency of individuals of each haplotype among population.
      2 PM = peak milk yield; MY270 = 270-d total milk yield; FY270 = 270-d fat yield; FP270 = 270-d fat percentage; PY270 = 270-d protein yield; PP270 = 270-d protein percentage.
      * P < 0.05
      ** P < 0.01
      *** P < 0.001: significant association after Bonferroni multiple test.

      Candidate Genes Within Identified Genomic Regions

      Within the candidate region, we discovered a total of 13 genes and 1 small nucleolar RNA, as listed in Table 4. Among the genes, MFSD14A, SLC35A3, and PALMD in the 3_43 region, and RGS22 and VPS13B in 14_66 region, were considered the most interesting candidate genes for buffalo milk production traits.
      Table 4Details of candidate genes within the 3_43 and 14_66 genomic regions
      Region
      Regions: 3_43, the region on BTA3 spanning 43.3 to 43.8 Mb; 14_66, the region on BTA14 spanning 66.02 to 66.15 Mb.
      Within gene
      Within genes, genes within candidate region are based on Bos taurus UMD 3.1 genome assembly (http://bovinegenome.org/?q=node/61); those in bold are the most interesting candidate genes based on their effect on milk production performance in dairy cattle or identified as harboring significant SNP associated with milk production traits of buffalo.
      Gene start (bp)Gene end (bp)Description
      3_43RTCA43,157,18843,185,261RNA 3′-terminal phosphate cyclase
      DBT43,189,90543,230,271Dihydrolipoamide branched chain transacylase E2
      LRRC3943,251,25143,272,503Leucine rich repeat containing 39
      TRMT1343,270,87443,291,189tRNA methyltransferase 13 homolog
      SASS643,291,33343,331,499Spindle assembly 6 homolog
      MFSD14A43,332,80143,386,388Major facilitator superfamily domain containing 14A
      SLC35A343,400,34643,444,844Solute carrier family 35 member A3
      AGL43,504,60143,585,149Amylo-α-1, 6-glucosidase, 4-α-glucanotransferase
      FRRS143,619,57043,629,406Ferric chelate reductase 1
      PALMD43,688,10343,748,131Palmdelphin
      14_66RGS2266,472,37966,593,433Regulator of G-protein signaling 22
      COX6C66,637,80166,647,721Cytochrome c oxidase subunit 6C
      VPS13B66,648,39567,461,111Vacuolar protein sorting 13 homolog B
      SNORA7066,777,37866,777,512Small nucleolar RNA
      1 Regions: 3_43, the region on BTA3 spanning 43.3 to 43.8 Mb; 14_66, the region on BTA14 spanning 66.02 to 66.15 Mb.
      2 Within genes, genes within candidate region are based on Bos taurus UMD 3.1 genome assembly (http://bovinegenome.org/?q=node/61); those in bold are the most interesting candidate genes based on their effect on milk production performance in dairy cattle or identified as harboring significant SNP associated with milk production traits of buffalo.
      On the 3_44 region, MFSD14A was mapped to 43,332,801–43,386,388 bp on BTA3; MFSD14A is a protein-coding gene with product homology to the solute carrier protein family (SLC;
      • Sreedharan S.
      • Stephansson O.
      • Schiöth H.B.
      • Fredriksson R.
      Long evolutionary conservation and considerable tissue specificity of several atypical solute carrier transporters.
      ;
      • Doran J.
      • Walters C.
      • Kyle V.
      • Wooding P.
      • Hammett-Burke R.
      • Colledge W.H.
      Mfsd14a (Hiat1) gene disruption causes globozoospermia and infertility in male mice.
      ). Gene ontology analysis showed that MFSD14A plays roles in the biological process of transmembrane transport (GO:0055085), molecular function of transporter activity (GO:0005215) and cellular components including GO:0016020 and GO:0016021.
      • Whitworth K.M.
      • Zhao J.
      • Spate L.D.
      • Li R.
      • Prather R.S.
      Scriptaid corrects gene expression of a few aberrantly reprogrammed transcripts in nuclear transfer pig blastocyst stage embryos.
      found that MFSD14A downregulated transcripts of the nuclear transfer in pig blastocyst stage, and
      • Doran J.
      • Walters C.
      • Kyle V.
      • Wooding P.
      • Hammett-Burke R.
      • Colledge W.H.
      Mfsd14a (Hiat1) gene disruption causes globozoospermia and infertility in male mice.
      reported that the gene might be involved in mice spermatogenesis. Although the functional study of MFSD14A is limited, we identified 2 significant SNP (AX-85106096 and AX-85148558) located in the gene, which implied MFSD14A may be a novel candidate gene and had genetic effects on buffalo milk fat yield.
      The SLC35A3 gene has been identified as the causative gene of vertebral malformation (
      • Thomsen B.
      • Horn P.
      • Panitz F.
      • Bendixen E.
      • Petersen A.H.
      • Holm L.-E.
      • Nielsen V.H.
      • Agerholm J.S.
      • Arnbjerg J.
      • Bendixen C.
      A missense mutation in the bovine SLC35A3 gene, encoding a UDP-N-acetylglucosamine transporter, causes complex vertebral malformation.
      ), which was reported to negatively affect milk production in dairy cattle (
      • Kadri N.K.
      • Sahana G.
      • Charlier C.
      • Iso-Touru T.
      • Guldbrandtsen B.
      • Karim L.
      • Nielsen U.S.
      • Panitz F.
      • Aamand G.P.
      • Schulman N.
      A 660-Kb deletion with antagonistic effects on fertility and milk production segregates at high frequency in Nordic Red cattle: Additional evidence for the common occurrence of balancing selection in livestock.
      ).
      • Chu Q.
      • Zhang Y.
      • Sun D.
      • Yu Y.
      • Wang Y.
      Identification of the complex vertebral malformation gene in Chinese Holstein and its association with dairy performance traits.
      studied association analyses between dairy performance traits and polymorphism of SLC35A3 and provided evidence that SLC35A3 could influence milk performance. The SLC35A3 gene encodes a protein that participates in molecular pathways involved mainly in transport of nucleotide sugars (R-HSA-727802) and transport of vitamins, nucleosides, and related molecules (R-BTA-425397). The GO annotations related to this gene including sugar:proton symporter activity (GO:0005351) and UDP-N-acetylglucosamine transmembrane transporter activity (GO:0015788). Hence, SLC35A3 should be an interesting candidate gene for milk performance.
      The palmdelphin (PALMD) gene, also called HIAT1, codes a cytosolic protein implicated in p53 phosphorylation (
      • Moioli B.
      • Scatà M.
      • Catillo G.
      • Napolitano F.
      • Pilla F.
      Whole-genome genotyping for detecting mutations in the genes that are directly responsible for the trait variation.
      ). The PALMD gene is thought to be involved in the variation of milk yield. It is differentially expressed between high-yielding and low-yielding Holsten cattle (
      • Seo M.
      • Lee H.-J.
      • Kim K.
      • Caetano-Anolles K.
      • Jeong J.Y.
      • Park S.
      • Oh Y.K.
      • Cho S.
      • Kim H.
      Characterizing milk production related genes in Holstein using RNA-seq.
      ) and
      • Moioli B.
      • Scatà M.C.
      • Steri R.
      • Napolitano F.
      • Catillo G.
      Signatures of selection identify loci associated with milk yield in sheep.
      showed that a mutation in exon 5 (c.547 T>A) of PALMD caused a reduction in the milk-producing ability of sheep. In a study to identify potential causal mutations related to pig production traits, PALMD polymorphisms had significant effects on back fat thickness (
      • Martínez-Montes A.M.
      • Fernandez A.
      • Pérez-Montarelo D.
      • Alves E.
      • Benitez R.
      • Nunez Y.
      • Ovilo C.
      • Ibañez-Escriche N.
      • Folch J.
      • Fernandez A.
      Using RNA-Seq SNP data to reveal potential causal mutations related to pig production traits and RNA editing.
      ). These functional studies of PALMD indicate it could be an important target gene with an important role in the regulation of milk production.
      On the 14_66 region, RGS22 was identified as one of the most interesting genes because it was previously reported as a candidate gene influencing protein yield in Holstein cattle (
      • Marques E.
      • Grant J.
      • Wang Z.
      • Kolbehdari D.
      • Stothard P.
      • Plastow G.
      • Moore S.
      Identification of candidate markers on bovine chromosome 14 (BTA14) under milk production trait quantitative trait loci in Holstein.
      ). From the GeneCards database, we discovered that RGS22 was more overexpressed in breast (42.6) than in other tissues. Similarly, some functional genes controlling of milk properties, likes DGAT1 (
      • Cases S.
      • Smith S.J.
      • Zheng Y.-W.
      • Myers H.M.
      • Lear S.R.
      • Sande E.
      • Novak S.
      • Collins C.
      • Welch C.B.
      • Lusis A.J.
      Identification of a gene encoding an acyl CoA:diacylglycerol acyltransferase, a key enzyme in triacylglycerol synthesis.
      ), ABCG2 (
      • Bionaz M.
      • Loor J.J.
      Gene networks driving bovine milk fat synthesis during the lactation cycle.
      ), and SCD1 (
      • Lengi A.J.
      • Corl B.A.
      Identification and characterization of a novel bovine stearoyl-CoA desaturase isoform with homology to human SCD5.
      ), also show high expression patterns in lactating mammary tissue. Among its related pathways are signaling by GPCR (R-HSA-372790) and peptide ligand-binding receptors (R-HSA-418594). The GO annotations related to this gene include GTPase activator activity (GO:0005096).
      The VPS13B gene on BTA14 was located at 66,648,395 to 67,461,111 bp. One significant SNP (AX-85063131) and 4 suggestive SNP (AX-85063132, AX-85104008, AX-85104009, and AX-85079883) detected for FY270 were harbored in this gene. The VPS13B gene encodes a potential transmembrane protein, which plays a role in the development and function of the eye, hematological system, and central nervous system. The main biological process described for VPS13B is protein transport (GO:0015031), and it may function in vesicle-mediated transport and protein sorting within cells. In a previous study, VPS13B was detected within a QTL associated with leg morphology in dairy cattle (
      • van den Berg I.
      • Fritz S.
      • Rodriguez S.
      • Rocha D.
      • Boussaha M.
      • Lund M.S.
      • Boichard D.
      Concordance analysis for QTL detection in dairy cattle: A case study of leg morphology.
      ). In the review by
      • Capitan A.
      • Michot P.
      • Baur A.
      • Saintilan R.
      • Hoze C.
      • Valour D.
      • Guillaume F.
      • Boichon D.
      • Barbat A.
      • Boichard D.
      Genetic tools to improve reproduction traits in dairy cattle.
      , it was suggested that the genomic region on BTA14 covering RGS22, COX6C, and VPS13B was associated with female fertility as well as milk production traits in Holstein cattle. Accordingly, all findings indicate that the 14_66 genomic region may be an important QTL, with RGS22 as a key candidate gene related to milk production or milk content in dairy buffalo.

      Genomic Prediction

      For the identified regions and markers affecting buffalo milk performance, it is essential to verify that the association can be generalized to other populations. Hence, we used a 5-fold cross validation to evaluate the reliability of GEBV based on genotype data from validation sets (Table 5). The average reliability of GEBV of 6 buffalo milk production traits ranged from 0.06 (PY270) to 0.22 (PP270). The accuracy of GEBV was also assessed by the strength of correlation between GEBV and EBV (
      • Gondro C.
      • Van der Werf J.
      • Hayes B.
      Genome-wide association studies and genomic prediction.
      ), which varied from 0.23 to 0.35 in the present validation set. The observed reliabilities of GEBV in our population were at a similar level to those of Nordic Red (0.19–0.23) using a single-country reference population (
      • Brøndum R.F.
      • Rius-Vilarrasa E.
      • Strandén I.
      • Su G.
      • Guldbrandtsen B.
      • Fikse W.
      • Lund M.S.
      Reliabilities of genomic prediction using combined reference data of the Nordic Red dairy cattle populations.
      ). The correlations between GEBV and published EBV ranged from 0.25 to 0.70 in Holstein, which were higher than those in our study (
      • Su G.
      • Guldbrandtsen B.
      • Gregersen V.
      • Lund M.
      Preliminary investigation on reliability of genomic estimated breeding values in the Danish Holstein population.
      ). The number of available phenotypes, pedigree information, and genotypes in training sets contribute to the accuracy of genomic prediction (
      • Pryce J.E.
      • Arias J.
      • Bowman P.
      • Davis S.
      • Macdonald K.
      • Waghorn G.
      • Wales W.
      • Williams Y.
      • Spelman R.
      • Hayes B.
      Accuracy of genomic predictions of residual feed intake and 250-day body weight in growing heifers using 625,000 single nucleotide polymorphism markers.
      ). Generally, reliability is improved by including a greater number of SNP and increasing the reference population size as reported in dairy (
      • Goddard M.E.
      • Hayes B.J.
      Mapping genes for complex traits in domestic animals and their use in breeding programmes.
      ) and beef (
      • Brito F.V.
      • Neto J.B.
      • Sargolzaei M.
      • Cobuci J.A.
      • Schenkel F.S.
      Accuracy of genomic selection in simulated populations mimicking the extent of linkage disequilibrium in beef cattle.
      ). The heritability of the trait is another factor that influences genomic prediction accuracy (
      • Hayes B.J.
      • Bowman P.
      • Chamberlain A.
      • Goddard M.
      Invited review: Genomic selection in dairy cattle: Progress and challenges.
      ,
      • Brito F.V.
      • Neto J.B.
      • Sargolzaei M.
      • Cobuci J.A.
      • Schenkel F.S.
      Accuracy of genomic selection in simulated populations mimicking the extent of linkage disequilibrium in beef cattle.
      ). Indeed, traits with higher heritability had higher reliability of GEBV in our prediction. The best accuracy was observed for PP270, which had the highest heritability (0.38) among the studied traits. In addition, the method used to calculate the reliability of GEBV also influenced results. Therefore, a more precise buffalo genomic map, a larger sample size, and higher SNP density are required to improve the detection power and prediction accuracy.
      Table 5The reliability of genomic EBV (GEBV) and correlation between EBV and GEBV for the validation sets (mean ± SEM)
      Trait
      PM = peak milk yield; MY270 = 270-d total milk yield; FY270 = 270-d fat yield; FP270 = 270-d fat percentage; PY270 = 270-d protein yield; PP270 = 270-d protein percentage.
      ReliabilityCorrelation
      PM0.14 ± 0.0030.23 ± 0.05
      MY2700.08 ± 0.0020.26 ± 0.03
      FY2700.12 ± 0.0030.30 ± 0.03
      FP2700.19 ± 0.0040.32 ± 0.04
      PY2700.06 ± 0.0020.27 ± 0.03
      PP2700.22 ± 0.0040.35 ± 0.03
      1 PM = peak milk yield; MY270 = 270-d total milk yield; FY270 = 270-d fat yield; FP270 = 270-d fat percentage; PY270 = 270-d protein yield; PP270 = 270-d protein percentage.

      CONCLUSIONS

      The aim of the present study was to investigate important genomic regions and genes associated with milk production traits in water buffalo. Genomic prediction was conducted to explore the feasibility of genomic selection as a potential strategy in buffalo breeding programs. Four SNP (AX-85148558, AX-85106096, AX-85073877, and AX-85063131) and 2 genomic region (3_43 and 14_66) were found to be associated with buffalo milk performance. Five genes, including MFSD14A, SLC35A3, PALMD, RGS22, and VPS13B, were identified as novel candidate genes for milk production in buffaloes. The reliability of genomic prediction for all traits ranged from 0.06 to 0.22. These findings provide useful preliminary information for buffalo breeding programs to adopt genomic selection. However, a more precise buffalo genomic map is required to identify the genes affecting the traits, and a larger sample size may help to improve the detection power and prediction accuracy. Furthermore, studies are needed to evaluate other interesting milk traits, such as coagulation ability and titratable acidity, which are important technological traits in producing cheese.

      ACKNOWLEDGMENTS

      This project was supported by 948 Key Project of China (No. 2011-G26) and International Cooperation Key Project of China (No. 2011 DFA32250). The authors acknowledge the contributions of the 2 (Italy and Canada) research institutions and all field staff and participants in our research. Authors L. G. Yang, G. Campanile, and G. Plastow conceived and designed the experiments; A. X. Liang, A. Salzano, and B. Gasparrini collected samples and records; C. Zhang and J. J. Liu performed the experiments; Z. Wang, C. Zhang, and J. J. Liu analyzed the data; J. J. Liu, C. Zhang, A. X. Liang, and M. Cassandro wrote the paper; and L. G. Yang, G. Campanile, and G. Plastow critically reviewed the manuscript. The authors declare no competing financial interests.

      REFERENCES

        • Amaral M.E.J.
        • Grant J.R.
        • Riggs P.K.
        • Stafuzza N.B.
        • Edson Filho A.
        • Goldammer T.
        • Weikard R.
        • Brunner R.M.
        • Kochan K.J.
        • Greco A.J.
        A first generation whole genome RH map of the river buffalo with comparison to domestic cattle.
        BMC Genomics. 2008; 9 (19108729): 631
        • Baldi F.
        • Laureano M.M.M.
        • Gordo D.G.M.
        • Bignardi A.B.
        • Borquis R.R.A.
        • Albuquerque L.G.d.
        • Tonhati H.
        Effect of lactation length adjustment procedures on genetic parameter estimates for buffalo milk yield.
        Genet. Mol. Biol. 2011; 34 (21637545): 62-67
        • Barrett J.C.
        • Fry B.
        • Maller J.
        • Daly M.J.
        Haploview: Analysis and visualization of LD and haplotype maps.
        Bioinformatics. 2005; 21 (15297300): 263-265
        • Bionaz M.
        • Loor J.J.
        Gene networks driving bovine milk fat synthesis during the lactation cycle.
        BMC Genomics. 2008; 9 (18671863): 366
        • Brito F.V.
        • Neto J.B.
        • Sargolzaei M.
        • Cobuci J.A.
        • Schenkel F.S.
        Accuracy of genomic selection in simulated populations mimicking the extent of linkage disequilibrium in beef cattle.
        BMC Genet. 2011; 12 (21933416): 80
        • Brøndum R.F.
        • Rius-Vilarrasa E.
        • Strandén I.
        • Su G.
        • Guldbrandtsen B.
        • Fikse W.
        • Lund M.S.
        Reliabilities of genomic prediction using combined reference data of the Nordic Red dairy cattle populations.
        J. Dairy Sci. 2011; 94 (21854944): 4700-4707
        • Capitan A.
        • Michot P.
        • Baur A.
        • Saintilan R.
        • Hoze C.
        • Valour D.
        • Guillaume F.
        • Boichon D.
        • Barbat A.
        • Boichard D.
        Genetic tools to improve reproduction traits in dairy cattle.
        Reprod. Fertil. Dev. 2014; 27 (25472040): 14-21
        • Cases S.
        • Smith S.J.
        • Zheng Y.-W.
        • Myers H.M.
        • Lear S.R.
        • Sande E.
        • Novak S.
        • Collins C.
        • Welch C.B.
        • Lusis A.J.
        Identification of a gene encoding an acyl CoA:diacylglycerol acyltransferase, a key enzyme in triacylglycerol synthesis.
        Proc. Natl. Acad. Sci. USA. 1998; 95 (9789033): 13018-13023
        • Chu Q.
        • Zhang Y.
        • Sun D.
        • Yu Y.
        • Wang Y.
        Identification of the complex vertebral malformation gene in Chinese Holstein and its association with dairy performance traits.
        Hereditas. 2010; 32 (20650855): 732-736
        • de Camargo G.M.
        • Aspilcueta-Borquis R.
        • Fortes M.
        • Porto-Neto R.
        • Cardoso D.
        • Santos D.
        • Lehnert S.
        • Reverter A.
        • Moore S.
        • Tonhati H.
        Prospecting major genes in dairy buffaloes.
        BMC Genomics. 2015; 16 (26510479): 872
        • Di Meo G.P.
        • Perucatti A.
        • Floriot S.
        • Hayes H.
        • Schibler L.
        • Incarnato D.
        • Di Berardino D.
        • Williams J.
        • Cribiu E.
        • Eggen A.
        An extended river buffalo (Bubalus bubalis, 2n= 50) cytogenetic map: assignment of 68 autosomal loci by FISH-mapping and R-banding and comparison with human chromosomes.
        Chromosome Res. 2008; 16 (18685962): 827-837
        • Do D.N.
        • Janss L.L.G.
        • Jensen J.
        • Kadarmideen H.N.
        SNP annotation-based whole genomic prediction and selection: An application to feed efficiency and its component traits in pigs.
        J. Anim. Sci. 2015; 93 (26020301): 2056-2063
        • Doran J.
        • Walters C.
        • Kyle V.
        • Wooding P.
        • Hammett-Burke R.
        • Colledge W.H.
        Mfsd14a (Hiat1) gene disruption causes globozoospermia and infertility in male mice.
        Reproduction. 2016; 152 (27107036): 91-99
        • El-Halawany N.
        • Abdel-Shafy H.
        • Abd-El-Monsif A.S.
        • Abdel-Latif M.A.
        • Al-Tohamy A.F.
        • El-Moneim O.M.A.
        Genome-wide association study for milk production in Egyptian buffalo.
        Livest. Sci. 2017; 198: 10-16
      1. El-Halawany, N., H. Abdel-Shafy, A. E. M. A. Shawky, M. A. Abdel-Latif, A. F. M. Al-Tohamy, and O. M. A. El-Moneim. 2015. Genome-wide association study for milk production in Egyptian buffalo. Pages 10–16 in Proc. Int. Symp. Animal Functional Genomics. The 6th International Symposium on Animal Functional Genomics, Piacenza, Italy.

        • El Nahas S.M.
        • de Hondt H.A.
        • Womack J.E.
        Current status of the river buffalo (Bubalus bubalis L.) gene map.
        J. Hered. 2001; 92 (11447236): 221-225
        • Elsik C.G.
        • Unni D.R.
        • Diesh C.M.
        • Tayal A.
        • Emery M.L.
        • Nguyen H.N.
        • Hagen D.E.
        Bovine Genome Database: New tools for gleaning function from the Bos taurus genome.
        Nucleic Acids Res. 2016; 44 (26481361): D834-D839
        • Endelman J.B.
        Ridge regression and other kernels for genomic selection with R package rrBLUP.
        Plant Genome. 2011; 4: 250-255
        • FAOSTAT
        Food and Agriculture Organization the United Nations Statistics Division.
        • Fulton J.E.
        Genomic selection for poultry breeding.
        Anim. Front. 2012; 2: 30-36
        • Garrick D.J.
        • Taylor J.F.
        • Fernando R.L.
        Deregressing estimated breeding values and weighting information for genomic regression analyses.
        Genet. Sel. Evol. 2009; 41 (20043827): 55
        • Genovese C.R.
        • Lazar N.A.
        • Nichols T.
        Thresholding of statistical maps in functional neuroimaging using the false discovery rate.
        Neuroimage. 2002; 15 (11906227): 870-878
        • Gilmour A.R.
        • Gogel B.
        • Cullis B.
        • Thompson R.
        • Butler D.
        ASReml user guide release 3.0.
        VSN International Ltd., Hemel Hempstead, UK2009
        • Goddard M.E.
        • Hayes B.J.
        Mapping genes for complex traits in domestic animals and their use in breeding programmes.
        Nat. Rev. Genet. 2009; 10 (19448663): 381-391
        • Gondro C.
        • Van der Werf J.
        • Hayes B.
        Genome-wide association studies and genomic prediction.
        Springer, Amsterdam, the Netherlands2013
        • Harder B.
        • Bennewitz J.
        • Reinsch N.
        • Thaller G.
        • Thomsen H.
        • Kühn C.
        • Schwerin M.
        • Erhardt G.
        • Förster M.
        • Reinhardt F.
        Mapping of quantitative trait loci for lactation persistency traits in German Holstein dairy cattle.
        J. Anim. Breed. Genet. 2006; 123 (16533362): 89-96
        • Hayes B.J.
        • Bowman P.
        • Chamberlain A.
        • Goddard M.
        Invited review: Genomic selection in dairy cattle: Progress and challenges.
        J. Dairy Sci. 2009; 92 (19164653): 433-443
        • Hu Z.-L.
        • Park C.A.
        • Reecy J.M.
        Developmental progress and current status of the Animal QTLdb.
        Nucleic Acids Res. 2016; 44: D827-D833
        • Iamartino D.
        • Williams J.L.
        • Sonstegard T.
        • Reecy J.
        • v. Tassell C.
        • Nicolazzi E.L.
        • Biffani S.
        • Biscarini F.
        • Schroeder S.
        • de Oliveira D.A.
        The buffalo genome and the application of genomics in animal management and improvement.
        in: Proc. Buffalo Bulletin. International Buffalo Information Center, Bangkok, Thailand2013: 151-158
        • Kadri N.K.
        • Sahana G.
        • Charlier C.
        • Iso-Touru T.
        • Guldbrandtsen B.
        • Karim L.
        • Nielsen U.S.
        • Panitz F.
        • Aamand G.P.
        • Schulman N.
        A 660-Kb deletion with antagonistic effects on fertility and milk production segregates at high frequency in Nordic Red cattle: Additional evidence for the common occurrence of balancing selection in livestock.
        PLoS Genet. 2014; 10 (24391517): e1004049
        • Khan M.
        • Hyder A.
        • Bajwa I.
        • Rehman M.
        • Hassan F.
        Prediction of lactation yield from last-record day and average daily yield in Nili-Ravi buffaloes.
        Pak. Vet. J. 2005; 25: 175
        • Khedkar C.
        • Kalyankar S.
        • Deosarkar S.
        Buffalo milk.
        in: Caballero B. Finglas P.M. Toldrá F. Encyclopedia of Food and Health. Academic Press, New York, NY2016: 522-528
        • Lengi A.J.
        • Corl B.A.
        Identification and characterization of a novel bovine stearoyl-CoA desaturase isoform with homology to human SCD5.
        Lipids. 2007; 42 (17468887): 499-508
        • Lenth R.V.
        Least squares means: The R package lsmeans.
        J. Stat. Softw. 2016; 69: 1-33
        • Lukić B.
        • Pong-Wong R.
        • Rowe S.
        • Koning D.
        • Velander I.
        • Haley C.
        • Archibald A.
        • Woolliams J.
        Efficiency of genomic prediction for boar taint reduction in Danish Landrace pigs.
        Anim. Genet. 2015; 46 (26449733): 607-616
        • Marques E.
        • Grant J.
        • Wang Z.
        • Kolbehdari D.
        • Stothard P.
        • Plastow G.
        • Moore S.
        Identification of candidate markers on bovine chromosome 14 (BTA14) under milk production trait quantitative trait loci in Holstein.
        J. Anim. Breed. Genet. 2011; 128 (21749477): 305-313
        • Martínez-Montes A.M.
        • Fernandez A.
        • Pérez-Montarelo D.
        • Alves E.
        • Benitez R.
        • Nunez Y.
        • Ovilo C.
        • Ibañez-Escriche N.
        • Folch J.
        • Fernandez A.
        Using RNA-Seq SNP data to reveal potential causal mutations related to pig production traits and RNA editing.
        Anim. Genet. 2017; 48 (27642173): 151-165
        • Miar Y.
        • Plastow G.
        • Moore S.
        • Manafiazar G.
        • Charagu P.
        • Kemp R.
        • Van Haandel B.
        • Huisman A.
        • Zhang C.
        • McKay R.
        Genetic and phenotypic parameters for carcass and meat quality traits in commercial crossbred pigs.
        J. Anim. Sci. 2014; 92 (24778330): 2869-2884
        • Michelizzi V.N.
        • Dodson M.V.
        • Pan Z.
        • Amaral M.E.J.
        • Michal J.J.
        • McLean D.J.
        • Womack J.E.
        • Jiang Z.
        Water buffalo genome science comes of age.
        Int. J. Biol. Sci. 2010; 6 (20582226): 333-349
        • Moioli B.
        • Scatà M.
        • Catillo G.
        • Napolitano F.
        • Pilla F.
        Whole-genome genotyping for detecting mutations in the genes that are directly responsible for the trait variation.
        in: Proc. Plant and Animal Genome XXI, San Diego, CA. 2013: P0621
        • Moioli B.
        • Scatà M.C.
        • Steri R.
        • Napolitano F.
        • Catillo G.
        Signatures of selection identify loci associated with milk yield in sheep.
        BMC Genet. 2013; 14 (24004915): 76
        • Moser G.
        • Tier B.
        • Crump R.E.
        • Khatkar M.S.
        • Raadsma H.W.
        A comparison of five methods to predict genomic breeding values of dairy bulls from genome-wide SNP markers.
        Genet. Sel. Evol. 2009; 41 (20043835): 56
        • Ostersen T.
        • Christensen O.F.
        • Henryon M.
        • Nielsen B.
        • Su G.
        • Madsen P.
        Deregressed EBV as the response variable yield more reliable genomic predictions than traditional EBV in pure-bred pigs.
        Genet. Sel. Evol. 2011; 43 (22070746): 38
        • Pryce J.E.
        • Arias J.
        • Bowman P.
        • Davis S.
        • Macdonald K.
        • Waghorn G.
        • Wales W.
        • Williams Y.
        • Spelman R.
        • Hayes B.
        Accuracy of genomic predictions of residual feed intake and 250-day body weight in growing heifers using 625,000 single nucleotide polymorphism markers.
        J. Dairy Sci. 2012; 95 (22459856): 2108-2119
        • Purcell S.
        • Chang C.
        PLINK 1.9.
        • Rosati A.
        • Van Vleck L.D.
        Estimation of genetic parameters for milk, fat, protein and mozzarella cheese production for the Italian river buffalo Bubalus bubalis population.
        Livest. Prod. Sci. 2002; 74: 185-190
        • Saatchi M.
        • Schnabel R.D.
        • Taylor J.F.
        • Garrick D.J.
        Large-effect pleiotropic or closely linked QTL segregate within and across ten US cattle breeds.
        BMC Genomics. 2014; 15 (24906442): 442
        • Sargolzaei M.
        • Schenkel F.S.
        • Vanraden P.M.
        gebv: Genomic breeding value estimator for livestock. Technical report to the Dairy Cattle Breeding and Genetics Committee.
        University of Guelph, 2009
        • Seo M.
        • Lee H.-J.
        • Kim K.
        • Caetano-Anolles K.
        • Jeong J.Y.
        • Park S.
        • Oh Y.K.
        • Cho S.
        • Kim H.
        Characterizing milk production related genes in Holstein using RNA-seq.
        Asian-australas. J. Anim. Sci. 2016; 29 (26950864): 343
        • Shi D.S.
        • Wang J.
        • Yang Y.
        • Lu F.H.
        • Li X.P.
        • Liu Q.Y.
        DGAT1, GH, GHR, PRL and PRLR polymorphism in water buffalo (Bubalus bubalis).
        Reprod. Domest. Anim. 2012; 47 (21883511): 328-334
        • Sreedharan S.
        • Stephansson O.
        • Schiöth H.B.
        • Fredriksson R.
        Long evolutionary conservation and considerable tissue specificity of several atypical solute carrier transporters.
        Gene. 2011; 478 (21044875): 11-18
        • Su G.
        • Guldbrandtsen B.
        • Gregersen V.
        • Lund M.
        Preliminary investigation on reliability of genomic estimated breeding values in the Danish Holstein population.
        J. Dairy Sci. 2010; 93 (20172238): 1175-1183
        • Tantia M.S.
        • Vijh R.K.
        • Mishra B.P.
        • Mishra B.
        • Kumar S.B.
        • Sodhi M.
        DGAT1 and ABCG2 polymorphism in Indian cattle (Bos indicus) and buffalo (Bubalus bubalis) breeds.
        BMC Vet. Res. 2006; 2 (17087837): 32
        • Thomsen B.
        • Horn P.
        • Panitz F.
        • Bendixen E.
        • Petersen A.H.
        • Holm L.-E.
        • Nielsen V.H.
        • Agerholm J.S.
        • Arnbjerg J.
        • Bendixen C.
        A missense mutation in the bovine SLC35A3 gene, encoding a UDP-N-acetylglucosamine transporter, causes complex vertebral malformation.
        Genome Res. 2006; 16 (16344554): 97-105
        • Tonhati H.
        • Cerón-Muñoz M.F.
        • Oliveira J.A.d.
        • El Faro L.
        • Lima A.L.F.
        • Albuquerque L.G.d.
        Test-day milk yield as a selection criterion for dairy buffaloes (Bubalus bubalis Artiodactyla, Bovidae).
        Genet. Mol. Biol. 2008; 31: 674-679
        • Tonhati H.
        • Muñoz M.
        • Duarte J.
        • Reichert R.
        • Oliveira J.
        • Lima A.
        Estimates of correction factors for lactation length and genetic parameters for milk yield in buffaloes.
        Arq. Bras. Med. Vet. Zootec. 2004; 56: 251-257
        • van den Berg I.
        • Fritz S.
        • Rodriguez S.
        • Rocha D.
        • Boussaha M.
        • Lund M.S.
        • Boichard D.
        Concordance analysis for QTL detection in dairy cattle: A case study of leg morphology.
        Genet. Sel. Evol. 2014; 46 (24884971): 31
        • Venturini G.C.
        • Cardoso D.
        • Baldi F.
        • Freitas A.
        • Aspilcueta-Borquis R.
        • Santos D.
        • Tonhati H.
        Association between single-nucleotide polymorphisms and milk production traits in buffalo.
        Genet. Mol. Res. 2014; 13 (25501237): 10256
        • Warriach H.M.
        • Mcgill D.M.
        • Bush R.D.
        • Wynn P.C.
        • Chohan K.R.
        A review of recent developments in buffalo reproduction - a review.
        Asian-australas. J. Anim. Sci. 2015; 28 (25656203): 451-455
        • Weller J.I.
        • Ron M.
        Invited review: quantitative trait nucleotide determination in the era of genomic selection.
        J. Dairy Sci. 2011; 94 (21338774): 1082-1090
        • Whitworth K.M.
        • Zhao J.
        • Spate L.D.
        • Li R.
        • Prather R.S.
        Scriptaid corrects gene expression of a few aberrantly reprogrammed transcripts in nuclear transfer pig blastocyst stage embryos.
        Cell. Reprogram. 2011; 13 (21548827): 191-204
        • Wibowo T.A.
        • Gaskins C.T.
        • Newberry R.C.
        • Thorgaard G.H.
        • Michal J.J.
        • Jiang Z.
        Genome assembly anchored QTL map of bovine chromosome 14.
        Int. J. Biol. Sci. 2008; 4 (19043607): 406
        • Wu J.J.
        • Song L.J.
        • Wu F.J.
        • Liang X.W.
        • Yang B.Z.
        • Wathes D.C.
        • Pollott G.E.
        • Cheng Z.
        • Shi D.S.
        • Liu Q.Y.
        Investigation of transferability of BovineSNP50 BeadChip from cattle to water buffalo for genome wide association study.
        Mol. Biol. Rep. 2013; 40 (23232712): 743-750
        • Yang S.-H.
        • Bi X.-J.
        • Xie Y.
        • Li C.
        • Zhang S.-L.
        • Zhang Q.
        • Sun D.-X.
        Validation of PDE9A gene identified in GWAS showing strong association with milk production traits in Chinese Holstein.
        Int. J. Mol. Sci. 2015; 16 (26556348): 26530-26542