Genome-Wide Association Studies in Indian Buffalo Revealed Genomic Regions for Lactation and Fertility

Abstract

Murrah breed of buffalo is an excellent dairy germplasm known for its superior milk quality in terms of milk fat and solids-not-fat (SNF); however, it is often reported that Indian buffaloes had lower lactation and fertility potential compared to the non-native cattle of the country. Recent techniques, particularly the genome-wide association studies (GWAS), to identify genomic variations associated with lactation and fertility traits offer prospects for systematic improvement of buffalo. DNA samples were sequenced using the double-digestion restriction-associated DNA (RAD) tag genotyping-by-sequencing. The bioinformatics pipeline was standardized to call the variants, and single-nucleotide polymorphisms (SNPs) qualifying the stringent quality check measures were retained for GWAS. Over 38,000 SNPs were used to perform GWAS on the first two principal components of test-day records of milk yields, fat percentages, and SNF percentages, separately. GWAS was also performed on 305 days’ milk yield; lactation persistency was estimated through the rate of decline after attaining the peak yield method, along with three other standard methods; and breeding efficiency, post-partum breeding interval, and age at sexual maturity were considered fertility traits. Significant association of SNPs was observed for the first principal component, explaining the maximum proportion of variation in milk yield. Furthermore, some potential genomic regions were identified to have a potential role in regulating milk yield and fertility in Murrah. Identification of such genomic regions shall help in carrying out an early selection of high-yielding persistent Murrah buffaloes and, in the long run, would be helpful in shaping their future genetic improvement programs.

Introduction

Buffalo (Bubalus bubalis) is an imperative livestock species and act as a key component for improving agricultural economy and supplying milk, meat, and draft power. The buffalo population across the world was recently estimated to be 194 million, 97% of which were present in Asia (FAO, 2014). 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 (Venturini et al., 2014; Khedkar et al., 2016). The concentration of these milk constituents offers a higher economic return of the buffalo milk and increases the demand of value-added products like mozzarella. Buffaloes are an integral and crucial genetic resource in Indian dairy industry, contributing 49% to the total milk produced according to the basic animal husbandry statistics, 2019¹. About 63% of global buffalo milk production and 95% of Asian buffalo milk production is contributed by Indian buffaloes (FAO, 2014). Buffaloes are more adaptable to harsh environments and often resist various bovine tropical diseases. However, the poor reproductive efficiency of buffaloes limits its potential. Buffaloes exhibit higher age at puberty and maturity, longer postpartum breeding interval, and low conception rates (Gordon, 1996; Warriach et al., 2015; Raina et al., 2016). Buffaloes in the field condition also suffer from short lactation of 252–270 days as compared to the standard 300 days (Ganguli, 1981; Griffiths, 2010). This underlines the scope for improving the milk production and fertility potential through implementation of systematic breeding programs for buffaloes in the country. However, in India, very few works have been done in order to identify the infinitesimally large number of underlying loci regulating the expression of complex traits such as milk production, lactation persistency and fertility.

Genome wide association studies (GWAS) in buffaloes for lactation traits have been mainly limited to the use of bovine single-nucleotide polymorphism (SNP) chip (Wu et al., 2013; Venturini et al., 2014) and Affymetrix’s buffalo 90K SNP chip (de Camargo et al., 2015; Iamartino et al., 2017; Liu et al., 2018). There have been very few GWAS conducted in India covering all aspects of production and reproduction performance due to constraints of cost incurred and organized large-scale genotyping programs. High-density SNP panels are a prerequisite for GWAS, which have led to developments of cost-effective and efficient next-generation sequencing (NGS) technologies such as reduced representation of genomic libraries (RRLs) (Van Tassell et al., 2008). The flexibility, robustness, and low cost of double-digestion restriction-associated DNA (RAD) tag genotyping-by-sequencing (ddRAD-GBS) technique renders it suitable for identification of SNPs in any species for GWAS (Davey et al., 2011; Elshire et al., 2011).

Hence, the present study was conducted to identify novel SNPs associated with milk production, composition, lactation persistency, and fertility traits at the genomic level using the genotype-by-sequencing technique in Murrah buffalo, India’s major buffalo breed and milch animal of the nation.

Materials and Methods

Sampling, Data Recording, and Genotyping

A total of 672 test-day records on each trait, i.e., milk yield (TDMY), fat percentage (TDFP), and solids-not-fat (SNF) percentage (TDSNF) were collected from 96 female Murrah buffaloes reared at LRC, NDRI, Karnal, India (29.68°N and 76.99°E). Records of 96 buffaloes on 305 days’ milk yield (305DMY), birth weight (bwt), age at first calving (AFC), calving interval (CI), and age at sexual maturity (ASM) in months were collected. Other traits such as lactation persistency, postpartum breeding interval in days (PPBI), and breeding efficiency (BE) were derived from the primary phenotype records. Breeding efficiency was calculated for female buffaloes as described by Tomar (1965).

Lactation persistency was estimated by following four different methods:

• (I)

  Wood (1967) incomplete gamma function, where individuals were classified as persistent and non-persistent as described by Macciotta et al. (2005) considering a positive rate of incline as a favorable condition and a negative rate as an unfavorable condition for persistency

