Exploring genetic polymorphism in innate immune genes in Indian cattle (Bos indicus) and buffalo (Bubalus bubalis) using next generation sequencing technology.

Activation of innate immunity initiates various cascades of reactions that largely contribute to defense against physical, microbial or chemical damage, prompt for damage repair and removal of causative organisms as well as restoration of tissue homeostasis. Genetic polymorphism in innate immune genes plays prominent role in disease resistance capabilities in various breeds of cattle and buffalo. Here we studied single nucleotide variations (SNP/SNV) and haplotype structure in innate immune genes viz CHGA, CHGB, CHGC, NRAMP1, NRAMP2, DEFB1, BNBD4, BNBD5, TAP and LAP in Gir cattle and Murrah buffalo. Targeted sequencing of exonic regions of these genes was performed by Ion Torrent PGM sequencing platform. The sequence reads obtained corresponding to coding regions of these genes were mapped to reference genome of cattle BosTau7 by BWA program using genome analysis tool kit (GATK). Further variant analysis by Unified Genotyper revealed 54 and 224 SNPs in Gir and Murrah respectively and also 32 SNVs was identified. Among these SNPs 43, 36, 11,32,81,21 and 22 variations were in CHGA, CHGB, CHGC, NRAMP1, NRAMP2, DEFB1 and TAP genes respectively. Among these identified 278 SNPs, 24 were found to be reported in the dbSNP database. Variant analysis was followed by structure formation of haplotypes based on multiple SNPs using SAS software revealed a large number of haplotypes. The SNP discovery in innate immune genes in cattle and buffalo breeds of India would advance our understanding of role of these genes in determining the disease resistance/susceptibility in Indian breeds. The identified SNPs and haplotype data would also provide a wealth of sequence information for conservation studies, selective breeding and designing future strategies for identifying disease associations involving samples from distinct populations.

Introduction

Innate immunity is considered to be an outcome of natural selection and plays an important role against natural pathogens in livestock. Hot and humid atmosphere is very favorable to disease causing microorganisms and parasite proliferation. In spite of such natural habitat, it has long been observed and established that Indian cattle and buffaloes are less susceptible to bacterial, viral and parasitic diseases. Livestock breeds of Indian origin are well adapted to tropical environments because of their ability to thrive under extreme nutritional stress, resistance to diseases, heat tolerance potential and sturdiness (Pal et al., 2011). The Gir animals mainly inhabiting the tropical regions are renowned for their stress tolerance ability and resistance to various tropical diseases and bullocks of this breed are well known for their heavy load dragging capacity on all kinds of soil (Gaur et al., 2003). Murrah breed of Bubalus bubalis is one of the most diffused breeds in the world with high milk and meat productivity and stress tolerance capacity (Thiruvenkadan et al., 2013). This might be attributed to a very strong innate immune system of these animals. Genes mainly involved in bovine innate immunity involve beta defensins viz DEFB1, BNBD4, BNBD5, TAP (Tracheal Antimicrobial Peptide), LAP (Lingual Antimicrobial Peptide); cathelicidins, toll like receptors, chromatogranins mainly chromatogranin A, chromatogrann B, chromagranin C; and NRAMP (Natural Resistance Associated Macrophage Pathogen) formally called SLC11 (Solute Carrier Family 11) including NRAMP1, NRAMP2 etc. (Hancock and Diamond, 2000, Zanetti, 2005).

