In Silico Analysis of Genome Wide Non-Synonymous Single Nucleotide Polymorphisms in Indigenous Cattle Breeds of Pakistan

Livestock have traditionally played important part in development of human societies for a long time. According to a rough estimate, livestock supports the livelihoods of approximately 0.6 billion farmers in underdeveloped nations (Thornton et al., 2006). Domestication of cattle has a major role in development of a country’s agriculture economy. They contribute more than milk, meat, fiber skin and fuel productions. Based on their physical appearances nearly 800 cattle breeds are divided into two categories, Bos taurus (taurine) and Bos indicus (indicine) (Felius et al., 2011). Beside many other reported differences, the appearance of a hump over the shoulders is one of the major phenotypic characteristics used to distinguish indicine breeds form taurine breeds (Zeder et al., 2006). Indicine cattle grow in tropical and subtropical climates such as Africa, Southeast Asia, Brazil, northern Australia, southern China, and sections of the United States thanks to its physiological characteristics. Taurine, on the other hand, is typically found in more developed European countries, North America, and Australia, where it has a greater metabolic rate and dietary requirements. There are collectively 35 recognized cattle breeds in Pakistan and India (Felius et al., 2014) which make considerable contribution to the agriculture economy. After goats, cows are the second most common animal species in Pakistan with an estimated 51.5 million animals (Finance Division, 2021), accounting for 3.2% of world cattle population (Singh et al., 2014). According to environmental conditions, most cattle breeds are created as drought breeds like Bhagnari, while Red Sindhi and Sahiwal are milk breeds and Tharparkar, Achai, Gabrali, and Cholistani are dual-purpose cattle breeds (Felius et al., 2014).

The process of natural and human driven selection has changed the cattle genotypes leading to diverse genetic and phenotypic profiles, and adaptation to temperature and tropical environments (MacHugh et al., 1997; Porto-Neto et al., 2013). Since the completion of the taurine cattle genome (Elsik et al., 2009), thousand bull genome project (Hayes and Daetwyler, 2019), the worldwide bovine genome sequencing, and the HapMap project (Gibbs et al., 2009), a considerable number of genetic variants including single nucleotide polymorphisms (SNPs) and InDels (Insertions/ Deletions) that have been recorded in a publicly available database (Sherry et al., 2001; Eck et al., 2009; Elsik et al., 2009; Gibbs et al., 2009; Iqbal et al., 2019).

Due to the advancements and cheap cost of next-generation sequencing (NGS), it is now feasible to simultaneously generate sequence data and estimate SNP allele frequencies in a range of reference populations (Kumar et al., 2012). Whole genome SNP genotyping analysis detected genomic regions targeted by selection, which, for example, contain immune and environmental adaptation related genes. This has opened new doors for studying the genetic variation and processes which help in adaptation to different biogeographic regions (Iso-Touru et al., 2016; Weldenegodguad et al., 2019). SNPs are the most prevalent type of genetic variation that are stable, have a low mutation rate, are inherited in a mendelian manner, and are relatively inexpensive to genotype (Snelling et al., 2005). As a result, they are routinely utilized as biomarkers in genetic research. Advances in bioinformatics and statistical tools have also improved our understanding of demographic evolution, the possible role of genomic structural variations, adaptation during domestication and selection, and the biological functions of these genomic variations in livestock breed (Gutenkunst et al., 2009; Li and Durbin, 2011; Jacob et al., 2020).

Regardless of these recent breakthroughs, there is a severe gap in whole-genome investigations of Asian cow breeds and European bovine breeds that are commonly employed for meat and dairy production (Kawahara-Miki et al., 2011; Choi et al., 2014). In this study, we analyzed the whole genome sequence data for four different purpose cattle breeds of Pakistan. Our main goal was to identify and analyze the missense variants resulting from nsSNP and predict their effect on resulting amino acid sequence, their structure, functionality, and stability. Missense variations not only modify the tertiary structure of proteins, but also lead to formation of deleterious phenotypes. Results of this study will be useful in furthering research into the genetic pathways underpinning features of interest in Indian cattle.

The whole genome sequences of Pakistani Bos indicus cattle samples included in this study were recently made available in public domain (Iqbal et al., 2019). Sequence data is available in China National GenBank Sequence Archive under project accession number CNP0000189. Total genomic DNA was extracted form whole blood (10 mL) samples and were sequenced using BGISEQ-500. For this study, we choose four different purpose breeds including Sahiwal (CNS0014619), Bhagnari (CNS0014606), Cholistani (CNS0014604), and Red Sindhi (CNS0014620).