• (II)

  Johansson and Hansson (1940) method

• (III)

  Ludwick and Petersen (1943) method

• (IV)

  Mahadevan (1951) method

DNA was isolated from 96 Murrah buffaloes selected for the study following the standard phenol-chloroform method (Sambrook and Russell, 2006). Samples were further processed using the standard RAD protocol as described by Peterson et al. (2012). DNA double digestion was carried out with SphI and MluCI restriction enzymes. Adapters (P1 and P2) were prepared as per standard Illumina read multiplexing protocol using an inline barcode along with Illumina index for library preparation. After adapter ligation and size selection, samples were sequenced on a Illumina Hi-Seq 2000 platform.

Variant Calling

The NGS pipeline was standardized after incorporating a few modifications in the standard mpileup variant calling pipeline (Li, 2011), to call variants present in the Murrah population. The Mediterranean buffalo genome, having accession ID GCF_003121395.1, was retrieved from the NCBI dataset and used as a reference genome. Index and sequence dictionary files were created using the Burrows–Wheeler algorithm (BWA) (Li and Durbin, 2010) and PicardTools,² respectively. The quality of paired-end raw FASTQ files generated after sequencing was checked using FastQC (Andrews, 2010), and each report was combined through MultiQC (Ewels et al., 2016). Adapters were marked and trimmed using bbmap (Brian, 2014).³ The BWA-MEM algorithm was used to align the trimmed FASTQ sequences with the reference genome. Aligned files were coordinate-sorted, and duplicate reads were removed. Read group identifiers were updated using PicardTools. The quality of aligned BAM files was checked using qualimap (García-Alcalde et al., 2012). Variants were called using bcftools-mpileup (Li, 2011).

Quality Control of Variants

Only biallelic SNPs having more than 95% genotyping rate were retained for further GWAS. SNPs with a MAF < 0.05 and deviating from the Hardy–Weinberg equilibrium at p < 0.0001 were removed from the dataset. SNPs in LD with r² > 0.8 were also removed. Only autosomal and X chromosome SNPs were retained for the final analysis. All the quality control operations and data preprocessing were performed using PLINK v1.9 (Chang et al., 2015).

Statistical Model

GWAS for milk yield, fat percentage, and SNF percentage were conducted on the principal components (PCs) instead of direct traits. Principal component analysis (PCA) was performed in the R programming environment (v4.0.3) on the 672 test-day records of each trait separately. As the first two PCs cumulatively account for most of the variation in the dataset, they were selected to be GWAS traits. The traits on which GWAS were performed were PCs (PC₁ and PC₂) of TDMYs, TDFPs, TDSNFs, and 305DMY; lactation persistency calculated by four different methods; age at sexual maturity in months; postpartum breeding interval (in days); and breeding efficiency. Genome-wide identity-by-state (IBS) for all pairs of individuals was checked. Multidimensional scaling (MDS) based on SNP information was done to check for the presence of any population stratification (Supplementary Figure 1) and was corrected by incorporating the first two MDS components as covariates in the model for GWAS. A genome-wide scan for significant SNPs considering only additive effects was accomplished through a simple regression model using PLINK v1.9 as described by Marees et al. (2018), where residuals were assumed to be normally and independently distributed. A linear regression model was fitted for determining the association between SNPs and continuous traits (Bush and Moore, 2012), while logistic regression was fitted for the binary trait (lactation persistent/non-persistent based on incomplete gamma function). The threshold for genome-wide significance was determined by correcting the p-values of the SNP association test with Benjamini–Hochberg’s false discovery rate (FDR) at 5 and 10% levels (Benjamini and Hochberg, 1995) using the “R” package fuzzySim v3.0 (Barbosa, 2020). A genome-wide significant threshold was set at FDR 5% and a suggestive threshold at 10% by calculating a nominal p-value for the largest index i for which P_(i) ≤ (i/m) × q, where i = rank of the SNPs, m = no. of individual tests performed, and q = either 0.05 or 0.1 (Glickman et al., 2014). The results were plotted as Manhattan plots and Q-Q plots using the “qqman” package of R.

Linear regression model used for GWAS:

where y = trait, x = additive effect of SNPs, C₁ = first component of MDS, C₂ = second component of MDS, β₀ = intercept term, β₁ = regression coefficient representing the strength of association between SNP x and trait y, β₂ = regression coefficient of C₁, β₃ = regression coefficient of C₂, and e = residuals or noise not explained by SNPs.

The model used for GWAS on 305DMY, however, included four covariates consisting of the first component of MDS, birth weight, age at first calving, and calving interval.

SNP Mapping and Pathway Enrichment

SNPs were identified as genic if present within the genes or intergenic if present within a range of 20 kb from the 5′ and 3′ ends of the gene (da Costa Barros et al., 2018). ARS-UCD1.2/bosTau9 cow assembly was used as the reference genome to identify regions around significant SNPs in the UCSC genome browser (Zimin et al., 2009). Gene ontology (GO) was carried out using gProfiler,⁴ and the GO classifications significant over the Benjamini–Hochberg FDR were selected for pathway enrichment through Cytoscape v3.8.2 (Shannon et al., 2003).