Beta defensins are a major family of cationic antibacterial peptides of 29–45 amino acid residues having broad spectrum of antimicrobial activity and hence as effector molecules of innate immunity in many livestock species including cattle, buffalo, and sheep (Mahoney et al., 1995, Zhao et al., 1999). The cationic antimicrobial peptides expressed at the site of infection attach to the anionic microbial membrane lead to permeabilization of lipid bilayer and subsequent release of cellular contents and hence demonstrating microbicidal activity. Potent antibacterial activity of neutrophil beta defensins against Staphylococcus aureus, Escherichia coli, Mycobacterium paratuberculosis and Pseudomonas aeruginosa was observed (Bals, 2000, Kaiser and Diamond, 2000, Selsted et al., 1993). Another major proteins involved are natural resistance associated macrophage proteins i.e., NRAMP. NRAMP are integral membrane proteins that mediate antimicrobial activity of macrophages (both cytocidal and cytostatic) against intracellular parasites during early infection (Vidal et al., 1993). NRAMP are divalent cation transporters which are located in the phagolysosome membrane of macrophages promoting macrophage mediated intracellular pathogen killing and hence play a critical role in bovine innate immunity. In cattle, they were studied especially for their association with resistance against Brucella abortus infection (Feng et al., 1996), Mycobacterium bovis, Salmonella typhimurium, and Leishmania donovani in the mouse, Mus musculus (Ables et al., 2002). NRAMP2 are Fe2 + transporters having much broader range of tissue expression and mutations in NRAMP2 lead iron deficiency thus conferring their role in iron metabolism (Gruenheid et al., 1999). The chromogranins are proteins of 27–100 kDa having at least 10% acidic (glutamic/aspartic acid) residues and multiple single and dibasic amino acid residues found as principal soluble components of the secretory granules in neuroendocrine cells. They are secreted from these cells in a physiologically regulated manner having autocrine and paracrine activities (Feldman and Eiden, 2003). The chromatogranin family includes chromatogranin A (CHGA, Secretory protein 1/SP-1), chromatogranin B (CHGB, Secretogranin I/SCG I), and chromatogranin C (CHGC, Secretogranin II/SCG II) (Natori and Huttner, 1994).

Genetic polymorphism in some innate immune genes in Indian livestock has been associated with resistance against natural pathogens. Apart from polymorphisms, haplotyping data provides information about the recombination of chromosome which is informative for locating disease causing mutations by linkage methods using linkage disequilibrium. Haplotypes are a set of SNPs or a combination of alleles along the same chromosome that are inherited together as a unit (Crawford and Nickerson, 2005).

The present study specifically focused on detecting SNPs among 10 candidate immune genes and its haplotype structure in Indian Gir cattle (Bos indicus) and Murrah buffalo (B. bubalis). Targeted amplification and sequencing of exonic regions of these candidate genes using Ion Torrent PGM sequencing platform is used to explore polymorphism and haplotypes in indigenous breed of cattle and buffalo.

Materials and methods

Animals, sample collection and total DNA isolation

Gir breed of cattle (B. indicus) and Murrah breed of buffalo (B. bubalis) each involving ten individuals were selected. Blood was collected aseptically using 2.7% ethylenediamine tetraacetic acid as an anti-coagulant agent from healthy old aged cattle and buffalo. DNA was isolated using John's non-idet method and was quality-checked on 0.8% agarose gel and quantified using a microvolume spectrophotometer (NanoDrop 1000; Thermo Fisher Scientific Inc., USA).

PCR amplification

A total of 10 genes were selected (CHGA, CHGB, CHGC, NRAMP1, NRAMP2, DEFB1, BNBD4, BNBD5, LAP and TAP). Gene sequence of each gene was downloaded from http://www.ncbi.nlm.nih.gov/ and primers targeting exonic regions were designed using Clone manager (version 9) software. A total of 96 primer pairs were designed for 10 genes having expected amplicon size ranging from 380 to 420 bp. Primer sequences along with expected amplicon sizes are given in supplementary Table 2. PCR-reactions were carried out with cycling conditions as, initial denaturation at 98 °C for 2 min; denaturation at 98 °C for 10 s, annealing at 60 °C for 10 s, and extension at 72 °C for 25 s for 32 cycles followed by final extension at 72 °C for 2 min. Reactions were performed in 15 μl total volume and contained 90 ng of bovine genomic DNA, 10 pmol of gene-specific primers and other reagents were added as recommended for use with Emerald PCR master mix (Takara, US) in MicroAmp® 96-Well Plate. Amplified fragments were visualized on 2% agarose gel. After confirmation of the expected size amplification on an agarose gel, the PCR products were pooled sample wise followed by visualization of pooled amplicons on 2% Agarose gel. Purification of pooled amplicons was performed using standard Gel cut method using the QIAquick gel extraction kit (Qiagen, CA). DNA concentration was estimated with the Qubit 2.0 instrument applying the Qubit dsDNA HS Assay (Life Technologies, Invitrogen division, Darmstadt, Germany).

