Comparative analysis of SNP candidates in disparate milk yielding river buffaloes using targeted sequencing

Sample collection and genomic DNA extraction

We recruited 12 water buffaloes, which were biological replicates from three breeds namely Banni, Mehsani and Jafrabadi, belonging to either high or low-milk yield category from the Gujarat state of India (Table 1). The gDNA was isolated from blood samples using a Qiagen DNeasy Blood and Tissue kit (Qiagen Corp., CA, USA) and the resultant DNA was quantified using Qubit® dsDNA BR Assay (Invitrogen Corp., CA, USA) and integrity was confirmed by agarose gels. Blood drawings were conducted in accordance with regulations and prior approval by the institutional animal ethics committee of Anand Agricultural University, Gujarat, India.

Probe design, target enrichment and exome sequencing

We obtained the intended targets (all coding exons, 3′UTR and 5′UTR exons) of cattle (Btau_4.6.1/bos Tau7) from RefGene tables of UCSC genome browser. Further, the custom probe design was performed by NimbleGen (Roche, Germany) which are compatible with Roche GS-FLX Titanium sequencer.

Rapid library for each sample was prepared from ∼1 μg of gDNA separately and multiplexed according to manufacturer’s protocol (Roche, Germany) using high quality DNA. Final libraries were used to set up hybridization reaction at isothermal temperature of 47 °C for 68–72 h in thermal cycler, with custom designed probes as per manufacturer’s protocol (NimbleGen). Captured DNA libraries were quantified spectrophotometrically, and evaluated electrophoretically with high sensitivity DNA assay on Agilent Bioanalyzer 2100 (Agilent). Finally, the libraries were sequenced on Roche 454 GS-FLX Titanium instrument according to manufacturer’s protocol (Roche, Germany). The raw sequence data will be deposited in a public repository.

Raw reads filtering, alignment and variant detection

Raw reads were filtered based on phred quality score ≥ 20 and length ≥ 50 bp using QTrim tool (Schmieder & Edwards, 2011). The resultant sequence reads were mapped against Bos taurus genome build 4.6.1 using ‘bwa-mem’ module of BWA v 0.7.5a (Li & Durbin, 2010). Potential PCR duplicates were removed using ‘MarkDuplicate’ module of Picard Tools. SNPs passing the criteria depth ≥ 5 and quality score ≥ 30 were identified using SAMtools, FreeBayes and VarScan tools, through pooled variant calling strategy (Li et al., 2009; Garrison & Marth, 2012; Koboldt et al., 2012). To deduce the shared (MultiCom) and specific SNPs from each tool, we used VCFtools package v0.1.11 and basic Linux commands (Danecek et al., 2011). The SNPs belonging to high-yield and low-yield were compared using VCFtools package v0.1.11 (Danecek et al., 2011). The specific SNPs (high or low yield) refers to the SNPs that were present in any of the sample in the given group which are not present or did not pass the SNP calling threshold in the other group.

Annotation and QTL dataset

The resulting sequence ontology (SO) and candidate gene annotation of SNPs specific for high and low yield groups was performed using SnpEff v 3.4 (Cingolani et al., 2012). The candidate gene enrichment analysis was performed in the GeneCodis tool and the enriched pathways and gene ontology (GO) terms were identified using Hypergeometric test followed by Benjamini-Hochberg’s correction (p-value < 0.01). Cattle QTL candidate genes belonging to milk-yield and milk-fat categories were obtained from the AnimalQTLdb online resource (Hu et al., 2013).

Mammary gland transcriptome analysis

From NCBI-GEO database, we obtained the whole gene expression data of 8 samples derived from mammary gland of Holstein Friesen cattle belonging to high or low milk yield (Accession: GSE33680). The gene expression was measured using Agilent-023647 v2 microarray, comprising of 45,220 oligo-nucleotide probes representing cattle genes. In the dataset, there were four samples in each for high-milk yield and low-milk yield groups. GEO2R implementation of LIMMA ebayes test was used to identify differentially regulated genes between high and low milk yield groups (p-value < 0.01).

Samples and overview of sequence data