Results

PCA on Test-Day Records

PCA on test-day records proved that the first two PCs cumulatively explain 77.72, 40.77, and 48.06% of the total variation in TDMYs, TDFPs, and TDSNFs, respectively. Eigen values of PC₁ for TDMYs, TDFPs, and TDSNFs were found to be 4.45, 1.67, and 2.30, respectively, while for PC₂ they were 0.98, 1.20, and 1.05, respectively. New synthetic variable PC₁ and PC₂ for the above three traits were constructed for each buffalo utilizing original variables and variable loadings of PCA.

Genotyping

An average of 1.3 million each of forward and reverse reads were obtained per sample. The total number of forward and reverse end reads was 252.56 million with the average read length being 151 base pairs (bp). After quality check of raw data, it was observed that the average GC content was 50.54%. The rate of duplication was quite high, i.e., 80.11%, which was later checked during variant calling. The quality score (Q) throughout the dataset varied from 35 to 40.

Variant Calling and Quality Control

Raw files that were quality checked in FastQC were combined using MultiQC, and the report is provided as Supplementary Document. After alignment, BAM files were obtained, and quality was checked. Of the raw reads, 98.11% were mapped accurately to the Mediterranean buffalo reference genome. After removing the duplicate reads, the rate of duplication was reduced to 15.85%. Average mapping quality after alignment was found to be 29.05. A single VCF file was obtained, with 3,854,990 SNPs, out of which 3,792,469 SNPs which were strictly biallelic were retained for further analysis. After applying quality control constraints, 38,560 SNPs present on autosomes and X-chromosomes were retained for further downstream analysis with the final genotyping rate of 98.24%. All the SNPs present on autosome no. 24 failed to pass quality control constraints; hence, in the final analysis, the 24th autosome has been removed.

GWAS Results

GWAS were performed with 38,560 SNPs, and the genome-wide significant threshold was set at the 5% FDR level. GWAS were performed on the two foremost variation explaining PCs (PC₁ and PC₂) of TDMY, TDFP, and TDSNF.

GWAS on Lactation Traits

Three SNPs present on chromosomes X (7588358 bp), 9 (62699319 bp), and 16 (74362861 bp) were significantly associated with PC₁ of TDMYs with FDR-corrected p-values of 0.007369, 0.04579, and 0.04579, respectively. No SNPs were found to be significantly associated with PC₂ TDMYs. GWAS on PC₁ and PC₂ of TDFPs and TDSNFs also could not establish any significant association between the traits and SNPs. However, the top 10 SNPs having the lowest p-values in the test of association with different lactation traits, along with their position and genes within a ± 20-kb region, are listed in Table 1. The Manhattan and Q-Q plots (panels A and B, respectively) for association results of PC₁–TDMYs, PC₂–TDMYs, PC₁–TDFPs, PC₂–TDFPs, PC₁–TDSNFs, and PC₂–TDSNFs are given in Figures 1–6, respectively.

FIGURE 1

FIGURE 2

FIGURE 3

FIGURE 4

FIGURE 5

FIGURE 6

GWAS on 305DMY revealed that there was no significant association between SNPs and the trait. It was observed that five SNPs that appeared in the list with GWAS for PC₁ and PC₂ of the test day’s milk yield were also in the top SNP list for 305DMY. The Manhattan and Q-Q plots for association results of 305DMY are given in Figures 7A,B, respectively.

FIGURE 7

GWAS on Lactation Persistency

GWAS on lactation persistency estimated by four different methods revealed no significant association between SNPs and the trait generated. The top 10 SNPs having the lowest p-values along with their genomic position and genes within a ± 20-kb region are listed in Table 1 for all four methods. It was observed that four SNPs were found in common between the top SNPs listed for persistency estimated by methods II and III. Among the four SNPs, one was on chromosome 9 at 62472303 bp, and three were present on chromosome 18 at 46100792, 46100619, and 17868492 bp. The Manhattan and Q-Q plots (panels A and B, respectively) for association results of persistency estimated by four different methods (I, II, III, and IV) are given in Figures 8–11, respectively.

FIGURE 8

Macciotta et al. (2005). (A) Manhattan plot of genome-wide SNPs. The red line indicates the p-value threshold (expressed as –log₁₀P) corresponding to FDR-corrected p-values or q = 0.05, above which the SNPs are considered to be significantly associated with the trait, while the blue line indicates the genome-wide suggestive threshold at q = 0.1. The SNPs identified on chromosome no. 24 were screened out from analysis because of stringent quality control; hence, it is not represented in the Manhattan plot. (B) Q-Q plot of the p-values.

FIGURE 9

FIGURE 10

FIGURE 11

GWAS on Fertility Traits

Upon performing GWAS on age at sexual maturity, postpartum breeding interval, and breeding efficiency, no significant association of any SNPs were observed for the traits; however, the top 10 SNPs with the lowest p-values and genes within a ± 20-kb region are listed in Table 1. The Manhattan and Q-Q plots (panels A and B, respectively) of GWAS on age at sexual maturity, postpartum breeding interval, and breeding efficiency are given in Figures 12–14, respectively.