Table 2

Data obtained after quality filtration by Prinseq lite.

Sample	Input Sequences	Good Sequences	Bad Sequences
Gir1	40,926	32,658 (79.80%)	8268 (20.20%)
Gir2	179,697	133,640 (74.37%)	46,057 (25.63%)
Gir3	46,772	37,963 (81.17%)	8809 (18.83%)
Gir4	231,909	176,262 (76.00%)	55,647 (24.00%)
Gir5	159,749	123,624 (77.39%)	36,125 (22.61%)
Gir6	312,477	217,764 (69.69%)	94,713 (30.31%)
Gir7	213,969	162,828 (76.10%)	51,141 (23.90%)
Gir8	154,613	115,661 (74.81%)	38,952 (25.19%)
Gir9	148,807	115,542 (77.65%)	33,265 (22.35%)
Gir10	182,220	139,652 (76.64%)	42,568 (23.36%)
Murrah1	187,635	147,105 (78.40%)	40,530 (21.60%)
Murrah 2	194,958	147,893 (75.86%)	47,065 (24.14%)
Murrah 3	95,211	75,764 (79.57%)	19,447 (20.43%)
Murrah 4	196,234	152,713 (77.82%)	43,521 (22.18%)
Murrah 5	183,165	140,540 (76.73%)	42,625 (23.27%)
Murrah 6	168,935	128,166 (75.87%)	40,769 (24.13%)
Murrah 7	139,400	107,093 (76.82%)	32,307 (23.18%)
Murrah 8	161,186	116,401 (72.22%)	44,785 (27.78%)
Murrah 9	188,480	141,113 (74.87%)	47,367 (25.13%)
Murrah 10	200,916	146,907 (73.12%)	54,009 (26.88%)

Sequencing

Amplicon library preparation

Libraries were generated using 150 ng of the amplicon DNA and an Ion Xpress™ Plus Fragment Library Kit comprising the Ion Shear chemistry following manufacturer's instruction. From the concentration and the average size of each amplicon library, the amount of DNA fragments per micro liter was calculated and libraries were diluted to 2.8 × 108 DNA molecules per micro liter prior to clonal amplification. Libraries were quantified on Qubit 2.0 instrument and concentrations for each library of breed were given in Supplementary Table 3 and quality checked on Agilent 2100 bioanalyzer (Agilent Biotechnologies, CA, USA).

Emulsion PCR and sequencing

The emulsion PCR was carried out applying the Ion XPress Template kit V2.0 (Life Technologies) as described in the appropriate user guide provided by the manufacturer. Sequencing of the amplicon libraries was carried out on the Ion Torrent Personal Genome Machine (PGM) system using the Ion Sequencing 400 kit (Life Technologies) following the corresponding protocol. The 316v2 chip was used and chip check and calibration was performed.

Alignment and variant calling