In order to identify SNPs specific for high and low milk-yield in Indian buffalo breeds, we designed two balanced sample groups of healthy Indian buffaloes with prominent difference in milk yields. Animals with fertility issues, inconsistent diets or other known disorders were excluded from the study. The custom designed capture probes based on UCSC RefTable data of Bos taurus genome version 4.6.1/bosTau7 covered 125,679 exons, 14,084 3′UTRs and 16,574 5′UTRs. The sequencing experiment of 12 samples generated about 2.49 gigabases encompassed in ∼6.5 million sequence reads with average length of 380 bases (Table 2). We obtained ∼68% average on-target capture efficiency across samples and approximately 98% of filtered reads were mapped to Bos taurus 4.6.1 reference sequence (NCBI), with the average per sample depth ≥ 5× (Table S1).

Identification of SNPs specific for high and low milk-yield using MultiCom approach

We grouped the samples based on milk-yield (high and low milk-yield groups), and performed the SNP analysis separately for each group, using cattle as reference. Initially, we used three variant calling tools namely: SAMtools, FreeBayes and VarScan with the consistent thresholds of base depth ≥ 5 and quality score ≥ 30 for all the tools. In the high milk-yield group, SAMtools detected 1.097 million SNPs. FreeBayes and VarScan reported 1.101 million and 1.090 million SNPs respectively (Fig. 1A–1C). Even though high concordances were observed between the outcomes from the three tools, several thousands of SNPs (6.5–9.3% of output from each tool) were tool-specific or non-reproducible SNP calls (Fig. 1A). On the other hand, in the low milk-yield group, SAMtools identified 1.389 million SNPs, followed by 1.103 and 1.098 million SNPs respectively by VarScan and FreeBayes. In the low milk-yield group, the tool-specific SNP calls for SAMtools was the same as that of high-production dataset (6.5%), but a considerable reduction of tool–specific calls was observed for VarScan and FreeBayes (4.6 and 3.8%, respectively) outputs, even though the total outcome by these two tools were comparable with the high milk-yield group (Fig. 1A). We hypothesized that the non-reproducible SNPs were potential false positive calls. To test this, we compared the transition/transversion (ts/tv) ratios of SNPs which were tool-specific and shared, separately for the two groups. A highly consistent ts/tv range was observed for the shared SNPs, ranging from 2.61 to 2.65. On contrary, the tool-specific or non-reproducible SNPs showed a dynamic range of lower ts/tv ratios, ranging from 1.60 to 2.38 (Fig. 1D). Strikingly, the difference in the ts/tv ratios were highly significant (p-value < 0.005) (Fig. 1D). Next, using SO annotation by snpEff tool, we compared the regions in the reference genome in which the SNPs are located. Notably, it was observed that the inter-genic region SNPs were almost doubled in case of tool-specific SNP calls, consistently in three tools and the two groups (Fig. 1E). Therefore, we considered only those SNPs which were reproduced by at least two tools (hereon called MultiCom approach). There were a total of 1.203 million SNPs in high milk-yield group and 1.315 million SNPs for low milk-yield group (Fig. 1A–1C). Interestingly, in almost all cases, SNP selection through MultiCom approach facilitated an increased discovery rate (up to 20%) compared to the outcome from each tool. Finally, we compared MultiCom outcomes in high milk-yield and low milk-yield groups, and have found that 255,741 SNPs (21.3%) were specific for high yield group, and 367,674 SNPs (28.0%) were specific for low yield group. On the other hand, the 947,646 SNPs were common to the two groups which may be generic SNPs in buffalo, as the reference genome was cattle.

Candidate gene prediction and comparison with milk-yield and milk-fat QTL dataset