FIGURE 12

FIGURE 13

FIGURE 14

Discussion

The test-day model (TDM) as repeated measurements of milk yield traits is the method of choice for predicting 305 days’ milk yield. TDMs account for environmental effects on each test day and is useful to model individual lactation curves (Schaeffer et al., 2000; Bignardi et al., 2012). PCA being a powerful multivariate technique is often used to predict 305 days’ milk yield from the inter-related variables such as test-day records (Taggar et al., 2012). Wara et al. (2019) had performed GWAS on test day’s milk yield in Vrindavani cattle; however, reports of GWAS on PCs of test days are very few. In the present study, PCA was applied on seven test-day records of first lactation in order to perform GWAS on the PCs explaining maximum variation rather than the test-day records themselves. The first two PCs cumulatively explained 40–78% of total variation in test-day records of traits included, indicating their potential to be used as GWAS traits to identify novel SNPs for milk yield, fat percentages, and SNF percentages.

It is important to understand the reasons for considering 96 samples taken in the study as optimum, as in the buffalo farming scenario in Asia, particularly in India, the maximum number of buffaloes maintained at any large organized herd ranges from 250 to 500, with 100–200 breedable buffaloes having complete phenotype information. The National Dairy Research Institute has the second-largest herd in the Indian Council of Agricultural Research (ICAR) with a well-managed herd of 250 breedable buffaloes. Furthermore, there is no buffalo sequencing consortium/project in operation (India). In such a case, one could afford this sample size to genotype with complete pedigree and with sufficient genetic diversity. However, the sample size is less for conducting GWAS and raises the question of whether the results are reliable enough. It was observed that the results obtained under the present study are encouraging and important, as one of the genes identified in the present study, i.e., GRIA3, was reported by Cole et al. (2011) in a sample of 1,654 US Holstein cows. There are occasions where a similar sample size (96) has been used to understand body morphology traits through GWAS (Rahmatalla et al., 2018).

To identify novel SNPs associated with various economic traits, sequencing through the ddRAD approach was performed. Ruffalo et al. (2012) highlighted that a mapping quality of 15–40 is of the highest nature, and theoretical accuracy corresponds to ∼100%, while depending on the aligner, the actual accuracy of mapping or base call may vary from 40 to 60%. Similar mapping quality averaging 30 for all samples was obtained through BWA aligner, which usually implies good overall base quality of reads and few mismatches in the best possible alignment.

An earlier bovine SNP chip was used to study the traits of buffalo. Using a bovine SNP chip (Illumina BovineSNP50 BeadChip) for a GWAS in buffalo population, Wu et al. (2013) identified seven SNPs that were significantly associated with milk yield. Later on, Affymetrix’s 90K SNP chip commercialization led to several GWAS on milk production traits in buffalo. de Camargo et al. (2015) in a GWAS on buffaloes reported that LOC100847171, BCL6, RTP2, SST, PTGS2, LOC100295047, and LOC101908004 are associated with milk production; KCTD8, LOC782855, LOC101904777, ESRRG, TRNAY-AUA, and GPATCH2 are associated with fat percentage; LOC101903483, SART3, ISCU, CMKLR1, WSCD2, MFNG, CARD10, and USP18 are associated with protein percentage. Liu et al. (2018) identified several candidate genes, namely, MFSD14A, SLC35A3, PALMD, RGS22, and VPS13B, for milk production in the Italian Mediterranean buffalo through GWAS. In another GWAS in buffalo, Iamartino et al. (2017) reported that genes regulating the D-glucose level in the blood affect milk production in buffalo.

In the present study, GWAS performed on the milk yield PCs explaining maximum variation, which revealed that three SNPs were significantly associated with PC₁ of TDMY at the 5% FDR level. However, six SNPs present on chromosomes 1, 6, and 9 at 182,059,836, 35,648,660, and 48,559,942 bp, respectively, were above the suggestive threshold of 10% FDR and are presented in Table 1 as the top six SNPs associated with PC₁ of TDMY. However, for the rest of the traits studied, no SNPs were detected to be significantly associated (or as rejections) at 5% or even 10% FDR.

Upon scanning ± 20 kb around that SNP on the X chromosome, the GRIA3 gene was found. Glutamate ionotropic receptor AMPA type subunit 3 (GRIA3) is reported to be very significantly associated with daughter pregnancy rate in US Holstein cows (Cole et al., 2011). In a study, a positive selection signal has been observed for the loci containing the CNIH3 gene. CNIH3 interacts with GRIA1, GRIA2, GRIA3, GRIA4, and GRIK1 belonging to a class of α-amino-3-hydroxy-5-methyl-4-isoxazolepropionic acid (AMPA) glutamate receptors, thus regulating trafficking of AMPA receptors. AMPA receptors are known to participate in luteinizing hormone (LH) secretion (Utsunomiya et al., 2013). The nearest genes to the top 10 SNPs from PC1 GWAS for milk yield were ZNF292, LOC112444602, TIGD2, GRIK2, and LRRC34. Similarly, for PC2 of milk yield, GWAS revealed that TCERG1, ANGPT4, ANKRD44, LOC112442949, and AKAP6 were the nearest genes from the top 10 SNPs. It could be observed that five SNPs present in the list of top SNPs from GWAS results of PC₁ and PC₂ of the test day’s milk yield were also present in the top 10 SNP list for 305 days’ milk yield (Table 1), highlighting the efficacy of PCA in predicting 305 days’ milk yield. Additionally, MYNN and ACTRT3 were found nearest to the top SNPs of GWAS for 305 days’ milk yield. Cai et al. (2020) have also reported the ANKRD44 gene to be associated with milk yield in Nordic Holstein cattle.