Quality filtering was performed with Prinseq lite software tool with the parameters set as min quality mean 20, min length 60 and by removing duplicates. Burrows–Wheeler aligner (BWA; version 0.5.9) was used to align the sequencing reads to the Bos taurus reference genome (BosTau7). Alignments were converted from sequence alignment map (SAM) format to sorted, indexed binary alignment map (BAM) files. GATK software tools (version 2.8; http://www.broadinstitute.org) were used for genotype calling with recommended parameters. Genotypes were called by the GATK Unified Genotyper tool, and variants were filtered by depth 60.

Haplotype analysis

Haplotype analysis was performed with the SAS software version 9.2 (SAS Institute Inc., Nashville, TN, USA). The procedure uses Expectation Maximization (EM) algorithm to predict the maximum likelihood estimates of the haplotype frequencies assuming Hardy–Weinberg equilibrium. Linkage disequilibrium (LD) option in PROC HALPOTYPE was specified and the haplotype frequency estimates from the alternative hypothesis were reported. The initial haplotype frequencies used in the EM algorithm are identical to the frequencies that appear in the H0 frequency column. The frequencies in the H1 frequency column are those calculated from the final iteration of the EM algorithm, and these frequencies' standard errors and confidence limits are included as well. Haplotypes with lower frequency were not included in the subsequent analysis. Overall workflow is given in Fig. 1.

Fig. 1

Overall work-flow.

Results

Approximately 27.863 kb of DNA spanning 10 innate immunity genes i.e., CHGA, CHGB, CHGC, NRAMP1, NRAMP2, DEFB1, BNBD4, BNBD5, LAP and TAP were sequenced in 10 individuals of Gir breed of cattle and 10 individuals of Murrah breed of buffalo.

Upon amplicon sequencing using Ion Torrent PGM platform, 895.60 Mb data were generated with 3,413,789 of total number of reads and 301 bp mean read length and other pieces of information are given in Supplementary Table 4. Sample wise total number of bases, total number of reads and mean read length are given in Table 1.

Table 1

Total number of bases and reads obtained sample wise upon amplicon sequencing.

Sample name	Bases	≥ Q20 bases	Total reads	Mean read length
Gir 1	55,311,671	48,597,157	187,635	294 bp
Gir 2	55,295,488	48,838,583	194,958	283 bp
Gir 3	27,537,462	24,091,168	95,211	289 bp
Gir 4	58,969,828	51,418,252	196,234	300 bp
Gir 5	57,021,245	49,813,918	183,165	311 bp
Gir 6	51,567,567	44,733,052	168,935	305 bp
Gir 7	42,447,537	37,139,277	139,400	304 bp
Gir 8	49,046,863	42,640,608	161,186	304 bp
Gir 9	57,919,134	50,332,489	188,480	307 bp
Gir 10	60,988,685	52,735,483	200,916	303 bp
Murrah 1	11,449,315	9,956,188	40,926	279 bp
Murrah 2	55,690,520	48,310,747	179,697	309 bp
Murrah 3	13,702,949	11,939,237	46,772	292 bp
Murrah 4	72,317,420	63,121,844	231,909	311 bp
Murrah 5	48,112,507	41,875,243	159,749	301 bp
Murrah 6	95,013,059	82,605,555	312,477	304 bp
Murrah 7	65,426,759	56,851,610	213,969	305 bp
Murrah 8	44,155,685	38,426,589	154,613	285 bp
Murrah 9	43,655,953	37,850,066	148,807	293 bp
Murrah 10	55,000,392	48,136,977	182,220	301 bp

Variant calling

Quality filtration using Prinseq resulted in to ~ 22% data reduction. Sample-wise total number of input sequences and obtained good sequences are given in below Table 2.

Upon variant calling using GATK software tools, a total of 54 and 224 SNPs were observed in Gir and Murrah breed respectively and 32 SNVs were shared by both breeds. From the total variations, 10 and 21 SNPs were found to be reported in dbSNP in Gir and Murrah respectively. Total SNPs in Gir and Murrah breeds, shared variations (SNVs) in both breeds and SNPs reported in dbSNP were given in Supplementary Table 5. In Gir breed, 13 SNPs were identified as non-synonymous. Apart from these, 14 synonymous SNPs were obtained along with 26 SNPs in UTR regions. Similarly in Murrah, 50 SNPs were non-synonymous with 61 synonymous and 121 SNPs in UTR region.

For chromagranins, in CHGA gene, number of non-synonymous variations observed i.e.,8 was higher followed by 4 synonymous and 3 and 2 in UTR 3′ and UTR 5′ regions respectively in Gir. Similarly for the same gene in Murrah, 14 non-synonymous, 11 synonymous and 12 SNPs in UTR regions were observed. While for CHGB gene, only one SNP in Gir breed i.e., non-synonymous and in Murrah 35 SNPs was revealed with higher number of 13 synonymous SNPs and also one stop gained SNP which causes start codon to be mutated in non-start codon. In CHGC gene, two variations were revealed in Gir breed, one was in UTR region and another was synonymous in Gir. In Murrah, for the same gene, 8 and 3 SNPs in UTR region and synonymous were identified respectively. For NRAMP1 gene, in Gir breed, two SNPs (one synonymous and one non-synonymous) and in Murrah breed, higher number of synonymous variations was identified. In NRAMP2 gene and genes of beta defensin family, maximum variations were in UTR region in both breeds. Also start gained variations were identified in DEFB1 gene. These start-gain mutations cause variation in the UTR 5′ region resulting in a three base sequence that can be a start codon. Gene wise distribution of SNPs in both breeds was given in Fig. 2.

Fig. 2

Distributions of SNPs.

Haplotyping

Haplotypes for Gir and Murrah breeds were given in Supplementary Table 6, Supplementary Table 7 respectively. In Gir and Murrah breeds, a total of 28 and 85 haplotypes were observed with significant haplotype frequency respectively.

Discussion

Most interesting components of innate immunity are a repertoire of receptors that recognize features common to many pathogens. SNPs occurring within the molecular components of innate immunity might be used as markers for determining response of an individual to infectious agent or therapies (Boekholdt et al., 2003). The importance of genomic variation within innate immune genes/molecules is supported by a number of observations (Collins et al., 1998). First, variation affecting key innate immune genes is heritable and provides evolutionary benefits to livestock populations. Second, these genes are important in initiating a host defense mechanism and offer clues to the adaptive antigen-specific response. Last, in given host–environment interface functions, genetic variations in innate immune genes are likely to alter downstream responses extending to the memory function of the adaptive response.

Variations might alter receptor function, their existence is most important, especially when seen in the context of innate immunity's evolutionary heritage. During millions of decades, a finely tuned variable host response developed polymorphisms enabling functional diversity among individuals within the same species (Janeway and Medzhitov, 2002, Schuster and Nelson, 2000). An association between polymorphisms in some innate immune genes and disease resistance has been reported (Mucha et al., 2009, Vidal et al., 1993). An association between genetic polymorphisms in TLR4 gene and NRAMP1 gene with susceptibility to brucellosis was observed in B. bubalis and bovine, respectively (Martinez et al., 2008). Thus, innate immune gene polymorphisms can be used to predict resistance against microbial infection. Present study aimed at identification of SNPs and haplotype structure with their frequencies for Gir breed of cattle and Murrah breed of buffalo. Sample size is an important parameter when estimating haplotype frequencies but the key determinant for association is LD across population. The overall pattern of LD for cattle and buffalo breed in our study showed that majority of SNPs in Murrah breeds were in significant linkage disequilibrium. Some level of haplotype sharing was observed in buffalo breed and cattle breed which represents retention of conserved ancestral variation.

Chromagranins are secretory proteins which serve as precursors for biologically active peptides involved in innate immunity of bovine (Bera et al., 2007). As observed 20 and 83 SNPs were in chromagranins in Gir and Murrah respectively and from these, a total of 12 SNPs were observed to be reported in dbSNP. None of those were reported to be studied for its association with disease susceptibility or resistance.

As observed in our study, polymorphisms in NRAMP1 were obtained both in Gir and Muraah. One SNP (rs109453173) in NRAMP1 gene was checked earlier for its association with paratuberculosis in Holstein-Friesian cattle but no significant association of paratuberculosis was observed with that SNP (Lazarus et al., 2002). Similar studies were reported in water buffalo in which single amino acid substitution from Gly to Asp in the murine NRAMP1 at position 169 was reported to be associated with susceptibility to Mycobacteria, Salmonella and Leishmania (Malo et al., 1993). NRAMP1 was also studied for its association with infection by Mycobacterium avium subspecies paratuberculosis in 14 bovine breeds. A significant association was detected between 2 SNPs (c.1067C > G and the A allele of c.1157–91A > T) out of 57 SNPs and susceptibility to infection by MAP in Holstein-Friesian cattle (Ruiz-Larranaga et al., 2010). Such SNPs were not observed in our study. Similarly mutation in NRAMP2 gene was reported in anemic mice models. One mutation (G185R) in two animal models mk mouse and Belgrade rat showed iron deficiency. They displayed deficiency in intestinal iron uptake in mk mouse and microcytic anemia in Belgrade rat (Gruenheid et al., 1999). In our case, numbers of non-synonymous SNPs were obtained but not such non-synonymous SNP was observed in NRAMP2 gene.

In our study, two from selected five beta defensin genes showed significant number of polymorphisms in UTR 3′ and UTR 5′ regions upon sequencing. Six novel non-synonymous SNPs were detected which were spread in TAP gene. Earlier studies reported change in expression level of beta defensin family genes during infection in bovine. Bovine BNBD4 was sequenced earlier which showed ~ 70% homology with cattle BNBD4 and on the contrary more similarity with cattle BNBD5 (Bera et al., 2007). Ten SNPs have been observed earlier in BNBD4 gene and among these; two SNPs generated new endonuclease digestion sites and therefore could easily be genotyped using the restriction fragment length polymorphism (RFLP) technique (Bagnicka et al., 2010). While in our study, no such SNP was obtained in BNBD4 gene. In our study, 10 and 14 SNPs were identified in DEFB1 gene in Gir and Murrah breed respectively which were similar to previous SNPs reported in DEFB1 gene. Variant analysis in beta defensin genes in five human ethnic populations has earlier revealed 15 SNPs in DEFB1 and DEFB2 genes harboring 3 SNPs in coding regions, 3 in UTR 5′ and 6 in UTR 3′region in DEFB1 gene. Similarly in our study we also observed total SNVs in both the UTR regions with a much higher frequency (19 in UTR 5′ region).

Extensive variation observed in our study of cattle and buffalo innate immune genes is possibly an outcome of their inhabiting area. Densities of SNPs observed in our study were significantly higher. This corroborates the Parham's proposition that innate immune genes are continuously evolving due to natural selection pressure.

Conclusion

A large number of genetic variations were observed in 7 out of 10 candidate genes and our analysis of polymorphisms within these genes and haplotype structure for the same will enable studies to access the relationships between variations within these genes and their susceptibility to various livestock diseases especially bovine diseases. In addition, variations in specific gene will contribute to the design of markers related to specific health related traits in bovine. Relating SNPs/SNVs of innate immune genes to clinically defined diseases will move the discovery process forward. This advance will be aided by new scientific understandings, inexpensive molecular assays, large clinical studies, and the need for the provision of efficient health care.

This is the first high throughput sequencing based approach to identify variations in B. indicus (Gir cattle) and B. bubalis (Murrah buffalo) in Indian livestock breeds. The identified genetic variations and haplotype in this study will be useful resource for future studies to understand genetic variations in innate immunity genes and resultant phenotypic variation in defensive effectiveness of Indian cattle and buffalo. Amplicon generation using PCR and NGS based sequencing proves to be an efficient method for variation and haplotype identification in small scale experiments for limited budget and sample. Given the large numbers of SNPs/SNVs generated that meet these minimal requirements, more stringent parameters can be applied in picking SNP sets for different applications. A very large number of SNPs generated among innate immune genes in our study will enrich the dbSNP database of NCBI (http://www.ncbi.nlm. nih.gov/projects/SNP/) and will be useful resource for array designing.

The following are the supplementary data related to this article.

Supplementary Tables

List of innate immunity genes with their size and location on genome.

Supplementary Table 2

Sequence of primers used for amplification of each innate immunity gene and expected product size.

Supplementary Table 5

SNP found in Gir and Murrah with position, dbSNP ID and their effect.

Supplementary Table 6

Haplotype frequencies in innate immunity genes in Gir cattle.

Supplementary Table 7

Haplotype frequencies in innate immunity genes in Murrah buffalo.