Using snpEff prediction candidate genes were determined for high and low milk-yield specific SNPs (Table S2). The snpEff annotation mapped the 255,741 high yield specific SNPs in 7,212 genes, and the 367,674 low yield specific SNPs in 8,284 coding genes. Even though the SNPs were specific for the two groups, we have found that the 5,037 genes were common to high and low milk-yield groups. On the other hand, there were 2,175 genes specific for high yield and 3,247 for low milk-yield. In order to determine whether these genes were present among the reported milk-yield and milk-fat related genes in QTL region, we obtained the list of candidate genes in milk-yield and milk-fat categories of cattle resources of the animal genome data repository. Remarkably, of 74 milk-yield candidates, 51 genes were present in our results, among which 26 were common to high and low milk-yield group and 17 genes specific for low yield specific genes, and 8 were high specific genes (Fig. 2A). On the other hand, of 91 milk-fat candidates, 30 genes were common to high and low-milk yield group, 24 and 15 candidate genes were specific for low and high milk-yield groups respectively (Fig. 2B). On the whole, about 70% of the QTL candidate genes in respective categories belonging to cattle genome annotation were present among the genes the genes identified by our study in buffalo breeds.

Enrichment analysis of yield specific candidate genes

Next, we performed the GO enrichment analysis of the candidate genes, specific for high and low milk-yield groups with the threshold of FDR corrected p-value < 0.01. In the GO molecular functions category, eight GO terms were enriched in high milk-yield group. On the other hand, 21 terms were enriched in low-milk yield group, among which four terms were common to GO terms belonging to high milk-yield group (Fig. 3). The common terms were metal ion binding, nucleotide binding, zinc ion binding and cytokine activity. Interestingly, several GO terms in the low milk-yield category which were related to transcriptional regulation, such as sequence-specific DNA binding transcription factor activity (GO:0003700), sequence-specific DNA binding (GO:0043565) and transcription regulatory region DNA binding (GO:0044212) (Fig. 3A). Collectively, 140 genes were found to be contributing to the statistical enrichment of these three GO terms, which were specific for low-milk yield group (Fig. 3C).

In order to obtain the clues about the functional consequences of the candidate genes specific for high and low-yield specific groups, we performed the KEGG pathway enrichment analysis using Genecodis tool. There were 15 pathways enriched for the high yield specific group and 32 pathways for the low-yield specific group (Table S3). However, we have noticed that several of the pathways were related to human diseases. Further, after excluding human disease specific pathways, there were 8 pathways present in high-yield group and 19 pathways for the low milk yield group (Table 3). The most significant pathways in the low milk-yield group were Oxidative phosphorylation, Toll-like receptor signalling, cytokine-cytokine receptor interaction etc. On the other hand, the most significant among the high-yield group were Jak-STAT signalling pathway, Wnt signalling pathway, ErbB signalling pathway etc. Surprisingly, Toll-like receptor signalling and MAPK signalling pathway were enriched in both high and low milk-yield groups by distinct set of genes falling in the same pathway (Figs. S1–S3).

Comparison of buffalo milk-yield specific candidate genes with milk-yield DEG (Differentially Expressed Genes) in cattle

Next, we hypothesized that the SNPs in the candidate gene may exert changes in the gene expression that could have an effect on milk yield. In order to verify this, we obtained the mammary cell transcriptome dataset belonging to high and low milk-yielding cattle (see Materials and Methods). LIMMA analysis between the high milk-yield and low milk-yield group detected 1,056 annotated differentially expressed genes, in which 579 gene were up-regulated and 477 genes were down-regulated among the low milk-yield samples (Fig. 4; Table S4). Quite interestingly, it was found that 582 out of 1,056 DEG were buffalo SNP candidate genes in low yield specific, high yield specific or common to both groups (Fig. 4). Among these 582 overlapping DEG, 196 genes (118 up-regulated and 78 down-regulated) belonged to low yield specific group, 139 genes (68 up-regulated and 71 down regulated) were high yield specific group. On the other hand, 247 DEG (147 up-regulated and 100 down-regulated in low milk yield) were common to high and low milk-yield groups. However, we didn’t observe any trend in the direction of fold change with respect to milk-yield group. Next, in order to assess whether the overlap between the DEG and the SNP candidates are not due to random chance, we performed Hypergeometric tests on each category viz. low yield specific, high yield specific and common. Interestingly, in all categories p-value < 0.0001. Thus, considering that SNP candidate genes accounted more than half of the DEG in a cross-species comparison, we propose that buffalo SNP candidates potentially influence downstream gene expression changes to a large extend.