A genome-wide scan for novel genes for fat percentage outlined CCSER1, DAPK2, CAMTA1, ROR1, CCDC34, GSTA1, GSTA2, GRIK4, CACNG6, SH3BP5L, and ZNF672 as the nearest genes to the top SNPs obtained from GWAS. Kolbehdari et al. (2009) in a whole-genome study in Canadian Holstein cattle mapped the gene DAPK2 to several milk production-related QTLs. They also mapped DAPK2 to a milk fat percentage QTL which concords with the findings of our study. Genes other than DAPK2 were not reported earlier to be associated with milk fat percentage in buffalo; however, CAMTA1 may possibly have a role in fatty acid metabolism.

TICAM2, TXN, HNF4G, SYBU, LOC104974614, TRIQK, CFAP44, CDH18, LOC782977, and SEC63 were found near the top SNPs associated with PCs of TDSNF. We could not find any previous reports stating the role of these genes in regulating milk SNF percentage; however, a functional enrichment analysis study by Zhou et al. (2019) revealed HNF4G and SYBU as candidate genes that regulate milk fat synthesis, transport, and metabolism. All the genes involved with milk yield and its composition identified through the study were enriched for pathways in Cytoscape. A sub-network of the genes is depicted in Supplementary Figure 2. It can be observed that most of the genes identified belong to the glutamate receptor family, involved in the regulation of the neuronal system.

In Canadian Holstein cattle, Do et al. (2017) reported a significant genetic effect of MAN1C1, MAP3K5, HCN1, TSPAN9, MRPS30, TEX14, and CCL28 genes on lactation persistency. Deng et al. (2019) reported that the MAP3K5 gene regulates p38 MAPKs and Jun N-terminal kinases (JNKs) pathways involved in mammary gland development. In the present study, lactation persistency was calculated using four different methods, and it was observed that the top four SNPs identified by GWAS for persistency as estimated by the Johansson and Hansson (1940) method were in common with the GWAS results obtained for persistency as per the Ludwick and Petersen (1943) method. The genes present near the top SNPs were DENND3, PSTPIP1, ADAMTS20, PRODH2, NPHS1, RREL2, FAM19A1, CLIC5, GPR12, DEFB125, STPG2, CFAP206, RINT1, EFCAB10, TCERG1, and C18H16orf78. Upon network enrichment, it was observed that TCERG1, RINT1, CLIC5, and PSTPIP1 are co-expressed with several other genes in the breast mammary tissue of human, indicating their possible role in persistency of lactation in dairy animals.

GWAS for fertility traits were performed on breeding efficiency, age at sexual maturity, and postpartum breeding interval. Genes present near the top SNPs of breeding efficiency were APC, LOC112448352, NDFIP2, LOC107131404, and PDP1. Mohamed et al. (2019) reported that APC gene in the canonical WNT signaling pathway plays a critical role in the regulation of ovarian development. Mis-regulation of this key pathway in the adult ovary is associated with subfertility in mice. The APC gene has also been mapped to the QTL region for conception rate in Holstein cattle (Kolbehdari et al., 2009). Genes identified near the top SNPs for age at sexual maturity are PDE11A, SLAMF6, LOC104974658, GRID2, LOC104970778, EIF5A2, and RPL22L1. Coyral-Castel et al. (2018) mapped the SLAMF6 gene to the QTL region affecting female fertility located on the bovine chromosome three (QTL-F-Fert-BTA3). The GRID2 gene is however reported to be associated with growth traits in Simmental cattle and may have role in regulating age at sexual maturity (Braz et al., 2020). Christensen et al. (2005) have reported EIF5A2 as a candidate gene for infertility in human. Genes identified for postpartum breeding interval were ERICH2, MALSU1, IGF2BP3, TCF24, PPP1R42, CCDC93, GRID2, EDIL3, LOC112449367, LOC101904981, and LOC112447353. However, none of these genes was identified to be involved significantly in any pathway through the gProfiler.

Although we have already emphasized that this a preliminary GWAS on buffaloes covering a large set of economic traits (lactation, lactation persistency, and fertility), the results obtained in the present study are bias free, as indicated by the FDR. The findings of the present study in the Murrah population will encourage researchers to come forward for GWAS in buffalo in a holistic manner.

Conclusion