Quality check of sequence reads was performed using fast QC (Andrews, 2010) tool to identify low quality reads and adapter sequences. Raw sequences were first subjected to filtering process to identify and remove low quality reads and adapters using trimmomatic software (Bolger et al., 2014). Burrows-Wheeler Alignment (BWA) (Li and Durbin, 2009) algorithm was used with default parameters to create index files for reference genome ARS-UCD1.2 Btau5.0.1Y accessed from the 1000 Bull Genome project resources (Hayes and Daetwyler, 2019). Short reads of each selected cattle sample were mapped to reference genome using BWA with option “bwa-mem”. SAM files generated as result of alignment were converted to their binary equivalent (BAM) files and sorted using SAM tools (Li et al., 2009). PCR duplicates were removed using Picard tools “Mark Duplicate” command line utility from aligned reads (Picard tools by broad institute https://broadinstitute.github.io/picard/). Base Quality Score Recalibration (BQSR) was performed on mapped reads to resolve high or low base quality scores estimation. Final BAM files were used for downstream variant calling.

Single nucleotide polymorphisms (SNPs) along with insertion and deletions mutations (InDels) were identified using “Haplotype Caller” tool of genome analysis toolkit (GATK) which performs local de-novo assembly of haplotypes in genomic variation regions (McKenna et al., 2010). GATK was used according to guidelines of 1000 Bull Genome Project available at http://www.1000bullgenomes.com/doco/1000bullsGATK3.8pipelineSpecifications_Run8_Revision 20191101.docx (Fernandes Júnior et al., 2020). GATK “Select Variants” mode was used to separate SNPs and InDels into separate files. All variants including SNPs and InDels were discovered as difference from the reference genome ARS-UCD1.2 sequence. To remove the false positive calls, variants were filtered with GATK “Variant FIltration” argument using hard-filters with the following exclusion criteria: (1) quality by depth - QD<2.0; (2) Fisher Strand test - FS>60.0; (3) root mean square of the mapping quality score - MQ < 40.0; (4) and SOR > 9.0. SNPs and INDELS of autosomal chromosomes + X were left after filtration process.

Genetic variants were annotated using SnpEff program with annotation database ARS-UCD1.2.99 (Cingolani et al., 2012). SnpEff software assigns each SNPs and InDels to a functional class and provide several fields of information describing the affected transcripts and proteins, if applicable. SNPs and InDels were assigned to different functional classes including exonic, intronic, intergenic, splice site acceptor, splice site donor, splice site region, downstream, upstream, UTR 3 prime and UTR 5 prime region. InDels were also allocated to the in-frame deletion and insertion, disruptive in-frame deletion, and insertion functional classes whereas functional classes stop lost, stop retained, initiator codon were assigned to SNPs only. Variants were categorized into missense, nonsense, and silent variants. Missense or non-synonymous single base pair substitution result in amino acid substitution which affects the proteins phenotypically and functionally.

SIFT tool (Vaser et al., 2016) was used to predicts the genes with tolerated and deleterious SNPs from the missense variants. This program identifies the tolerated and deleterious SNPs to predict the impact of amino acid substitution on phenotypic and functional alterations in protein molecules. SIFT program to identifies tolerance score (TI) for each in range from 0.0 to 1.0. The SIFT score ≤ 0.05 labels the non-synonymous variant as deleterious to protein function (Ng and Henikoff, 2003; Sim et al., 2012). Gene enrichment analysis was performed to identify over-represented biological processes using genes with deleterious SNPs. This gene lists comprising was submitted to the ClueGO (Bindea et al., 2009), a Cytoscape (Shannon et al., 2003) plug-in that combines Gene Ontology and KEGG pathways to produce a well-organized GO/ pathway annotation network, allowing researchers to see if these genes are linked to key pathways.

Genomic sequences data of 4 different cattle breeds of Pakistan was analyzed which included dairy breeds Sahiwal and Red Sindhi, dual purpose breed Cholistani and drought tolerant breed Bhagnari. Individual sample sequence reads ranged from 754,986,078 to 1,605,207,524. Mapping of each sample to cattle reference genome ARS-UCD1.2 using BWA software yielded 83% to 99% alignment or sample reads (Table I). SAM alignment files were sorted and then converted to their binary format (BAM) files. Flow diagram depicts the bioinformatics analysis performed on the new generation sequence data (Fig. 1).

Alignment BAM files (including Bhagnari, Sahiwal, Cholistani and Red Sindhi) were pooled together to perform variant calling. Total ~25.46 M SNPs and ~3.56 M InDels were identified in the 29 autosomes and X chromosome against bovine reference genome ARS-UCD1.2 (Fig. 2 and Supplementary Table S1). Among the 3,563,505 InDels, 1,87,986 were deletions. Number of SNPs present per individual samples were 13.19, 18.96, 15.63, and 15.21 million in Bhagnari, Cholistani, Sahiwal, and Red Sindhi, respectively. The transition to transversion (Ts/Tv) ratio was pretty much similar across the samples, 2.27 in Bhagnari, 2.3 in Cholistani, 2.28 in Sahiwal and 2.3 in Red-Sindhi (Table II). This is indicative of high quality of our SNPs identified. Number of SNPs and Indels depend on the length of chromosomes. Out of total SNPs and InDels, only 104152 (0.41%) SNPs and 409 (0.01%) InDels matched to previously reported variations in dbSNP build 150.

Variant sets including SNPs and InDels were functionally annotated them to their attribute and genes using SnpEff tool. In our dataset of total SNPs identified, 16,750,695 (41.53%) were discovered in between genes. 1,785,969 (4.43%) SNPs were in thousand base pair (bp) upstream region of genes and 1,826,857 (4.53%) were in thousand base-pair downstream region of genes (Table II). Due to distinct isoforms or overlapping genes, multiple functional effects were identified when compared the total number of SNPs and InDels used for variant. 58,970,854 transitions and 25,313,817 transversions were detected with Ts/Tv ratio of 2.3296 in SNPs. SnpEff identified 2,909 high impact (disruptive), 122,494 moderate and 239,647 low impact SNPs in all sequenced samples (Supplementary Table SII). Functional Annotation of InDels predicted the presence of 2,331,013 (40.78%) and 2,779,680 (48.63%) InDels in intergenic and intronic regions, respectively. Moreover, 285,210 (5.00%) and 273, 913 (4.80%) indels were located within thousand base pair upstream and downstream genic regions (Table III). 5,155 high impact (disruptive), 2,359 moderate and 5,730 low impact InDels were identified in all sequenced samples (Supplementary Table SII). SnpEff categorized the identified SNPs as missense, nonsense and silent according to their functional class (Table V). Missense variations (nsSNPs) can possibly cause amino acid substitutions which phenotypically and functionally effect important traits of cattle. For further downstream analysis, missense variants were carried forward.

Phenotypic effects of the missense variants were identified by SIFT tool. SIFT uses sequence homology-based approach to characterize the effect of amino acid substitution on protein function. SIFT analysis predicted 154 deleterious missense variants overlapping 134 genes in coding sequence regions with a score between 0.00-0.04. Moreover, 325 missense variants in coding sequence regions were predicted as tolerated (Table VI and Supplementary Table S3). Even though a small number of deleterious SNPs were identified, all these variants overlapped protein coding genes. Gene list comprising deleterious SNPs was submitted in ClueGO to identify the functional gene ontology (GO) biological processes (BP). Gene ontology (GO) network is represented in Figure 3. Enrichment analysis revealed presence of these genes in several biological processes. Highlighted biological processes in gene interaction network included lipid transport (GO:0010876), transmembrane transport (GO:0055085), response to organic cyclic compounds (GO:0014070), response to endoplasmic reticulum stress (GO:0034976), negative regulation of transport (GO:0051051), detection of chemical stimulus (GO:0009593), skin development (GO:0043588), endomembrane system organization (GO:0010256), cytoskeleton organization (GO:0007010), response to growth factors (GO:0070848) and DNA conformation change (GO:0071103). List of all biological pathways and genes related is present in Supplementary Table S3.

In silico approach to characterize genetic variations can assist us in predicting the consequences of mutations and explain their affecting role in biological mechanisms. In this study, we investigated whole genome sequence data of four different purpose cattle breeds in Pakistan (Cholistani, Sahiwal, Bhagnari and Red Sindhi) to identify genetic variations which lead toward deleterious effect on protein structure and functions. This analysis led to the discovery of 25,469,157 SNPs and 3,563,505 InDels in these breeds.

Numerous SNPs, InDels and genes were found that have not been annotated yet. We annotated each variant’s possible functional role, allowing us to find numerous functionally significant candidate variants. Functional annotation revealed 122,943 missense (non-synonymous) SNPs in all samples 2,909 and 5,155 high impact (disruptive) SNPs and InDels in were also predicted by annotation. Analysis of missense variations revealed several genes with deleterious variation, involved in different biological and cellular functions. Go enrichment analysis revealed that genes harboring deleterious variations are significantly enriched in important biological processes, such as metabolic processes, development, transport, and homeostasis processes. Because this whole-genome sequencing study used just one animal of each breed, more research is needed to understand the exact dynamics of each gene–trait combination. The findings of this study can further be used as significant resource for further research on genomic characteristics to find variations in economically important traits and to develop precise genomic tools for cow breeding.

The authors are highly grateful to Director NIGAB and NARC for providing all the necessary facilities.

There is supplementary material associated with this article. Access the material online at: https://dx.doi.org/10.17582/journal.pjz/20211221081250

The authors have declared no conflict of interest.