Optimum production and reproduction in buffaloes are long-standing questions in the Indian subcontinent. We present the analysis and identification of genomic regions that play role in shaping the selection and breeding decisions of Murrah buffalo for persistency of production and fertility. The dataset information presented in the paper is also submitted so that it can be used to compare and evaluate other breeds of buffalo based on the genomic information generated. The putative identified regions have a potential to improve the existing breeding decisions in this important dairy germplasm.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

• 1

  AndrewsS. (2010). FastQC: A Quality Control Tool for High Throughput Sequence Data. Available online at: http://www.bioinformatics.babraham.ac.uk/projects/fastqc(accessed August 21, 2020).

  • Google Scholar
  • View reference in article

• 2

  BarbosaA. M. (2020). FuzzySim: applying fuzzy logic to binary similarity indices in ecology.Methods Ecol. Evol.6853–858. 10.1111/2041-210X.12372

  • CrossRef
  • Google Scholar
  • View reference in article

• 3

  BenjaminiY.HochbergY. (1995). Controlling the false discovery rate: a practical and powerful approach to multiple testing.J. R. Stat. Soc. Series B Methodol.57289–300. 10.1111/j.2517-6161.1995.tb02031.x

  • CrossRef
  • Google Scholar
  • View reference in article

• 4

  BignardiA. B.El FaroL.RosaG. J. M.CardosoV. L.MachadoP. F.AlbuquerqueL. G. D. (2012). Principal components and factor analytic models for test-day milk yield in Brazilian Holstein cattle.J. Dairy Sci.952157–2164. 10.3168/jds.2011-4494

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 5

  BrazC. U.RowanT. N.SchnabelR. D.DeckerJ. E. (2020). Extensive genome-wide association analyses identify genotype-by-environment interactions of growth traits in Simmental cattle.Biorxiv. [Preprint].10.1101/2020.01.09.900902

  • CrossRef
  • Google Scholar
  • View reference in article

• 6

  BrianB. (2014). “BBMap: a fast, accurate, splice-aware aligner,” in Proceedings of the 9th Annual Genomics of Energy & Environment Meeting, (Berkeley, CA: Lawrence Berkeley National Lab).

  • Google Scholar
  • View reference in article

• 7

  BushW. S.MooreJ. H. (2012). Genome-wide association studies.PLoS Comput. Biol.8:e1002822. 10.1371/journal.pcbi.1002822

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 8

  CaiZ.DuszaM.GuldbrandtsenB.LundM. S.SahanaG. (2020). Distinguishing pleiotropy from linked QTL between milk production traits and mastitis resistance in Nordic Holstein cattle.Genet. Sel. Evol.52:19. 10.1186/s12711-020-00538-6

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 9

  ChangC. C.ChowC. C.TellierL. C.VattikutiS.PurcellS. M.LeeJ. J. (2015). Second-generation PLINK: rising to the challenge of larger and richer datasets.Gigascience4:7. 10.1186/s13742-015-0047-8

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 10

  ChristensenG. L.IvanovI. P.AtkinsJ. F.MielnikA.SchlegelP. N.CarrellD. T. (2005). Screening the SPO11 and EIF5A2 genes in a population of infertile men.Fertil. Steril.84758–760. 10.1016/j.fertnstert.2005.03.053

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 11

  ColeJ. B.WiggansG. R.MaL.SonstegardT. S.LawlorT. J.CrookerB. A.et al (2011). Genome-wide association analysis of thirty one production, health, reproduction and body conformation traits in contemporary US Holstein cows.BMC Genomics12:408. 10.1186/1471-2164-12-408

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 12

  Coyral-CastelS.RaméC.CognieJ.LecardonnelJ.MartheyS.EsquerréD.et al (2018). KIRREL is differentially expressed in adipose tissue from ‘fertil+’and ‘fertil-’cows: in vitro role in ovary?Reproduction155181–196. 10.1530/REP-17-0649

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 13

  da Costa BarrosC.de Abreu SantosD. J.Aspilcueta-BorquisR. R.De CamargoG. M. F.de Araújo NetoF. R.TonhatiH. (2018). Use of single-step genome-wide association studies for prospecting genomic regions related to milk production and milk quality of buffalo.J. Dairy Res.85402–406. 10.1017/S0022029918000766

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 14

  DaveyJ. W.HohenloheP. A.EtterP. D.BooneJ. Q.CatchenJ. M.BlaxterM. L. (2011). Genome-wide genetic marker discovery and genotyping using next-generation sequencing.Nat. Rev. Genet.12499–510. 10.1038/nrg3012

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 15

  de CamargoG. M. F.Aspilcueta-BorquisR. R.FortesM. R. S.Porto-NetoR.CardosoD. F.SantosD. J. A.et al (2015). Prospecting major genes in dairy buffaloes.BMC Genomics16:872. 10.1186/s12864-015-1986-2

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 16

  DengT.LiangA.LiangS.MaX.LuX.DuanA.et al (2019). Integrative analysis of transcriptome and GWAS data to identify the hub genes associated with milk yield trait in buffalo.Front. Genet.10:36. 10.3389/fgene.2019.00036

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 17

  DoD. N.BissonnetteN.LacasseP.MigliorF.SargolzaeiM.ZhaoX.et al (2017). Genome-wide association analysis and pathways enrichment for lactation persistency in Canadian Holstein cattle.J. Dairy Sci.1001955–1970. 10.3168/jds.2016-11910

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 18

  ElshireR. J.GlaubitzJ. C.SunQ.PolandJ. A.KawamotoK.BucklerE. S.et al (2011). A robust, simple genotyping-by-sequencing (GBS) approach for high diversity species.PLoS One6:e19379. 10.1371/journal.pone.0019379

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 19

  EwelsP.MagnussonM.LundinS.KällerM. (2016). MultiQC: summarize analysis results for multiple tools and samples in a single report.Bioinformatics323047–3048. 10.1093/bioinformatics/btw354

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 20

  FAO (2014). Food and Agriculture Organization the United Nations Statistics Division. Available online at: http://faostat3.fao.org/compare/E(accessed November 12, 2020).

  • Google Scholar
  • View reference in article

• 21

  GanguliN. C. (1981). Buffalo as a candidate for milk production.Int. Dairy Fed. Bull.13748–56.

  • Google Scholar
  • View reference in article

• 22

  García-AlcaldeF.OkonechnikovK.CarbonellJ.CruzL. M.GötzS.TarazonaS.et al (2012) Qualimap: evaluating next-generation sequencing alignment data.Bioinformatics282678–2679. 10.1093/bioinformatics/bts503

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 23

  GlickmanM. E.RaoS. R.SchultzM. R. (2014). False discovery rate control is a recommended alternative to Bonferroni-type adjustments in health studies.J. Clin. Epidemiol.67850–857. 10.1016/j.jclinepi.2014.03.012

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 24

  GordonI. (1996). Controlled Reproduction in Cattle and Buffaloes, Vol. 1. Wallingford: Centre for Agriculture and Bioscience International.

  • Google Scholar
  • View reference in article

• 25

  GriffithsM. W. (2010). Improving the Safety and Quality of Milk: Improving Buffalo Milk.Amsterdam: Elsevier. 10.1533/9781845699437

  • CrossRef
  • Google Scholar
  • View reference in article

• 26

  IamartinoD.NicolazziE. L.Van TassellC. P.ReecyJ. M.Fritz-WatersE. R.KoltesJ. E.et al (2017). Design and validation of a 90K SNP genotyping assay for the water buffalo (Bubalus bubalis).PLoS One12:e0185220. 10.1371/journal.pone.0185220

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 27

  JohanssonI.HanssonA. (1940). Causes of variation in milk and butterfat yield of dairy cows.Kungliga Lantbruksakademiens Handlingar79:127.

  • Google Scholar
  • View reference in article

• 28

  KhedkarC.KalyankarS.DeosarkarS. (2016). “Buffalo milk,” in Encyclopedia of Food and Health, edsCaballeroB.FinglasM.ToldráF. (New York, NY: Academic Press), 522–528. 10.1016/B978-0-12-384947-2.00093-3

  • CrossRef
  • Google Scholar
  • View reference in article

• 29

  KolbehdariD.WangZ.GrantJ. R.MurdochB.PrasadA.XiuZ.et al (2009). A whole genome scan to map QTL for milk production traits and somatic cell score in Canadian Holstein bulls.J. Animal Breed. Genet.126216–227. 10.1111/j.1439-0388.2008.00793.x

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 30

  LiH. (2011). A statistical framework for SNP calling, mutation discovery, association mapping and population genetical parameter estimation from sequencing data.Bioinformatics272987–2993. 10.1093/bioinformatics/btr509

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 31

  LiH.DurbinR. (2010). Fast and accurate long-read alignment with burrows–wheeler transform.Bioinformatics26589–595. 10.1093/bioinformatics/btp698

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 32

  LiuJ. J.LiangA. X.CampanileG.PlastowG.ZhangC.WangZ.et al (2018). Genome-wide association studies to identify quantitative trait loci affecting milk production traits in water buffalo.J. Dairy Sci.101433–444. 10.3168/jds.2017-13246

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 33

  LudwickT. M.PetersenW. E. (1943). A measure of persistency of lactation in dairy cattle.J. Dairy Sci.26439–445. 10.3168/jds.S0022-0302(43)92739-0

  • CrossRef
  • Google Scholar
  • View reference in article

• 34

  MacciottaN. P. P.VicarioD.Cappio-BorlinoA. (2005). Detection of different shapes of lactation curve for milk yield in dairy cattle by empirical mathematical models.J. Dairy Sci.881178–1191. 10.3168/jds.S0022-0302(05)72784-3

  • CrossRef
  • Google Scholar
  • View reference in article

• 35

  MahadevanP. (1951). The effect of environment and heredity on lactation. II. persistency of lactation.J. Agric. Sci.4189–93. 10.1017/S0021859600058573

  • CrossRef
  • Google Scholar
  • View reference in article

• 36

  MareesA. T.de KluiverH.StringerS.VorspanF.CurisE.Marie-ClaireC.et al (2018). A tutorial on conducting genome-wide association studies: quality control and statistical analysis.Int. J. Methods Psychiatr. Res.27:e1608. 10.1002/mpr.1608

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 37

  MohamedN. E.HayT.ReedK. R.SmalleyM. J.ClarkeA. R. (2019). APC2 is critical for ovarian WNT signalling control, fertility and tumour suppression.BMC Cancer19:677. 10.1186/s12885-019-5867-y

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 38

  PetersonB. K.WeberJ. N.KayE. H.FisherH. S.HoekstraH. E. (2012). Double digest RADseq: an inexpensive method for de novo SNP discovery and genotyping in model and non-model species.PloS One7:e37135. 10.1371/journal.pone.0037135

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 39

  RahmatallaS. A.ArendsD.ReissmannM.WimmersK.ReyerH.BrockmannG. A. (2018). Genome-wide association study of body morphological traits in Sudanese goats.Anim. Genet.49478–482. 10.1111/age.12686

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 40

  RainaV.NarangR.MalhotraP.KaurS.DubeyP. P.TekamS.et al (2016). Breeding efficiency of crossbred cattle and Murrah buffaloes at organized dairy farm.Indian J. Anim. Res.50867–871. 10.18805/ijar.v0iOF.6663

  • CrossRef
  • Google Scholar
  • View reference in article

• 41

  RuffaloM.KoyutürkM.RayS.LaFramboiseT. (2012). Accurate estimation of short read mapping quality for next-generation genome sequencing.Bioinformatics28i349–i355. 10.1093/bioinformatics/bts408

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 42

  SambrookJ.RussellD. W. (2006). Purification of nucleic acids by extraction with phenol: chloroform.Cold Spring Harb. Protoc.2006:4455. 10.1101/pdb.prot4455

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 43

  SchaefferL. R.JamrozikJ.Van DorpR.KeltonD. F.LazenbyD. W. (2000). Estimating daily yields of cows from different milking schemes.Livest. Prod. Sci.65219–227. 10.1016/S0301-6226(00)00153-6

  • CrossRef
  • Google Scholar
  • View reference in article

• 44

  ShannonP.MarkielA.OzierO.BaligaN. S.WangJ. T.RamageD.et al (2003). Cytoscape: a software environment for integrated models of biomolecular interaction networks.Genome Res.132498–2504. 10.1101/gr.1239303

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 45

  TaggarR. K.DasA. K.KumarD.MahajanV. (2012). Prediction of milk yield in Jersey cows using principal component analysis.Progr. Res.7272–274.

  • Google Scholar
  • View reference in article

• 46

  TomarN. S. (1965). A note on the method of working out breeding efficiency in Zebu cows and buffaloes.Indian Dairyman17389–390.

  • Google Scholar
  • View reference in article

• 47

  UtsunomiyaY. T.O’BrienA. M. P.SonstegardT. S.Van TassellC. P.do CarmoA. S.MészárosG.et al (2013). Detecting loci under recent positive selection in dairy and beef cattle by combining different genome-wide scan methods.PloS One8:e64280. 10.1371/journal.pone.0064280

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 48

  Van TassellC. P.SmithT. P.MatukumalliL. K.TaylorJ. F.SchnabelR. D.LawleyC. T.et al (2008). SNP discovery and allele frequency estimation by deep sequencing of reduced representation libraries.Nat. Methods5247–252. 10.1038/nmeth.1185

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 49

  VenturiniG. C.CardosoD. F.BaldiF.FreitasA. C.Aspilcueta-BorquisR. R.SantosD. J. A.et al (2014). Association between single-nucleotide polymorphisms and milk production traits in buffalo.Genet. Mol. Res.1310256–10268. 10.4238/2014.December.4.20

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 50

  WaraA. B.KumarA.SinghA.ArthikeyanA. K.DuttT.MishraB. P. (2019). Genome wide association study of test day’s and 305 days milk yield in crossbred cattle.Indian J. Anim. Sci.89861–865.

  • Google Scholar
  • View reference in article

• 51

  WarriachH. M.McGillD. M.BushR. D.WynnP. C.ChohanK. R. (2015). A review of recent developments in buffalo reproduction—a review.Asian-Australas. J. Anim. Sci.28:451. 10.5713/ajas.14.0259

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 52

  WoodP. D. P. (1967). Algebraic model of the lactation curve in cattle.Nature216164–165. 10.1038/216164a0

  • CrossRef
  • Google Scholar
  • View reference in article

• 53

  WuJ. J.SongL. J.WuF. J.LiangX. W.YangB. Z.WathesD. C.et al (2013). Investigation of transferability of BovineSNP50 BeadChip from cattle to water buffalo for genome wide association study.Mol. Biol. Rep.40743–750. 10.1007/s11033-012-1932-1

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 54

  ZhouC.ShenD.LiC.CaiW.LiuS.YinH.et al (2019). Comparative transcriptomic and proteomic analyses identify key genes associated with milk fat traits in Chinese Holstein cows.Front. Genet.10:672. 10.3389/fgene.2019.00672

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 55

  ZiminA. V.DelcherA. L.FloreaL.KelleyD. R.SchatzM. C.PuiuD.et al (2009). A whole-genome assembly of the domestic cow, Bos taurus.Genome Biol.10:R42. 10.1186/gb-2009-10-4-r42

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article